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

常微分方程数值解实验报告 学院：数学与信息科学 专业：信息与计算科学 姓名：郑思义 学号：201216524 课程：常微分方程数值解实验一：常微分方程的数值解法1、 分别用Euler法、改进的Euler法（预报校正格式）和SK法求解初值问题。（h=0.1）并与真解作比较。1.1实验代码：%欧拉法function x,y=naeuler(dyfun,xspan,y0,h)%dyfun是常微分方程，xspan是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的取值范围，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的取值范围，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 1.00001.00001.00001.00000.1 1.00481.00001.00001.00501.00480.2 1.01871.01001.01451.01901.01870.3 1.04081.02901.03711.04121.04080.4 1.07031.05611.06711.07081.07030.5 1.10651.09051.10371.10711.10650.6 1.14881.13141.14641.14941.14880.7 1.19661.17831.19451.19721.19660.8 1.24931.23051.24751.25001.24930.9 1.30661.28741.30501.30721.30661.0 1.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.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取值范围，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法function vu=Ritz_Galerkin(x0,y0,x1,y1,h)%x0,x1为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、若