常微分方程数值解实验报告

1、常微分方程数值解实验报告 学院：数学与信息科学 专业：信息与计算科学 姓名：郑思义 学号:201216524 课程:常微分方程数值解 实验一：常微分方程的数值解法 1、分别用Euler法、改进的Euler法(预报校正格式)与S K法求解 初值问题。(h=0、1)并与真解作比较。 y y x 1 y(0) 1 1、1实验代码： %欧拉法 fun cti onx,y =n aeuler(dyfu n, xspa n,y 0,h) %dyfun就是常微分方程，xspan就是x的取值范围，y0就是初值,h就是步长 x=xspa n( 1):h:xspa n( 2); y(1)=y0; for n=1:

2、le ngth(x)-1 y( n+1)=y( n)+h*feval(dyfu n, x( n),y( n); end %改进的欧拉法 fun cti onx,m, y=n aeuler2(dyfu n, xspa n,y 0,h) %dyfun就是常微分方程，xspan就是x的取值范围，y0就是初值,h就是步长。 %返回值x为x取值,m为预报解，y为校正解 x=xspa n( 1):h:xspa n( 2); y(1)=y0; m=zeros(le ngth(x)-1,1); for n=1:le ngth(x)-1 k1=feval(dyfu n,x(n ),y (n); y( n+1)=

3、y( n)+h*k1; m(n)=y(n+1); k2=feval(dyfu n,x (n+1),y( n+1); y( n+1)=y( n)+h*(k1+k2)/2; end %四阶SK法 fun cti onx,y=rk(dyfu n, xspa n,y 0,h) %dyfun就是常微分方程，xspan 就是x的取值范围，y0就是初值,h就是步长。 x=xspa n( 1):h:xspa n( 2); y(1)=y0; for n=1:le ngth(x)-1 k1=feval(dyfu n,x(n ),y (n); k2=feval(dyfu n,x( n)+h/2,y( n)+(h*k

4、1)/2); k3=feval(dyfu n,x( n)+h/2,y( n)+(h*k2)/2); k4=feval(dyfu n,x( n)+h,y( n)+h*k3); y(n+1)=y( n)+(h/6)*(k1+2*k2+2*k3+k4); end %主程序 x=0:0、1:1;y=exp(-x)+x; dyfu n=in li ne(-y+x+1); x1,y1=naeuler(dyfu n,0,1,1,0、1); x2,m,y2=naeuler2(dyfu n,0,1,1,0、1); x3,y3=rk(dyfu n,0,1,1,0、1); plot(x,y, r,x1,y1,+,x

5、2,y2,* ,x3,y3,o: ); xlabel( x );ylabel(y); lege nd( y 为真解,y1为欧拉解,y2 为改进欧拉解 ,y3 为 S K 解 ,Location ,NorthWest); 1、2实验结果: x 真解y 欧拉解y1 预报值m 校正值y2 S K 解 y3 0、0 1、 0000 1、 0000 1、 0000 1、 0000 0、1 1、 0048 1、 0000 1、 0000 1、 0050 1、 0048 0、2 1、 0187 1、 0100 1、 0145 1、 0190 1、 0187 0、3 1、 0408 1、 0290 1、 03

6、71 1、 0412 1、 0408 0、4 1、 0703 1、 0561 1、 0671 1、 0708 1、 0703 0、5 1、 1065 1、 0905 1、 1037 1、 1071 1、 1065 0、6 1、 1488 1、 1314 1、 1464 1、 1494 1、 1488 0、7 1、 1966 1、 1783 1、 1945 1、 1972 1、 1966 0、8 1、 2493 1、 2305 1、 2475 1、 2500 1、 2493 0、9 1、 3066 1、 2874 1、 3050 1、 3072 1、 3066 1、01、36791、34871、

7、36651、36851、3679 2、选取一种理论上收敛但就是不稳定的算法对问题 1进行计算，并 与真解作比较。（选改进的欧拉法） 2、1实验思路：算法的稳定性就是与步长h密切相关的。而对于问 题一而言,取定步长h=0、1不论就是单步法或低阶多步法都就是稳定 的算法。所以考虑改变h取值范围，借此分析不同步长会对结果造成 什么影响。故依次采用h=2、0、2、2、2、4、2、6的改进欧拉法 2、2实验代码： %主程序 x=0:3:30;y=exp(-x)+x; dyfu n=in li ne(-y+x+1); x1,m1,y1=naeuler2(dyfu n,0,20,1,2); 、2); 、4)

8、; 、6); x2,m2,y2=naeuler2(dyfu n,0,22,1,2 x3,m3,y3=naeuler2(dyfu n,0,24,1,2 x4,m4,y4=naeuler2(dyfu n,0,26,1,2 subplot(2,2,1) plot(x,y, r,x1,y1,+ );xlabel( h=2 、0); plot(x,y,r ,x2,y2, + );xlabel( h=2 、2 subplot(2,2,3) plot(x,y,r ,x3,y3, + );xlabel( h=2 、4 subplot(2,2,4) plot(x,y,r ,x4,y4. + );xlabel(

9、h=2 、6 subplot(2,2,2) 2、3实验结果： x h=2、0 h=2、2 h=2、4 h=2、6 0、 0 1、 0000 1、 0000 1、 0000 1、 0000 0、 1 3、 0000 3、4200 3、 8800 4、 3800 0、 2 5、 0000 5、8884 d=zeros(1,N-1); for i=1:N x(i)=xspan(1)+i*h; q(i)=1; f(i)=x(i); end for i=1:N-1 d(i)=q(i)*h*h+2; end a=diag(d); b=zeros(N-1); c=zeros(N-1); for i=1:N-

10、2 b(i+1,i)=-1; end for i=1:N-2 c(i,i+1)=-1; end A=a+b+c; for i=2:N-2 B(i,1)=f(i)*h*h; end B(1,1)=f(1)*h*h+y0; B(N-1,1)=f(N-1)*h*h+y1; U= in v(A)*B; %主程序 x=0:0、1:1; y=x+(exp(1)*exp(-x)/(exp(2)-1)-(exp(1)*exp(x)/(exp(2)-1); y1=fdm(0,1,0,0,0 、1); y1=0,y1,0; plot(x,y,r,x,y1, + ) xlabel(x );ylabel( ) /);

11、 legend(y 真解，y1 中心差分法 ,Location, NorthWest ) 1、2实验结果 : x y 真解 y1中心差分法 0、0 0、 0000 0、 0000 0、1 0、 0148 0、 0148 0、2 0、 0287 0、 0287 0、3 0、 0409 0、 0408 0、4 0、 0505 0、 0504 0、5 0、 0566 0、 0565 0、6 0、 0583 0、 0582 0、7 0、 0545 0、 0545 0、8 0、 0443 0、 0443 0、9 0、 0265 0、 0265 1、0 0、 0000 0、 0000 2、用Ritz-Ga

12、lerkin方法求解上述问题，并与真值作比较，列表画图 2、1实验代码： %Ritz_Galerki n法 fun cti onvu=Ritz_Galerki n(x0,y0,x1,y1,h) %x0,x1为x取值范围，yO,y1 为边界条件,h为步长 N=1/h; syms x; for i=1:N fai(i)=x*(1-x)*(xA(i-1); dfai(i)=diff(x*(1-x)*(xA(i-1); end for i=1:N for j=1:N fun=dfai(i)*dfai(j)+fai(i)*fai(j); A(i,j)=i nt(fun, x,0,1); end fun=

13、x*fai(i)+dfai(i); f(i)=i nt(fu n, x,0,1); end c=i nv(A)*f; product=c 、*fai; sum=0; for i=1:N sum=sum+product(i); end VU=; for y=0:h:1 v=subs(sum,x,y); vu=vu,v; end y=0:h:1; yy=0:0、1:1; u=si n(yy)/si n(1)-yy; u=vpa(u,5); vu=vpa(vu,5); %主程序 x=0:0、1:1; y=x+(exp(1)*exp(-x)/(exp(2)-1)-(exp(1)*exp(x)/(exp

14、(2)-1); y1=Ritz_Galerki n(0,0,1,0,0、1); y1=double(y1); plot(x,y, r,x,y1,+) xlabel( x );ylabel( y); legend( y 为真解,y1 2、2实验结果： 为 R G 法，Location ,NorthWest); x y真解 y1R G 法 0、0 0、 0000 0、 0000 0、1 0、 0148 0、 0148 0、2 0、 0287 0、 0287 0、3 0、 0409 0、 0409 0、4 0、 0505 0、 0505 0、5 0、 0566 0、 0566 0、6 0、 0583 0、 0583 0、7 0、 0545 0、 0545 0、8 0、 0443 0、 0443 0、9 0、 0265 0、 0265 1、0 0、 0000 0、 0000 3、若边值条件为y(0)=0,y(1)=1;则