微分方程组的龙格库塔公式求解matlab版

下面给出经典四阶龙格-库塔格式解常微分方程组的 matlab 通用程序： %marunge4s.m %用途：4 阶经典龙格库塔格式解常微分方程组 y'=f(x,y),y(x0)=y0 %格式：[x,y]=marunge4s(dyfun,xspan,y0,h) %dyfun 为向量函数 f(x,y),xspan 为求解区间[x0,xn], %y0 为初值向量，h 为步长，x 返回节点，y 返回数值解向量 function [x,y]=marunge4s(dyfun,xspan,y0,h) x=xspan(1):h:xspan(2); y=zeros(length(y0),length(x)); y(:,1)=y0(:); for n=1:(length(x)-1) k1=feval(dyfun,x(n),y(:,n)); k2=feval(dyfun,x(n)+h/2,y(:,n)+h/2*k1); k3=feval(dyfun,x(n)+h/2,y(:,n)+h/2*k2); k4=feval(dyfun,x(n+1),y(:,n)+h*k3); y(:,n+1)=y(:,n)+(h/6).*(k1+2*k2+3*k3+k4); end
再编写运行函数： clear all; close all; clc; [t,y]=marunge4s(@mafun,[0 500],[0 1 2],0.005); plot(y(1,:),y(3,:),'r'); 得到如下图所示的结果：
60
50
40
30
20
10
0 -25
-20
-15
-10
-5
0
5
如下为例题： 例 1：取 h=0.02，利用程序 marunge4s.m 求刚性微分方程组 y' 0.01 y 99.99 z , y 0 2 , z' 100 z , z 0 1 的数值解，其解析解为： y e 0.01x e 100 x , z e 100 x 解：首先编写 M 函数 dyfun.m %dyfun.m function f=dyfun(t,y) f(1)=-0.01*y(1)-99.99*y(2); f(2)=-100*y(2); f=f(:); 然后编写一个执行函数： close all; clear all; clc; [x,y]=marunge4s(@dyfun,[0 500],[2 1],0.002); plot(x,y); axis([-50 500 -0.5 2]); text(120,0.4,'y(x)'); text(70,0.1,'z(x)'); title('数值解'); figure; y1=exp(-0.01*x)+exp(-100*x);%解析解，用来做对比的。 z1=exp(-100*x); plot(x,y1,'r'); hold on; plot(x,z1,'g'); text(120,0.4,'y(x)'); text(70,0.1,'z(x)'); axis([-50 500 -0.5 2]); title('解析解'); 可以看出数值解和解析解的运算结果误差很小：
10
15
20
25
2. 高阶微分方程组 对于高阶微分方程，总是可以化成方程组的形式。例如，二阶方程 y" g x , y , y' y x0 y0 , y' x0 y'0 总是可以化为一阶方程组： y' z z' g x , y , z y x y , z x y' z 0 0 0 0 0 因此不需要对于高阶方程给出计算公式。 例 3：求二阶方程 y" 2 y 3 , 1 x 1.5 y 1 y' 1 1
1 x2 解：首先将二阶方程写为一阶方程组的形式：
的数值解，其解析解为： y
y' z , y 1 1 3 z' 2 y , z 1 1 再编写函数的 m 文件： %这是一个二阶微分方程 function f=dyfun1(t,y) f(1)=y(2);%y(1)是 y,y(2)是 z。 f(2)=2*(y(1)*y(1)*y(1)); f=f(:);
直接使用 ode45 来计算： close all; clear all; clc; [x,y]=ode45('dyfun1',[1 1.5]',[-1 -1]') plot(x,(y(:,1))); 其结果如图：
-0.8
-1
-1.2
-1.4
-1.6
-1.8
-2
-2.2
1
1.05
1.1
1.15
1.2
2
解析解
1.5
1
0.5
y(x) z(x)
0
-0.5 -50
0
50
100
150
200
250
300
350
400
450
500
2
数值解
1.5ቤተ መጻሕፍቲ ባይዱ
1
0.5
y(x) z(x)
0
-0.5 -50
0
50
100
150
200
250
300
350
400
450
500
例 2：考虑下面的 Lorenz 方程组
dx dt x y dy x y xz dt dz xy z dt 参数 , , 适当的取值会使得系统趋于混沌状态。取 =30, =2.8, =12，利 用经典四阶龙格-库塔法求其数值解，并绘制 z 随 x 变化的曲线。 解：首先编写函数的 m 文件： %mafun.m function ff=mafun(t,y) b=2.8;r=30;sigma=12; ff(1)=-sigma*y(1)+sigma*y(2); ff(2)=r*y(1)-y(2)-y(1)*y(3); ff(3)=y(1)*y(2)-b*y(3); ff=ff(:);
微分方程组的龙格-库塔公式求解 matlab 版
南京大学 王寻
1. 一阶常微分方程组 考虑方程组 y' f x , y , z , y x0 y0 z' g x , y , z , z x0 z0 其经典四阶龙格-库塔格式如下： 对于 n = 0, 1, 2,...,计算 h yn 1 yn 6 K1 2 K 2 2 K 3 K 4 h z n 1 z n L1 2 L2 2 L3 L4 6 其中 K1 K 2 K 3 K 4
f xn , yn , z n , L1 g xn , yn , z n h hK hL h hK hL f xn , yn 1 , z n 1 , L2 g xn , yn 1 , z n 1 2 2 2 2 2 2 h hK 2 hL2 h hK 2 hL2 f xn , y n , zn , zn , L3 g xn , yn 2 2 2 2 2 2 f xn h , yn hK 3 , z n hL3 , L4 g xn h , yn hK 3 , z n hL3
1.25
1.3
1.35
1.4
1.45
1.5
这与解析解作出来的图像很接近。