有限差分的欧拉法.docVIP

西北农林科技大学实习报告 学院：理学院 专业年级：信计061 姓名：袁金龙 学号课程：微分方程数值解 报告日期：2008-11-26 实习一、一维问题的有限差分方法-----Euler法 一）实习问题： 用欧拉法，龙格库塔法，米尔恩法求解下面的初值问题： 二）算法描述： ⑴欧拉法： ⑵龙格库塔法的中点法： ⑶米尔恩法： ⑷初始值的确定： 泰勒级数法 ⑸雅可比迭代法 三）matlab程序： ⑴问题函数： function [f1]=f(t,u) f1=-2*u+exp(t); *************************************************** ⑵欧拉法： function []=oula(a,b,u0,n) %[a,b]表示t的取值区间 %u0表示初值 %n表示将[0，1]区间分成的分数 h=(b-a)/n; t0=a; u(1)=u0+h*(f(t0,u0)); for i=1:n t(i)=a+i*h; end t for i=2:n u(i)=u(i-1)+h*f(t(i-1),u(i-1)); end u %精确解的求法 for i=1:n u1(i)=(2/3)*exp(-2*t(i))+(1/3)*exp(t(i)); end u1 plot(t,u,t,u1) title(欧拉法中的预测值与真实值的比较); xlabel(采样点);ylabel(幅度);grid; legend(预测值,真实值); *********************************************************** ⑶龙格库塔法的中点法： function []=zhongdian(a,b,u0,n) %[a,b]表示t的取值区间 %u0表示初值 %n表示将[0，1]区间分成的分数 h=(b-a)/n; t0=a; u(1)=u0+h*(f((t0+(1/2)*h),(u0+(1/2)*h*f(t0,u0)))); for i=1:n t(i)=a+i*h; end t for i=2:n u(i)=u(i-1)+h*(f((t(i-1)+(1/2)*h),(u(i-1)+(1/2)*h*f(t(i-1),u(i-1))))); end u %精确解的求法 for i=1:n u1(i)=(2/3)*exp(-2*t(i))+(1/3)*exp(t(i)); end u1 plot(t,u,t,u1) title(中点法中的预测值与真实值的比较); xlabel(采样点);ylabel(幅度);grid; legend(预测值,真实值); ************************************************************** ⑷米尔恩法： function []=Milne(n,e) %n表示区间[0,1]要分的份数 %对于初值精度我们选择q=3阶的 %初值的求解 %经过计算运用泰勒级数法求解初值 %q(i)表示第i次迭代满足精度的最终值 t0=0; for i=1:n t(i)=(1/n)*i; end u0=1; h=1/n; u2(1,1)=1-h+2.5*h*h %Jacobi迭代求解隐式法 u2(2,1)=u2(1,1)+h*f(t(1),u2(1,1)) for i=2:1000000 u2(2,i)=h*(1/3)*f(t(2),u2(2,i-1))+u0+h*(1/3)*f(t0,u0)+h*(4/3)*f(t(1),u2(1,1)); if abs(u2(2,i)-u2(2,(i-1)))e u2(3,1)=u2(2,i)+h*f(t(2),u2(2,i)); break; end end for i=3:n for j=2:1000000 u2(i,j)=h*(1/3)*f(t(i),u2(i,j-1))+u2(i-2,1)+h*(1/3)*f(t(i-2),u2(i-2,1))+h*(4/3)*f(t(i-1),u2(i-1,1)); if abs(u2(2,j)-u2(2,(j-1)))e u2(i+1,1)=u2(i,j)+h*f(t(i),u2(i,j)); break; end end end q(1)=u2(1