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;
4
5 %系统龙格库塔法
6 [t,h] = ode45(@test_fun,[040],[1240]);
7 plot3(h(:,1),h(:,2),h(:,3));
8 grid on;
9
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;
6
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)];
记住我
得到很经典的洛伦兹吸引⼦，结果如下：