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

微分方程组的龙格-库塔公式求解matlab版 南京大学 王寻 1. 一阶常微分方程组 考虑方程组     yf x,y,z , y x y  0 0     zg x,y,z , z x z  0 0 其经典四阶龙格-库塔格式如下： 对于n=0,1,2,...,计算  h y y  K 2K 2K K   n1 n 6 1 2 3 4  h z z  L 2L 2L L  n1 n 1 2 3 4  6 其中 K f x ,y ,z , L g x ,y ,z    1 n n n 1 n n n   h hK1 hL1  h hK1 hL1 K f x  ,y  ,z 2  n n n ,L g x  ,y  ,z 2  n n n    2 2 2   2 2 2   K f x  ,y 3  n h n hK2 ,z n hL2 ,L g x  ,y 3  n h n hK2 ,z n hL2    2 2 2   2 2 2  K f x h,y hK ,z hL ,L g x h,y hK ,z hL     4 n n 3 n 3 4 n n 3 n 3 下面给出经典四阶龙格-库塔格式解常微分方程组 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求刚性微分方程组   y0.01y99.99z, y 0 2,    z100z, z 0 1 的数值解，其解析解为： 0.01x 100x 100x y e e ,ze 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(:); 然后编