微分方程作业.doc

P10习题1.用Euler法和改进的Euler法求u=-5u (0t1),u(0)=1的数值解，步长h=0.1,0.05;并比较两个算法的精度。解：function du=Euler_fun1(t,u)du=-5*u;clear;h=0.1;tend=1;N=1/h;t(1)=0;u(1)=1;t=h.*(0:N);for n=1:N u(n+1)=u(n)+h*Euler_fun1(t(n),u(n);endplot(t,u,*);hold onfor n=1:N v(1)=u(n)+h*Euler_fun1(t(n),u(n); for k=1:6v(k+1)=u(n)+h/2*(Euler_fun1(t(n),u(n)+Euler_fun1(t(n+1),v(k); end u(n+1)=v(k+1);endplot(t,u,o);sol=dsolve(Du=-5*u,u(0)=1);u_real=eval(sol);plot(t,u_real,r);将上述 h 换为0.05得：由图像知道:显然改进的Euler法要比Euler法精确度要高；3.将u=-u（0t1)，u(0)=0,u(0)=1化为一阶方程组，并用Euler法和改进的的Euler法求解，步长h=0.1,0.05；并比较两个算法的精度。解：function du=fun31(y)du=y;function dy=fun32(u)dy=-u;clear;h=0.1;tend=1;N=1/h;t(1)=0;u(1)=0;y(1)=;t=h.*(0:N);for n=1:N u(n+1)=u(n)+h*y(n); y(n+1)=y(n)+h*(-u(n);endplot(t,u,*);hold onfor n=1:N v(1)=u(n)+h*fun31(y(n); w(1)=y(n)+h*fun32(u(n);for k=1:6v(k+1)=u(n)+h/2*(fun31(y(n)+fun31(.w(k);w(k+1)=y(n)+h/2*(fun32(u(n)+fun32(.v(k); end u(n+1)=v(k+1); y(n+1)=w(k+1);endplot(t,u,o);sol=dsolve(D2u=-u,u(0)=0,Du(0)=1;u_real=eval(sol);plot(t,u_real,r);将上述 h 换为0.05得：由图像可以知道:显然改进的Euler法要比Euler法精确度要高；实习题（二）1.取步长 ，分别用Euler 法和改进的Euler 法求下列初值问题的解，并与真解相比较.（1）真解 ；解：function du=fun1(x,u)du=u-2*x/u;clear;h=0.1;xend=1;N=1/h;x(1)=0;u(1)=1;x=h.*(0:N);%Eluer法%for n=1:N u(n+1)=u(n)+h*fun1(x(n),u(n);endplot(x,u,*);hold on%改进的Eluer法%for n=1:N v(1)=u(n)+h*fun1(x(n),u(n); for k=1:6 v(k+1)=u(n)+h/2*(fun1(x(n),u(n)+fun1(x(n+1),v(k); end u(n+1)=v(k+1);endplot(x,u,o);hold on%真解%u_real=sqrt(1+2*x);plot(x,u_real,r);由图像可以知道:显然改进的Euler法要比Euler法精确度要高；（2）真解 ；解：function du=fun2(x,u)du=(u/x)-x.2/u.2;clear;h=0.1;N=1/h;x=1:h:2;x(1)=1;u(1)=2;for n=1:N u(n+1)=u(n)+h*fun2(x(n),u(n);endplot(x,u,*);hold onfor n=1:N v(1)=u(n)+h*fun2(x(n),u(n); for k=1:6 v(k+1)=u(n)+h/2*(fun2(x(n),u(n)+fun2(x(n+1),v(k); end u(n+1)=v(k+1);endplot(x,u,o);hold onu_real=x.*(8-3.*log(x).(1/3);plot(x,u_real,r);由图像可知：改进的Euler法和Euler法都很接近真值。（3）真解 .解：function du=fun3(x,u)du=u/(2*x)-x/(2*u2);clear;h=0.1;N=0.5/h;x=1:h:1.5;x(1)=1;u(1)=1;for n=1:N u(n+1)=u(n)+h*fun3(x(n),u(n);endplot(x,u,*);hold onfor n=1:N v(1)=u(n)+h*fun3(x(n),u(n); for k=1:6 v(k+1)=u(n)+h/2*(fun3(x(n),u(n)+fun3(x(n+1),v(k); end u(n+1)=v(k+1);endplot(x,u,o);hold onu_real=(4*x.(3/2)-3*x.2).(1/3);plot(x,u_real,r);由图像可以知道:显然改进的Euler法要比Euler法精确度要高；2.试用预报校正格式(1.20)解初值问题并与Euler格式比较精度，取h=0.1。作业要求：写出程序，列表或用图形显示结果，并给出图或表所说明的结果。解：function du=Euler_fun2(t,u)du=-u+t+1;clear;h=0.1;tend=1;N=1/h;t(1)=0;u(1)=1;t=h.*(0:N);for n=1:N u(n+1)=u(n)+h*Euler_fun2(t(n),u(n);endplot(t,u,*);hold onfor n=1:N u0(n+1)=u(n)+h*Euler_fun2(t(n),u(n); u(n+1)=u(n)+h/2*(Euler_fun2(t(n),u(n)+Euler_fun2(t(n+1),u0(n+1);endplot(t,u,o);hold onsol=dsolve(Du=-u+t+1,u(0)=1);u_real=eval(sol);plot(t,u_real,r);由图像可以知道:显然预报校正格式要比Euler法精确度要高；P37 例4.1 用四级四阶Runge-Kutta法计算初值问题： u=4tu0.5,0t2, u(0)=1.取h=0.1,0.5,1.精确解为 u(t)=(1+t2)2作业要求：写出程序，列表或用图形显示结果，并给出图或表所说明的结果.解：function du=fun4(t,u)du=4*t*u.(1/2);clear;h=0.1;N=2/h;t=0:h:2;t(1)=0;u(1)=1;for n=1:N k1=fun4(t(n),u(n); k2=fun4(t(n)+0.5*h,u(n)+0.5*h*k1); k3=fun4(t(n)+0.5*h,u(n)+0.5*h*k2); k4=fun4(t(n)+h,u(n)+h*k3); u(n+1)=u(n)+h*(k