常微分方程数值解实验报告材料.doc
实验一：常微分方程的数值解法1、 分别用Euler法、改良的Euler法预报校正格式和SK法求解初值问题。h=0.1并与真解作比拟。1.1实验代码：%欧拉法function x,y=naeuler(dyfun,xspan,y0,h)%dyfun是常微分方程，xspan是x的取值X围，y0是初值，h是步长x=xspan(1):h:xspan(2);y(1)=y0; for n=1:length(x)-1 y(n+1)=y(n)+h*feval(dyfun,x(n),y(n);end%改良的欧拉法function x,m,y=naeuler2(dyfun,xspan,y0,h)%dyfun是常微分方程，xspan是x的取值X围，y0是初值，h是步长。%返回值x为x取值，m为预报解，y为校正解x=xspan(1):h:xspan(2);y(1)=y0; m=zeros(length(x)-1,1);for n=1:length(x)-1 k1=feval(dyfun,x(n),y(n);y(n+1)=y(n)+h*k1; m(n)=y(n+1);k2=feval(dyfun,x(n+1),y(n+1);y(n+1)=y(n)+h*(k1+k2)/2;end%四阶SK法function x,y=rk(dyfun,xspan,y0,h)%dyfun是常微分方程，xspan是x的取值X围，y0是初值，h是步长。x=xspan(1):h:xspan(2);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*k1)/2); k3=feval(dyfun,x(n)+h/2,y(n)+(h*k2)/2); k4=feval(dyfun,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;dyfun=inline('-y+x+1'); x1,y1=naeuler(dyfun,0,1,1,0.1);x2,m,y2=naeuler2(dyfun,0,1,1,0.1);x3,y3=rk(dyfun,0,1,1,0.1);plot(x,y,'r',x1,y1,'+',x2,y2,'*',x3,y3,'o');xlabel('x');ylabel('y');legend('y为真解','y1为欧拉解','y2为改良欧拉解','y3为SK解','Location','NorthWest');1.2实验结果：x真解y欧拉解y1预报值m校正值y2SK解y30.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 2、 选取一种理论上收敛但是不稳定的算法对问题1进展计算，并与真解作比拟。选改良的欧拉法2.1实验思路：算法的稳定性是与步长h密切相关的。而对于问题一而言，取定步长h=0.1不论是单步法或低阶多步法都是稳定的算法。所以考虑改变h取值X围，借此分析不同步长会对结果造成什么影响。故依次采用h=2.0、2.2、2.4、2.6的改良欧拉法。2.2实验代码：%主程序x=0:3:30;y=exp(-x)+x;dyfun=inline('-y+x+1'); x1,m1,y1=naeuler2(dyfun,0,20,1,2);x2,m2,y2=naeuler2(dyfun,0,22,1,2.2);x3,m3,y3=naeuler2(dyfun,0,24,1,2.4);x4,m4,y4=naeuler2(dyfun,0,26,1,2.6);subplot(2,2,1)plot(x,y,'r',x1,y1,'+');xlabel('h=2.0');subplot(2,2,2)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('h=2.6');2.3实验结果：x0.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 6.9904 8.3684 0.3 7.0000 8.4158 10.4418 13.4398 0.4 9.0000 11.0153 14.3979 20.4388 0.5 11.0000 13.7027 19.1008 30.8690 0.6 13.0000 16.4973 24.9092 47.4068 0.7 15.0000 19.4227 32.3536 74.8161 0.8 17.0000 22.5077 42.2194 121.5767 0.9 19.0000 25.7874 55.6687 202.7825 1.0 21.0000 29.3046 74.4217 345.3008 实验结果分析：从实验1结果可以看出，在算法满足收敛性和稳定性的前提下，Eluer法虽然计算并不复杂，但凡精度不足，反观改良的Eluer法和SK法虽然计算略微复杂但是结果很准确。实验2改变了步长，导致算法理论上收敛但是不满足稳定性。结果表示步长h越大，结果越失真。对于同一个问题，步长h的选取变得尤为重要，这三种单步法的绝对稳定区间并不一样，所以并没有一种方法是万能的，我们应该根据不同的步长来选取适宜的方法。实验二：Ritz-Galerkin方法与有限差分法1、 用中心差分格式求解边值问题取步长h=0.1，并与真解作比拟。1.1实验代码：%中心差分法function U=fdm(xspan,y0,y1,h)%xspan为x取值X围，y0,y1为边界条件，h为步长N=1/h;d=zeros(1,N-1);for i=1:N x(i)=xspan(1)+i*h; q(i)=1; f(i)=x(i);endfor 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-2 b(i+1,i)=-1;endfor i=1:N-2 c(i,i+1)=-1;endA=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= inv(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('y');legend('y真解','y1中心差分法','Location','NorthWest');1.2实验结果：xy真解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-Galerkin方法求解上述问题，并与真值作比拟，列表画图。2.1实验代码：%Ritz_Galerkin法functionvu=Ritz_Galerkin(x0,y0,x1,y1,h)%x0,x1为x取值X围，y0,y1为边界条件，h为步长N=1/h;syms x;for i=1:N fai(i)=x*(1-x)*(x(i-1); dfai(i)=diff(x*(1-x)*(x(i-1); endfor i=1:N for j=1:N fun=dfai(i)*dfai(j)+fai(i)*fai(j); A(i,j)=int(fun,x,0,1); end fun=x*fai(i)+dfai(i); f(i)=int(fun,x,0,1);endc=inv(A)*f'product=c.*fai' sum=0; for i=1:N sum=sum+product(i);endvu=;for y=0:h:1 v=subs(sum,x,y); vu=vu,v; endy=0:h:1;yy=0:0.1:1; u=sin(yy)/sin(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(2)-1);y1=Ritz_Galerkin(0,0,1,0,0.1);y1=double(y1);plot(x,y,'r',x,y1,'+')xlabel('x');ylabel('y');legend('y为真解','y1为RG法','Location','NorthWest');2.2实验结果：xy真解y1RG法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;如此上述问题的数值解法怎样变化？结果如何？程序运算出来真解与数值解完全一样。其值为y=x。具体运算不再赘述。实验结果分析：对于实验1、2，我们可以看出不论是有限差分法还是Ritz-Galerkin法都非常稳定，结果非常准确误差小于。对于实验3，编程中确定系数矩阵和常数项是最重要的。确定过程中，要注意matlab中循环是从1开始的，而我们推导的公式中循环是从0开始的。所以要区分开来慎重对待，不然会产生极大地误差。
