用四阶RungeKutta方程解二阶常微分方程，并计算船舶在规则波中的横摇角（附Matlab代码）

前几天接到师姐分派的任务，让我求解一艘船模的横摇角的时间历程曲线，为后期的减摇控制准备。

1 首先冷静分析一下，原方程如下：

我们要求解的就是theta角和时间t之间的关系曲线，这是一道典型的二阶常微分方程的求解，用四阶龙格库塔方程可以求解。

2 古典龙格库塔算法公式：

高等数值计算课本（清华大学出版社，186页）

但是古典龙格库塔方法解决的是一阶常微分方程，也就是类似于这样的方程组

但是我们现在面对的是

查了一下书，没有教二阶常微分该怎么解，于是自己网上找答案：

3 网上的一道例题提供了解决办法

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

%横摇角时间历程曲线计算主程序
global a b c d h t y z D h_
a = 18.4;
b = 0.3;
c = 214.88;
d = 0;
h = 0.1;
y = 0;
z = 0;
D