四阶龙格库塔方程（Rungekutta）解二阶常微分方程组并计算船舶在迎浪下的纵摇埀荡耦合运动方程-附Matlab代码

今年年初的时候给师姐做了DDPG算法的船舶减横摇控制算法，师姐还有想法要让我把纵摇-埀荡两个自由度的减摇也做出来，这个任务归我了。实际上不管是多少个自由度的减摇，其实都需要解运动方程，当初做单自由度横摇的时候就是用四阶龙格库塔求解二阶微分方程，具体代码见：

https://blog.csdn.net/Mezikov/article/details/103898387

0 写在前面

代码内容在第二节，第一节是运动方程各种参数求解。

1 冷静分析一下方程组

这个方程组是典型的二阶常微分方程组，一堆的水动力参数，但都是可以求出来的已知量（求解方法看下图，Z和M分别是方程1等号右边项和方程2等号右边项），未知量是埀荡量ZG和纵摇角theta及其一阶二阶导。

2 现在可以开始编程了，我没有用Matlab自带的Ode45函数，我是自己编写的四阶Rungekutta方程，代码如下：

首先是纵摇运动方程 pitch equation

function w2 = pitch_equation_test3(y1,z1,y2,z2)
% global B d delta V v g L rho A_w C_w C_p C_b W_e x_b GM_l b_zz fi
fi = pi/3;       % 波浪遭遇角
B = 1;           % width of boat
d = 0.341;         % draft
delta = 1.76*10^3;  % displacement
V = 1.76;           % displacement volume
g = 9.8;
L = 6;             % length
lamda = 6;         % 波长
rho = 1*10^3;
A_w = 5.448;        % area of water plane
A_m = 0.346;        % '中横剖面面积'
C_w = A_w / (L * B);  % coefficient of water plane
C_p = V/(L*A_m);
C_b = delta/(1000*L*B*d);
v = 4.44;              %船速
T = 0.8*sqrt(lamda);
W_0 = 2*pi/T  ;      % 波浪圆频率
W_e = W_0*(1+v*W_0*cos(fi)/g) ;   % 波浪遭遇频率'%%%%%%%%%%%%%%%%%
x_b = -0.23 ;     % 'center of float' 这个很可能是错的，因为我自己变得,按照大尺度船的计算是0.23，但是根据参考文献，这里应该是负的
GM_l = (L^2*(5.55*C_w+1)^3/3450)/(C_b*d);
M = 1.6*10^3;      % weight
H_0 = B/(2*d);
a_zz = 0.8*H_0*C_w*M;
b_zz = (5.4*(C_w/C_p)*(H_0)^(0.5)-4.7)*delta/((g*L)^0.5);
c_zz = rho*g*A_w;
a_z0 = (v/W_e^2)*b_zz-a_zz*x_b;
b_z0 = -b_zz*x_b+v*a_zz;     % 0其实就是theta
c_z0 = -rho*g*A_w*x_b;
I_00 = 0.07*C_w*M*L^2;
a_00 = 0.83*H_0*C_p^2*(0.25*L)^2*(delta/g)^0.5;
b_00 = 0.08*H_0*delta*L^2/(g*L)^0.5;
c_00 = rho*g*(V*GM_l);
a_0z = -(a_zz*x_b);
b_0z = -(b_zz*x_b);
c_0z = -rho*g*A_w*x_b-v*b_zz;
fi = pi/3;
a = 0.5*0.17*lamda^(3/4);
k = 2*pi/lamda;
n1 = (I_00+a_00)-(a_0z*a_z0)/(M+a_zz);
n2 = rho*g*a*k*exp(-k*d/2)*B*d*sin(k*B*L*sin(fi)/2)/(k*B*L*sin(fi)/2);
n3 = (I_00+a_00)*(M+a_zz)-a_0z*a_z0;
  
w2 =  (b_00*(M+a_zz)-a_0z*b_z0)/n3*z2+...
      (c_00*(M+a_zz)-a_0z*c_z0)/n3*y2+...
      (b_0z*(M+a_zz)-a_0z*b_zz)/n3*z1+...
      (c_0z*(M+a_zz)-a_0z*c_zz)/n3*y1;
             
end

然后是埀荡运动方程 heave equation

function w1 = heave_equation_test3(y1,z1,y2,z2)
% global B d delta V v g L rho A_w C_w C_p C_b W_e x_b GM_l b_zz fi
fi = pi/3;       % 波浪遭遇角2*pi/3
B = 1;           % width of boat
d = 0.341;         % draft
delta = 1.76*10^