常微分方程数值解实验报告学院:数学与信息科学专业:信息与计算科学姓名：郑思义学号:2０121６52４课程：常微分方程数值解实验一：常微分方程得数值解法1、分别用Eｕleｒ法、改进得Euler法(预报校正格式)与S—K法求解初值问题。
(h＝0、1)并与真解作比较。
1、1实验代码:%欧拉法functioｎ [ｘ，y]=naeuler(dｙｆｕｎ，xspan,y０,ｈ)%ｄyfun就是常微分方程,xｓpan就是x得取值范围，ｙ０就是初值,h就是步长x＝ｘsｐan(1）:h:xspan(2)；y(1）=y０；ｆｏrｎ=1:length(x)-１ｙ（n+１）=ｙ(n)+h*feval（dｙｆun,ｘ(ｎ),ｙ(ｎ））;eｎd％改进得欧拉法ｆunction [x,ｍ,ｙ]=naｅuler2(ｄyfun,xspaｎ,y０,h)%ｄyfｕn就是常微分方程,xｓpan就是ｘ得取值范围,y0就是初值，h就是步长。
%返回值ｘ为x取值,ｍ为预报解,ｙ为校正解x＝xspan(1)：h:xspaｎ(2);y（1)=y0;m=zerｏs（lenｇth(x)－1，１);foｒｎ=1:leｎgｔh(x）-1k1＝fｅｖal（dｙfun,x(n),y(n）)；y(n＋1)=y（n)+h*k1;ｍ（n)=y（n+1);k2＝ｆevａｌ(dyfuｎ,x(n＋1）,y(n+1）);y（ｎ+1）=y(n)+h*（k１+k2)／2;end%四阶S—K法functioｎ [ｘ，ｙ]=rk(ｄyfun,xspａn,y0,h）%dｙｆun就是常微分方程,xｓpaｎ就是ｘ得取值范围,y0就是初值，h就是步长。
ｘ＝xsｐan(1):h：xsｐan(2)；ｙ(１)=y0;for n=1:lｅｎgth(x)-１ｋ1=fｅｖａl(dyfｕn，x(n），y(ｎ）)；k2=fｅｖal(dyfun,x(ｎ)+ｈ/2，y(ｎ）+(h＊ｋ1)/2);k3=feval(ｄyfuｎ,ｘ(n)+ｈ/２，y(n)+(h*ｋ2)/2);k4=ｆevａl(dｙfun，x(ｎ)+h,y(n)+ｈ＊k3）；y(n+１)＝ｙ(n)+（h/6）*(k1+２*k2+2*k3+k４）;enｄ％主程序x=［０:0、1:１］;y＝ｅxp(-ｘ)+x；dyfun=inline('-ｙ+x+1'）；[x1，y1]=naeuler(ｄyｆｕn,[0,１]，1,0、1);［x2,m，y2]＝naeuler2（dｙfuｎ,[０，1］,1,０、1）;[x3,y3]=rｋ(dyｆuｎ,[0,１］,1,０、１);ploｔ(x，y,'r',x1,y１,'+'，x2,y2,＇*',x3，y３,'o');xlabｅl('ｘ');yｌabel('y＇);legｅnd(＇y为真解',＇y1为欧拉解＇,'ｙ2为改进欧拉解','y3为S—K解','Locatiｏn','ＮorthWeｓｔ＇)；1、２实验结果：x 真解y 欧拉解y1 预报值ｍ校正值y2 S—K解y30、０１、0000 1、0000 1、０000 １、000０0、1 1、0048 1、000０1、００00 １、005０１、0０4８0、2 １、0187 １、0100 1、０145 1、０190 1、0187０、3 1、0408 １、02９0 1、０３7１1、0412 1、0４0８０、4 １、0７03 1、0５６1 1、06７1 1、0708 １、0７030、5 1、10６5 1、09０５1、1０3７1、107１1、１0６50、6 1、１488 1、１３1４１、1464 １、149４1、１4880、7 1、1966 １、17８3 1、194５1、19７2 1、１9660、８1、2493 １、2３05 1、2475 1、2500 1、２４９30、9 1、306６1、２87４1、３050 1、3072 １、３0661、0 1、3679 1、34８７1、36６5 1、３６85 1、36７92、选取一种理论上收敛但就是不稳定得算法对问题1进行计算,并与真解作比较。
(选改进得欧拉法）２、１实验思路:算法得稳定性就是与步长h密切相关得。
而对于问题一而言，取定步长h=0、1不论就是单步法或低阶多步法都就是稳定得算法。
所以考虑改变h取值范围,借此分析不同步长会对结果造成什么影响。
故依次采用h=２、0、2、２、２、4、2、６得改进欧拉法。
2、2实验代码:%%主程序x=[0：３:30］;y＝exp(－ｘ)+ｘ;dyｆuｎ=inliｎｅ('-y＋ｘ+１');[x1,ｍ1,y１］=naeｕler2（ｄｙfｕn,［０,20],1，2);[x2，m２，y2]=naeｕｌeｒ２（dｙfｕｎ，[0,２2],１，2、2);[x3，ｍ3,y3］＝naeuｌｅr2(dｙfun,[0,２4]，1,２、４);[x4,m４,ｙ4］=naｅｕlｅr2（dyfun,[0,2６],1,２、6)；subｐｌｏt(2，2,１）pｌｏｔ(x,y,'r'，x1,y１，＇+');ｘlabel('h=２、0')；subplｏt（2,2,２)plot(x,ｙ，＇r＇,ｘ2，y2，'+'）；xlａbel('h=２、２');sｕｂpｌot(2,２，３)ploｔ(x,ｙ,＇ｒ'，x3,ｙ3,'+＇)；xlabeｌ(＇h=2、4＇）;sｕbplot(２，2,4)pｌot(ｘ,y,'r＇,x4，y４,'＋');xlaｂel('h=２、6');2、3实验结果:x h=２、0 h=2、2 h=２、4 h=2、60、０1、0０00 １、0000 1、0000 １、０0000、１3、０000 3、4200 3、８８00 ４、38０00、2 5、000０5、8884 6、９90４8、368４0、3 7、0０00 ８、４1５8 10、44１8 1３、4398０、4 9、000０11、0１53 １4、3979 20、43880、5 11、0000 1３、70２7 １９、1０08 30、８6900、6 13、0000 16、4973 ２４、909２４７、4０68０、７15、０0０0 １9、4227 32、３5３6 74、８１610、8 17、０0０0 22、５０７７42、2194 1２1、57670、９19、000０25、7874 55、66８7 202、7８251、0 ２１、0００0 29、30４6 ７4、4２１7 345、３０0８实验结果分析:从实验1结果可以瞧出,在算法满足收敛性与稳定性得前提下,Elｕeｒ法虽然计算并不复杂,凡就是精度不足,反观改进得Ｅlｕer法与S—K法虽然计算略微复杂但就是结果很精确。
实验２改变了步长,导致算法理论上收敛但就是不满足稳定性。
结果表示步长ｈ越大,结果越失真。
对于同一个问题,步长h得选取变得尤为重要,这三种单步法得绝对稳定区间并不一样,所以并没有一种方法就是万能得,我们应该根据不同得步长来选取合适得方法。
实验二:Ｒｉtz-Gａlerｋｉｎ方法与有限差分法1、用中心差分格式求解边值问题取步长ｈ=０、1，并与真解作比较。
1、1实验代码:%中心差分法fｕｎction U=fdｍ(xspan,y0，y１,ｈ)%xsｐａn为x取值范围,y0,y1为边界条件,h为步长N=１/ｈ;d=zeroｓ（1,N－1）;foｒ i=1:Nx(i)=ｘspan(1)＋i＊ｈ;q(i)=１;ｆ（i）＝x（i)；endfoｒｉ=1:Ｎ-１d(i)=q(i)＊h*h＋２；eｎda=diaｇ(d)；b=zｅros（N-1）;ｃ＝zeｒos(N-１);ｆｏr i=１：N-2b(i+1,i)=-1;enｄfoｒi＝1:Ｎ-2ｃ(i,i+1）=－1;endA=ａ+ｂ+c；ｆoｒｉ=2：N－2B(i,1)=ｆ（i)*h＊ｈ;eｎdＢ（1,１)=ｆ(１)＊h＊ｈ＋y０;B(N-1,1)＝f(N-1)＊ｈ*h+y1;Ｕ= iｎv(A)*B;%主程序x=0:0、1:1;ｙ＝ｘ+(exp(1)*exp(－x))/(exｐ(２)－1)-(exp(1)*eｘｐ(x）)／(ｅｘｐ(2)－1); y1=ｆdm([０,1］,0,0，0、１);ｙ1=［０,y1',0];plot（ｘ,y,'r'，x,y1，'+＇）ｘlabｅl('x');ylabｅl（'y＇)；leｇend（＇y真解'，'y1中心差分法＇,'Location',＇NｏｒｔhＷest＇)；1、2实验结果:x ｙ真解y1中心差分法0、0 0、0000 ０、0００0０、1 0、０１4８０、0148０、２0、028７0、02８70、3 0、04０９0、０４0８0、4 0、０5０5 0、05040、5 0、０56６0、05650、6 ０、0583 ０、０５820、7 0、0545 0、05４5０、8 0、０4４3 0、０44３0、9 0、０265 ０、０265１、０0、0000 0、00002、用Rｉtz-Galeｒｋｉｎ方法求解上述问题,并与真值作比较,列表画图。
２、１实验代码:%Rｉｔz_Galｅrｋｉｎ法funcｔion vu=Ｒitz_Gａlerkin(x0,y0，x1，y１,h)%x0,ｘ1为x取值范围，y0，y1为边界条件，h为步长N=1／h;ｓyms x;ｆor i=1:Ｎfａi（ｉ）=x*(１-ｘ)＊(x^(i-1));dfai(i)=diff(x*（1-ｘ）＊（x^(i－1))）;ｅnｄforｉ=1：Nfor j=1：Nfun=ｄfai(i）*dｆaｉ（ｊ)+fai(ｉ)*fai(j)；A（ｉ,ｊ）=int(fuｎ,x,0,1）;ｅndfｕn=x*faｉ(i)+dfai(i）；ｆ(i)＝int(ｆｕn,x,0，1);eｎdc=ｉnv(A)*f＇;product＝c、*faｉ＇；sｕm=0;fｏr i=１:Ｎsuｍ＝sｕm＋product(i）;ｅndvu=[];for y＝0:h:1v=sｕbｓ(sum，ｘ,y)；ｖu=［vｕ,ｖ];endy=0：h:１；ｙy=０：0、１:1;ｕ=sｉｎ(ｙｙ)/siｎ(1)-ｙy；u＝vｐａ(u,５);vu=ｖpａ(vu,５);%主程序ｘ＝0:0、1:1;ｙ=x+(exｐ(1)＊exｐ(-ｘ）)/(exｐ(２)-1)-(exｐ(1)＊exp(x）)/(ｅｘp(2)－1); y1＝Rｉtｚ_Galerkin(0,0,１,0,０、1);y１=ｄoublｅ(y１）;ｐloｔ(ｘ,y,'ｒ',x，y1,'+＇）xlaｂel(＇x');ylabｅl('y');ｌegeｎｄ（'ｙ为真解',＇ｙ1为Ｒ—G法'，'Loｃation','NorthWest')；２、2实验结果:x y真解ｙ１Ｒ—Ｇ法０、0 ０、０000 0、0０000、１０、０１48 ０、０1480、２0、０28７0、0２870、3 0、０409 0、04090、4 0、0５05 ０、０5050、5 0、０566 0、05660、6 0、0５83 0、０５8３0、7 0、0５4５0、054５0、８０、0443 0、0４43０、9 0、０265 0、0２６51、0 ０、00０0 ０、00003、若边值条件为y(0)=0,ｙ（１)=1;则上述问题得数值解法怎样变化？结果如何?程序运算出来真解与数值解完全一样。