常微分方程數(shù)值解實(shí)驗(yàn)報(bào)告

1、.常微分方程數(shù)值解實(shí)驗(yàn)報(bào)告 學(xué)院：數(shù)學(xué)與信息科學(xué) 專業(yè)：信息與計(jì)算科學(xué)：鄭思義 學(xué)號(hào)：201216524 課程：常微分方程數(shù)值解實(shí)驗(yàn)一：常微分方程的數(shù)值解法1、 分別用Euler法、改進(jìn)的Euler法（預(yù)報(bào)校正格式）和SK法求解初值問題。（h=0.1）并與真解作比較。1.1實(shí)驗(yàn)代碼：%歐拉法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

2、),y(n);end%改進(jìn)的歐拉法function x,m,y=naeuler2(dyfun,xspan,y0,h)%dyfun是常微分方程，xspan是x的取值范圍，y0是初值，h是步長。%返回值x為x取值，m為預(yù)報(bào)解，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;en

3、d%四階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+

4、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為真解','

5、y1為歐拉解','y2為改進(jìn)歐拉解','y3為SK解','Location','NorthWest');1.2實(shí)驗(yàn)結(jié)果：x真解y歐拉解y1預(yù)報(bào)值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

6、.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、 選取一種理論上收斂但是不穩(wěn)定的算法對問題1進(jìn)行計(jì)算，并與真解作比較。（選改進(jìn)的歐拉法）2.1實(shí)驗(yàn)思路：算法的穩(wěn)定性是與步長h密切相關(guān)的。而對于問題一而言，取定步長h=0.1不論是單步法或低階多步法都是穩(wěn)定的算法。所以考慮

7、改變h取值范圍，借此分析不同步長會(huì)對結(jié)果造成什么影響。故依次采用h=2.0、2.2、2.4、2.6的改進(jìn)歐拉法。2.2實(shí)驗(yàn)代碼：%主程序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

8、,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實(shí)驗(yàn)結(jié)果：xh=2.0h=2.2h=2.4h

9、=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.

10、2194 121.5767 0.9 19.0000 25.7874 55.6687 202.7825 1.0 21.0000 29.3046 74.4217 345.3008 實(shí)驗(yàn)結(jié)果分析：從實(shí)驗(yàn)1結(jié)果可以看出，在算法滿足收斂性和穩(wěn)定性的前提下，Eluer法雖然計(jì)算并不復(fù)雜，凡是精度不足，反觀改進(jìn)的Eluer法和SK法雖然計(jì)算略微復(fù)雜但是結(jié)果很精確。實(shí)驗(yàn)2改變了步長，導(dǎo)致算法理論上收斂但是不滿足穩(wěn)定性。結(jié)果表示步長h越大，結(jié)果越失真。對于同一個(gè)問題，步長h的選取變得尤為重要，這三種單步法的絕對穩(wěn)定區(qū)間并不一樣，所以并沒有一種方法是萬能的，我們應(yīng)該根據(jù)不同的步長來選取合適的方法。實(shí)驗(yàn)二：Ritz

11、-Galerkin方法與有限差分法1、 用中心差分格式求解邊值問題取步長h=0.1，并與真解作比較。1.1實(shí)驗(yàn)代碼：%中心差分法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

12、 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');lege

13、nd('y真解','y1中心差分法','Location','NorthWest');1.2實(shí)驗(yàn)結(jié)果：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-

14、Galerkin方法求解上述問題，并與真值作比較，列表畫圖。2.1實(shí)驗(yàn)代碼：%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

15、(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

16、(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實(shí)驗(yàn)結(jié)果：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.