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

1、常微分方程数值解实验报告学院：数学与信息科学专业：信息与计算科学.jz*:思义学号：201216524课程：常微分方程数值解实验一：常微分方程的数值解法1、分别用Euler法、改良的Euler法预报校正格式和SK法求解初值问题。h=0.1并与真解作比拟。y'yx1y(0)11.1实验代码：%欧拉法functionx,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

2、(n),y(n);end%改良的欧拉法functionx,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;.jz*m(n)=y(n+1);k2=feval(dyfun,x(n+1),y(n+1);y(n+1)=y(n)+h*(k1+k2)/2;en

3、d%四阶SK法functionx,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%主

4、程序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为改

5、良欧拉解,'y3为SK解','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.1065.jz*0.61.14881

6、.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

7、.0、2.2、2.4、2.6的改良欧拉法。2.2 实验代码：.jz*%主程序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=

8、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.000

9、00.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.3

10、008.jz*实验结果分析：从实验1结果可以看出，在算法满足收敛性和稳定性的前提下，日uer法虽然计算并不复杂，但凡精度缺乏，反观改良的Eluer法和S-K法虽然计算略微复杂但是结果很准确。实验2改变了步长，导致算法理论上收敛但是不满足稳定性。结果表示步长h越大，结果越失真。对于同一个问题，步长h的选取变得尤为重要，这三种单步法的绝对稳定区间并不一样，所以并没有一种方法是万能的，我们应该根据不同的步长来选取适宜的方法。实验二：Ritz-Galerkin方法与有限差分法1、用中心差分格式求解边值问题.jz*y''yx0,0x1y(0)0,y(i)0取步长h=0.1,并与真解作比拟

11、。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)

12、*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');.jz*1

13、.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.0000.jz*2、用Ritz-Galerkin方法求解上述问题，并与真值作比拟，列表画图2.1 实验代码：%Ritz_Galerkin法functionvu=Ritz_Galerkin(x0,y0,x1,y1,h)%x0,x1为x取值围，y

14、0,y1为边界条件，h为步长N=1/h;symsx;fori=1:Nfai(i)=x*(1-x)*(xA(i-1);dfai(i)=diff(x*(1-x)*(xA(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

15、(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);%主程序.jz*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真解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.000