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

龙格库塔法（Runge-Kutta methods）是用于非线性常微分方程的解的紧张的一类隐式或显式迭代法。在工程范畴应用广泛，可用于求解复杂体系的运动方程等题目。
这里采取matlab步伐编写代码实现龙格库塔法对于二元二阶常微分方程的求解。
例

{x¨+x˙+x+2y¨=−cost,x¨+y¨=−sint−cost,x(0)=0,x˙(0)=1,y(0)=1,y˙(0)=0,\left\{\begin{array}{lr}\ddot{x}+\dot{x}+x+2\ddot{y}=-\rm cos\it t, &  \\\ddot{x}+\ddot{y}=-\rm sin\it t-\rm cos\it t, \\x(0)=0, \\\dot{x}(0)=1, \\y(0)=1, \\\dot{y}(0)=0, \\\end{array}\right.⎩⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎧​x¨+x˙+x+2y¨​=−cost,x¨+y¨​=−sint−cost,x(0)=0,x˙(0)=1,y(0)=1,y˙​(0)=0,​
求当t=(0,40)t=(0,40)t=(0,40)的时间xxx，yyy的值。
该方程有准确解
{x=sinty=cost\left\{\begin{array}{lr}x=\rm sin\it t \\y=\rm cos\it t \\\end{array}\right.{x=sinty=cost​
下面我们通过四阶龙格库塔法举行求解：
起首设{u1=x,u2=x˙w1=y,w2=y˙\left\{\begin{array}{lr}u1=x,u2=\dot{x}\\w1=y,w2=\dot{y}\end{array}\right.{u1=x,u2=x˙w1=y,w2=y˙​​
则上式变动为：
{f1(t)=u1˙=u2f2(t)=u2˙=x¨=x˙+x−2sint−costf3(t)=w1˙=w2f4(t)=w2˙=y¨=−x˙−x+sint\left\{\begin{array}{lr}f1(t)=\dot{u1}=u2 \\f2(t)=\dot{u2}=\ddot{x}=\dot{x}+x-2\rm sin\it t-\rm cos\it t \\f3(t)=\dot{w1}=w2\\f4(t)=\dot{w2}=\ddot{y}=-\dot{x}-x+\rm sin\it t\end{array}\right.⎩⎪⎪⎨⎪⎪⎧​f1(t)=u1˙=u2f2(t)=u2˙=x¨=x˙+x−2sint−costf3(t)=w1˙=w2f4(t)=w2˙=y¨​=−x˙−x+sint​
则有四个子函数

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

复制代码

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

复制代码

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

复制代码

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

复制代码

龙格库塔迭代函数

1. function [u1,u2,u3,w1,w2,w3] = RK4_2variable(u1 , u2 , w1 , w2 , h , a , b , P1 , P2 , T)
3. x = a:h:b;
5. for i = 1:length(x)-1
7. k11 = f1(x(i) , u1(i) , u2(i) , w1(i) , w2(i));
8. k21 = f2(x(i) , u1(i) , u2(i) , w1(i) , w2(i) , P1 , P2 , T);
9. L11 = f3(x(i) , u1(i) , u2(i) , w1(i) , w2(i));
10. L21 = f4(x(i) , u1(i) , u2(i) , w1(i) , w2(i) , P1 , P2 , T);
12. 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);
13. 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);
14. 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);
15. 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);
17. 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);
18. 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);
19. 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);
20. 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);
22. k14 = f1(x(i)+h , u1(i)+h*k13 , u2(i)+h*k23 , w1(i)+h*L13 , w2(i)+h*L23);
23. k24 = f2(x(i)+h , u1(i)+h*k13 , u2(i)+h*k23 , w1(i)+h*L13 , w2(i)+h*L23 , P1 , P2 , T);
24. L14 = f3(x(i)+h , u1(i)+h*k13 , u2(i)+h*k23 , w1(i)+h*L13 , w2(i)+h*L23);
25. L24 = f4(x(i)+h , u1(i)+h*k13 , u2(i)+h*k23 , w1(i)+h*L13 , w2(i)+h*L23 , P1 , P2 , T);
27. u1(i+1) = u1(i) + h/6 * (k11 + 2*k12 + 2*k13 + k14);
28. u2(i+1) = u2(i) + h/6 * (k21 + 2*k22 + 2*k23 + k24);
29. w1(i+1) = w1(i) + h/6 * (L11 + 2*L12 + 2*L13 + L14);
30. w2(i+1) = w2(i) + h/6 * (L21 + 2*L22 + 2*L23 + L24);
32. u3(i) = k11;
33. w3(i) = L11;
34. end
35. end

复制代码

主函数

1. % main func
2. clear;clc;
3. u1(1) = 0;
4. u2(1) = 1;
5. w1(1) = 1;
6. w2(1) = 0;
7. h=0.01;
8. a = 0;b=20;
9. t = a:h:b;
10. [shu2,suoyin2]=unique(t2);
11. t2_=shu2; PU_=PU(suoyin2);
12. P1=-cos(t);
13. P2=-sin(t)-cos(t);
14. [u1,u2,u3,w1,w2,w3] = RK4_2variable(u1,u2,w1,w2,h,a,b,P1,P2,t);

复制代码

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