微分方程作業(yè)

1、-作者xxxx-日期xxxx微分方程作業(yè)【精品文檔】P10習題1.用Euler法和改進的Euler法求u=-5u (0t1),u(0)=1的數(shù)值解，步長h=0.1,0.05;并比較兩個算法的精度?！揪肺臋n】解：function du=Euler_fun1(t,u)du=-5*u;clear;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

2、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法精確度要高；u=-u（0t1)，u(0)=0,u(0)=1化為一階方程組，并用Euler法和改進的的Euler法求解，步長h=0.1,0.05；并比較兩個算法的精度。解：function du=fun

3、31(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)+f

4、un32(.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

5、(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法要比Eule

6、r法精確度要高；（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

7、=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

8、(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。作業(yè)要求：寫出程序，列表或用圖形顯示結果，并給出圖或表所說明的結果。解：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

9、.*(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法精確度要高

10、；P37 例4.1 用四級四階Runge-Kutta法計算初值問題： u=4tu,0t2, u(0)=1.取h=0.1,0.5,1.精確解為 u(t)=(1+t2)2作業(yè)要求：寫出程序，列表或用圖形顯示結果，并給出圖或表所說明的結果.解：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*(