微分方程组的龙格-库塔公式求解 matlab 版 南京大学 王寻 1. 一阶常微分方程组 考虑方程组          00 00 zxz,z , y, xg' z yxy,z , y, xf' y 其经典四阶龙格-库塔格式如下： 对于 n = 0, 1, 2,.,计算           43211 43211 22 6 22 6 LLLL h zz KKKK h yy nn nn 其中                                        334334 22 3 22 3 11 2 11 2 11 222222 222222 hLz ,hKy, hxgL,hLz ,hKy, hxfK hL z , hK y, h xgL, hL z , hK y, h xfK hL z , hK y, h xgL, hL z , hK y, h xfK z ,y,xgL,z ,y,xfK nnnnnn nnnnnn nnnnnn nnnnnn 下面给出经典四阶龙格-库塔格式解常微分方程组的 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 如下为例题： 例 1：取 h=0.02，利用程序 marunge4s.m 求刚性微分方程组          10100 209999010 z, z' z ,y, z.y.' y 的数值解，其解析解为： xxx. ez,eey 100100010  解：首先编写 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('解析解'); 可以看出数值解和解析解的运算结果误差很小： -50050100150200250300350400450500 -0.5 0 0.5 1 1.5 2 y(x) z(x) 解 析 解 -50050100150200250300350400450500 -0.5 0 0.5 1 1.5 2 y(x) z(x) 数 值 解 例 2：考虑下面的 Lorenz 方程组             zxy dt dz xzyx dt dy yx dt dx    参数,,适当的取值会使得系统趋于混沌状态。取=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(:); 再编写运行函数： clear all; close all; clc; [t,y]=marunge4s(@mafun,[0 500],[0 1 2],0.005); plot(y(1,:),y(3,:),'r'); 得到如下图所示的结果： -25-20-15-10-50510152025 0 10 20 30 40 50 60 2. 高阶微分方程组 对于高阶微分方程，总是可以化成方程组的形式。例如，二阶方程          0000 ' yx' y,yxy ' y, y, xg“ y 总是可以化为一阶方程组：             00000 z' yxz ,yxy z , y, xg' z z' y 因此不需要对于高阶方程给出计算公式。 例 3：求二阶方程         111 5112 3 ' yy .x,y“ y 的数值解，其解析解为： 2 1   x y 解：首先将二阶方程写为一阶方程组的形式：          112 11 3 z,y' z y, z' y 再编写函数的 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))); 其结果如图： 11.051.11.151.21.251.31.351.41.451.5 -2.2 -2 -1.8 -1.6 -1.4 -1.2 -1 -0.8 这与解析解作出来的图像很接近。

展开阅读全文