第八章 常微分方程及初值问题1.ppt

function [x,y]=RK4(fun,x0,y0,xmax,h) %fun为显式一阶微分方程右端的函数 %x0,y0为初始条件,满足y(x0)=y0 %xmax为x的取值最大值 %h为将x等分后的步长 N=(xmax-x0)/h+1;%N为总的节点数 x(1)=x0;y(1)=y0; for n=1:N-1 x(n+1)=x(n)+h; k1=h*feval(fun,x(n),y(n)); k2=h*feval(fun,x(n)+1/2*h,y(n)+k1/2); k3=h*feval(fun,x(n)+1/2*h,y(n)+k2/2); k4=h*feval(fun,x(n+1),y(n)+k3); y(n+1)=y(n)+1/6*(k1+2*k2+2*k3+k4); end fun=inline(y+2*x/y^2,x,y); [x1,y1]=Euler_mend(fun,0,1,2,0.1); [x2,y2]=RK4(fun,0,1,2,0.1); x=0:0.1:2; y=1/3*(-18-54*x+45*exp(3*x)).^(1/3); plot(x1,y1,*,x2,y2,s,x,y) axis([0,2,0,9]) 再分别调用Euler_mend.m和RK4.m求解： 从图形上看，似乎二阶龙格—库塔法与四阶龙格-库塔法同样接近解析解，故再从数值结果比较看看哪种误差小 err=[abs(y-y1),abs(y-y2)] 结果为 err = 0 0 0.0009 0.0000 0.0016 0.0000 0.0022 0.0000 0.0026 0.0000 0.0031 0.0000 0.0035 0.0000 0.0040 0.0000 0.0047 0.0000 0.0054 0.0000 0.0062 0.0000 0.0073 0.0000 0.0084 0.0000 0.0098 0.0000 0.0115 0.0000 0.0133 0.0000 0.0155 0.0000 0.0181 0.0000 0.0210 0.0000 0.0243 0.0000 0.0281 0.0000 四、二阶ODE问题 二阶ODE的一般形式为： 其中 是常数或是 的函数，后两个方程为初 始条件。 如果 与u无关，则该方程称为线性ODE; 如果三者之中有一个是u或 的函数，称为非线性的。 下面着重研究用向前Euler法求解二阶ODE，及MATLAB程序。 所以二阶ODE分解为两个一阶ODE： 设： 定义 则上述两个一阶ODE写为： 其向前Euler法的格式为： 例 求解二阶ODE 解：设 令 则原方程的向量形式为： 向前Euler法的格式为： h = 0.05; t_max=5; n=1; y(:,1) = [0; 1];t(1) = 0; while t(n) t_max y(:,n+1) = y(:,n) + h*f_def(y(:,n),t); t(n+1) = t(n)+h ; n=n+1; end axis([0 5 -1 1]);plot(t, y(1,:),