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

常微分方程数值解实验报告学院：数学与信息科学专业：信息与计算科学姓名：郑思义学号：201216524课程：常微分方程数值解实验一：常微分方程的数值解法1、分别用Euler法、改进的Euler法（预报校正格式）和S—K法求解初值问题。（h=0.1）并与真解作比较。10(1y'）yxy1.1实验代码：%欧拉法function[x,y]=naeuler(dyfun,xspan,y0,h)%dyfun是常微分方程，xspan是x的取值范围，y0是初值，h是步长x=xspan(1):h:xspan(2);y(1)=y0;forn=1:length(x)-1y(n+1)=y(n)+h*feval(dyfun,x(n),y(n));end%改进的欧拉法function[x,m,y]=naeuler2(dyfun,xspan,y0,h)%dyfun是常微分方程，xspan是x的取值范围，y0是初值，h是步长。%返回值x为x取值，m为预报解，y为校正解x=xspan(1):h:xspan(2);y(1)=y0;m=zeros(length(x)-1,1);forn=1:length(x)-1k1=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%四阶S—K法function[x,y]=rk(dyfun,xspan,y0,h)%dyfun是常微分方程，xspan是x的取值范围，y0是初值，h是步长。x=xspan(1):h:xspan(2);y(1)=y0;forn=1:length(x)-1k1=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为S—K解','Location','NorthWest');1.2实验结果：x真解y欧拉解y1预报值m校正值y2S—K解y30.01.00001.00001.00001.00000.11.00481.00001.00001.00501.00480.21.01871.01001.01451.01901.01870.31.04081.02901.03711.04121.04080.41.07031.05611.06711.07081.07030.51.10651.09051.10371.10711.10650.61.14881.13141.14641.14941.14880.71.19661.17831.19451.19721.19660.81.24931.23051.24751.25001.24930.91.30661.28741.30501.30721.30661.01.36791.34871.36651.36851.36792、选取一种理论上收敛但是不稳定的算法对问题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;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实验结果：xh=2.0h=2.2h=2.4h=2.60.01.00001.00001.00001.00000.13.00003.42003.88004.38000.25.00005.88846.99048.36840.37.00008.415810.441813.43980.49.000011.015314.397920.43880.511.000013.702719.100830.86900.613.000016.497324.909247.40680.715.000019.422732.353674.81610.817.000022.507742.2194121.57670.919.000025.787455.6687202.78251.021.000029.304674.4217345.3008实验结果分析：从实验1结果可以看出，在算法满足收敛性和稳定性的前提下，Eluer法虽然计算并不复杂，凡是精度不足，反观改进的Eluer法和S—K法虽然计算略微复杂但是结果很精确。实验2改变了步长，导致算法理论上收敛但是不满足稳定性。结果表示步长h越大，结果越失真。对于同一个问题，步长h的选取变得尤为重要，这三种单步法的绝对稳定区间并不一样，所以并没有一种方法是万能的，我们应该根据不同的步长来选取合适的方法。实验二：Ritz-Galerkin方法与有限差分法1、用中心差分格式求解边值问题0)1(,0)0(10,0''yyxxyy取步长h=0.1，并与真解作比较。1.1实验代码：%中心差分法functionU=fdm(xspan,y0,y1,h)%xspan为x取值范围，y0,y1为边界条件，h为步长N=1/h;d=zeros(1,N-1);fori=1:Nx(i)=xspan(1)+i*h;q(i)=1;f(i)=x(i);endfori=1:N-1d(i)=q(i)*h*h+2;enda=diag(d);b=zeros(N-1);c=zeros(N-1);fori=1:N-2b(i+1,i)=-1;endfori=1:N-2c(i,i+1)=-1;endA=a+b+c;fori=2:N-2B(i,1)=f(i)*h*h;endB(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.00.00000.00000.10.01480.01480.20.02870.02870.30.04090.04080.40.05050.05040.50.05660.05650.60.05830.05820.70.05450.05450.80.04430.04430.90.02650.02651.00.00000.00002、用Ritz-Galerkin方法求解上述问题，并与真值作比较，列表画图。2.1实验代码：%Ritz_Galerkin法functionvu=Ritz_Galerkin(x0,y0,x1,y1,h)%x0,x1为x取值范围，y0,y1为边界条件，h为步长N=1/h;symsx;fori=1:Nfai(i)=x*(1-x)*(x^(i-1));dfai(i)=diff(x*(1-x)*(x^(i-1)));endfori=1:Nforj=1:Nfun=dfai(i)*dfai(j)+fai(i)*fai(j);A(i,j)=int(fun,x,0,1);endfun=x*fai(i)+dfai(i);f(i)=int(fun,x,0,1);endc=inv(A)*f';product=c.*fai';sum=0;fori=1:Nsum=sum+product(i);endvu=[];fory=0:h:1v=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为R—G法','Location','NorthWest');2.2实验结果：xy真解y1R—G法0.00.00000.00000.10.01480.01480.20.02870.02870.30.04090.04090.40.05050.05050.50.05660.05660.60.05830.05830.70.05450.05450.80.04430.04430.90.02650.02651.00.00000.00003、若边值条件为y(0)=0,y(1)=1;则上述问题的数值解法怎样变化？结果如何？程序运算出来真解与数值解完全一样。其值为y=x。（具体运算不再赘述）。实验结果分析：对于实验1、2，我们可以看出不论是有限差分法还是Ritz-Galerkin法都非常稳定，结果非常精确（误差小于0.0001）。对于实验3，编程中确定系数矩阵和常数项是最重要的。确定过程中，要注意matlab中循环是从1开始的，而我们推导的公式中循环是从0开始的。所以要区分开来谨慎对待，不然会产生极大地误差。