实验4微分方程数值解

《实验4微分方程数值解》由会员上传分享，免费在线阅读，更多相关内容在教育资源-天天文库。

1、数学实验常微分方程数值解实验4§1欧拉方法和龙格-库塔(Runge-Kutta)方法设有常微分方程初值问题：常微分方程数值解：计算（精确）解y(x)在一系列离散点x0

2、)=x(1)+n*h;y(n+1)=y(n)+h*feval(f,x(n),y(n));endz=[x'y'];>>z=foeula('ode121',0,1,1,0.1);>>y1=sqrt(1+2*x);>>sol=[zy1’]sol=01.00001.00000.10001.10001.09540.20001.19181.18320.30001.27741.26490.40001.35821.34160.50001.43511.41420.60001.50901.48320.70001.580

3、31.54920.80001.64981.61250.90001.71781.67331.00001.78481.73211欧拉方法1.1向前欧拉公式设初值为y(x0)=y0，yn+1=yn+hf(xn,yn),n=0,1,2,…这就叫向前欧拉公式。其精度为1阶（局部截断误差为O(h2))。p1p2p3p4xyox1x2x3x4functiondy=ode121(x,y)dy=y-2*x/y;return1.3梯形公式梯形公式:kn=f(xn,yn)kn+1=f(xn+1,yn+1)PnPn+1Qn

4、+1knkn+1xyoxnxn+1j(x)yn+1=yn+hky=adveula('ode121',0,1,1,10)y=01.00000.10001.09590.20001.18410.30001.26620.40001.34340.50001.41640.60001.48600.70001.55250.80001.61650.90001.67821.00001.7379functions=adveula(f,a,b,y0,m)x=zeros(1,m+1);y=zeros(1,m+1);h=(b-

5、a)/m;x(1)=a;y(1)=y0;forn=1:mx(n+1)=x(n)+h;k1=feval(f,x(n),y(n));k2=feval(f,x(n+1),y(n)+h*k1);y(n+1)=y(n)+h*(k1+k2)/2;ends=[x'y'];1.4改进的欧拉公式先由向前欧拉公式算出yn+1的预测值，再把它代替梯形公式的右端的yn+1，作为校正，即称为改进的欧拉公式，它可写作例用改进的欧拉公式解方程functiondy=ode121(x,y)dy=y-2*x/y;return改进的欧拉

6、公式可以从几何上加以解释如下：Q(xn+1,yn+hk1)k1k2yxnxn+1Pn(xn,yn)y=f(x)xoPn+1(xn+1,yn+1)以Pn,Q两点处的斜率之平均值为斜率，求得Pn+1欧拉方法可以推广到解高阶方程和方程组，如：向前的欧拉公式为：解线性方程组,可采用向量函数表示.对于高阶方程,需先降阶化为一阶常微分方程组,然后再按上方法求解.如对高阶线性方程常用的方法是:令转换为一阶常微分线性方程组:注意到这里l1,l2为加权系数a(0

7、组合l1,l2，a的选取标准：让精度更高2龙格-库塔方法由二元函数的泰勒展式（推导在教材P72）,得当满足l1+l2=1,2al2=1时，精度达到2阶。显然解不唯一。特别当l1=l2=1/2，a=1时就是改进的欧拉公式。满足l1+l2=1,2al2=1的公式（39）称为二阶龙格-库塔公式。可以证明，当在小区间[xn,xn+1]内只取2点的二阶龙格-库塔公式的最高精度为2阶（O(h3)）。类似地，在小区间选四个点，适当选取加权系数，得到精度为四阶O(h5)的四阶龙格-库塔公式：同样可推广到解微分方程组

8、和高阶方程。3四阶龙格-库塔方法的MATLAB实现先对所解方程或方程组（高阶方程）编写ODE文件（M函数文件，OrdinaryDifferentialEquation）保存在MATLAB设置的搜索路径，然后在解方程的脚本文件中调用：[t,x]=ode23(@f,ts,x0,options)[t,x]=ode45('f',tspan,x0,options)或options用于设置误差限（缺省时相对误差为10-3，绝对误差为10-6），用语句实现为：options=od