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

常微分方程数值解实验报告 学院:数学与信息科学 专业:信息与计算科学 姓名：郑思义 学号: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）-1 k1＝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真...