MATLAB龙格库塔法

⾮刚性常微分⽅程的数值解法通常会⽤四阶龙格库塔算法，其matlab函数对应ode45。

对于dy/dx = f(x,y)，y(0)=y0。

其四阶龙格库塔公式如下：

对于通常计算，四阶已经够⽤，四阶以上函数f(x,y)计算⼯作量⼤⼤增加⽽精度提⾼较慢。

下⾯以龙格库塔法解洛伦兹⽅程为例：

matlab代码如下：

main.m：

1 clear all;

2 close all;

3 clc;

5 %系统龙格库塔法

6 [t,h] = ode45(@test_fun,[040],[1240]);

7 plot3(h(:,1),h(:,2),h(:,3));

8 grid on;

10 %⾃定义龙格库塔法

11 [t1,h1]=runge_kutta(@test_fun,[1240],0.01,0,40);

12 figure;

13 plot3(h1(1,:),h1(2,:),h1(3,:),'r')

14 grid on;

runge_kutta.m（函数参考⽹络）：

1 %参数表顺序依次是微分⽅程组的函数名称，初始值向量，步长，时间起点，时间终点（参数形式参考了ode45函数）

2 function [x,y]=runge_kutta(ufunc,y0,h,a,b)

记住我3 n=floor((b-a)/h); %步数

4 x(1)=a; %时间起点

5 y(:,1)=y0; %赋初值，可以是向量，但是要注意维数

6for i=1:n %龙格库塔⽅法进⾏数值求解

7 x(i+1)=x(i)+h;

8 k1=ufunc(x(i),y(:,i));

9 k2=ufunc(x(i)+h/2,y(:,i)+h*k1/2);

10 k3=ufunc(x(i)+h/2,y(:,i)+h*k2/2);

11 k4=ufunc(x(i)+h,y(:,i)+h*k3);

12 y(:,i+1)=y(:,i)+h*(k1+2*k2+2*k3+k4)/6;

13 end

test_fun（洛伦兹⽅程）：

1 %构造微分⽅程

2 function dy=test_fun(t,y)

3 a = 16;

4 b = 4;

5 c = 45;

7 dy=[a*(y(2)-y(1));

8 c*y(1)-y(1)*y(3)-y(2);

9 y(1)*y(2)-b*y(3)];

得到很经典的洛伦兹吸引⼦，结果如下：