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

微分方程组的龙格-库塔公式求解 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(dyfu