四阶龙格库塔法求解二元二阶常微分方程

龙格库塔法（Runge-Kutta methods）是用于非线性常微分方程的解的重要的一类隐式或显式迭代法。在工程领域应用广泛，可用于求解复杂系统的运动方程等问题。
这里采用matlab程序编写代码实现龙格库塔法对于二元二阶常微分方程的求解。

例

$\{\begin{array}{lr}\ddot{x}+\dot{x}+x+2\ddot{y}=-\mathrm{c}\mathrm{o}\mathrm{s}t,& \\ \ddot{x}+\ddot{y}=-\mathrm{s}\mathrm{i}\mathrm{n}t-\mathrm{c}\mathrm{o}\mathrm{s}t,& \\ x(0)=0,& \\ \dot{x}(0)=1,& \\ y(0)=1,& \\ \dot{y}(0)=0,& \end{array}$
求当$t=(0,40)$的时候$x$，$y$的值。
该方程有精确解
$\{\begin{array}{l}x=\mathrm{s}\mathrm{i}\mathrm{n}t\\ y=\mathrm{c}\mathrm{o}\mathrm{s}t\end{array}$
下面我们通过四阶龙格库塔法进行求解：
首先设$\{\begin{array}{l}u1=x,u2=\dot{x}\\ w1=y,w2=\dot{y}\end{array}$
则上式变换为：
$\{\begin{array}{l}f1(t)=\dot{u1}=u2\\ f2(t)=\dot{u2}=\ddot{x}=\dot{x}+x-2\mathrm{s}\mathrm{i}\mathrm{n}t-\mathrm{c}\mathrm{o}\mathrm{s}t\\ f3(t)=\dot{w1}=w2\\ f4(t)=\dot{w2}=\ddot{y}=-\dot{x}-x+\mathrm{s}\mathrm{i}\mathrm{n}t\end{array}$
则有四个子函数

function output = f1(x,u1,u2,w1,w2)
output = u2;

function output = f2(x,u1,u2,w1,w2)
output = u2+u1-2*sin(x)-cos(x);

function output = f3(x,u1,u2,w1,w2)
output = w2;

function output = f4(x,u1,u2,w1,w2)
output = -u2-u1+*sin(x);

龙格库塔迭代函数

function [u1,u2,u3,w1,w2,w3] = RK4_2variable(u1 , u2 , w1 , w2 , h , a , b , P1 , P2 , T)

x = a:h:b;

for i = 1:length(x)-1

k11 = f1(x(i) , u1(i) , u2(i) , w1(i) , w2(i));
k21 = f2(x(i) , u1(i) , u2(i) , w1(i) , w2(i) , P1 , P2 , T);
L11 = f3(x(i) , u1(i) , u2(i) , w1(i) , w2(i));
L21 = f4(x(i) , u1(i) , u2(i) , w1(i) , w2(i) , P1 , P2 , T);

k12 = f1(x(i)+h/2 , u1(i)+h*k11/2 , u2(i)+h*k21/2 , w1(i)+h*L11/2, w2(i)+h*L21/2);
k22 = f2(x(i)+h/2 , u1(i)+h*k11/2 , u2(i)+h*k21/2 , w1(i)+h*L11/2, w2(i)+h*L21/2 , P1 , P2 , T);
L12 = f3(x(i)+h/2 , u1(i)+h*k11/2 , u2(i)+h*k21/2 , w1(i)+h*L11/2, w2(i)+h*L21/2);
L22 = f4(x(i)+h/2 , u1(i)+h*k11/2 , u2(i)+h*k21/2 , w1(i)+h*L11/2, w2(i)+h*L21/2 , P1 , P2 , T);

k13 = f1(x(i)+h/2 , u1(i)+h*k12/2 , u2(i)+h*k22/2 , w1(i)+h*L12/2 , w2(i)+h*L22/2);
k23 = f2(x(i)+h/2 , u1(i)+h*k12/2 , u2(i)+h*k22/2 , w1(i)+h*L12/2 , w2(i)+h*L22/2 , P1 , P2 , T);
L13 = f3(x(i)+h/2 , u1(i)+h*k12/2 , u2(i)+h*k22/2 , w1(i)+h*L12/2 , w2(i)+h*L22/2);
L23 = f4(x(i)+h/2 , u1(i)+h*k12/2 , u2(i)+h*k22/2 , w1(i)+h*L12/2 , w2(i)+h*L22/2 , P1 , P2 , T);

k14 = f1(x(i)+h , u1(i)+h*k13 , u2(i)+h*k23 , w1(i)+h*L13 , w2(i)+h*L23);
k24 = f2(x(i)+h , u1(i)+h*k13 , u2(i)+h*k23 , w1(i)+h*L13 , w2(i)+h*L23 , P1 , P2 , T);
L14 = f3(x(i)+h , u1(i)+h*k13 , u2(i)+h*k23 , w1(i)+h*L13 , w2(i)+h*L23);
L24 = f4(x(i)+h , u1(i)+h*k13 , u2(i)+h*k23 , w1(i)+h*L13 , w2(i)+h*L23 , P1 , P2 , T);

u1(i+1) = u1(i) + h/6 * (k11 + 2*k12 + 2*k13 + k14);
u2(i+1) = u2(i) + h/6 * (k21 + 2*k22 + 2*k23 + k24);
w1(i+1) = w1(i) + h/6 * (L11 + 2*L12 + 2*L13 + L14);
w2(i+1) = w2(i) + h/6 * (L21 + 2*L22 + 2*L23 + L24);

u3(i) = k11;
w3(i) = L11;
end
end

主函数

% main func
clear;clc;
u1(1) = 0;
u2(1) = 1;
w1(1) = 1;
w2(1) = 0;
h=0.01;
a = 0;b=20;
t = a:h:b;
[shu2,suoyin2]=unique(t2);
t2_=shu2; PU_=PU(suoyin2);
P1=-cos(t); 
P2=-sin(t)-cos(t); 
[u1,u2,u3,w1,w2,w3] = RK4_2variable(u1,u2,w1,w2,h,a,b,P1,P2,t);

运行主函数即可得到近似结果