常微分方程模型及其数值解.pptVIP

改进的Euler方法可写为构造差分公式…………其中??i,?i,?ij?为待定参数.若此公式的局部截断误差为由于yn+1=yn+h?1?n+h?2(?n+?h?xn+?h?n?yn)+O(h3)O(h3),称此公式为p阶Runge-kutta方法,简称p阶R-K方法.对于p=2的情形,应有=yn+h(?1+?2)?n+h2?2(??xn+??n?yn)+O(h3)所以,只要令?1+?2=1,??2=1/2,??2=1/2(8)一般地,参数由(8)确定的一族差分公式(7)统称为二阶R-K方法.称之为中点公式,或可写为若取?=1,则得?1=?2=1/2,?=1,此时公式(7)就是改进的Euler公式;若取?1=0,则得?2=1,?=?=1/2,公式(7)为高阶R-K公式可类似推导.下面列出常用的三阶、四阶R-K公式.四阶标准R-K公式三阶R-K公式解四阶标准R-K公式为例3用四阶标准R-K方法求初值问题y?=y-2x/y,0?x?1y(0)=1的数值解,取步长h=0.2.计算结果如下:MyRunge_Kutta.mfunction[x,y]=MyRunge_Kutta(fun,x0,xt,y0,PointNum,varargin)%Runge-Kutta方法解微分方程形为y’(t)=f(x,y(x))%此程序可解高阶的微分方程。只要将其形式写为上述微分方程的向量形式%x范围为[x0,xt],初值为y0,PointNum为离散点数,varargin为可选输入项可传适当参数给函数f(x,y)ifnargin4|PointNum=0PointNum=100;endifnargin3y0=0;endy(1,:)=y0(:);%初值存为行向量形式h=(xt-x0)/(PointNum-1);%计算步长x=x0+[0:PointNum]*h;%得x向量值fork=1:PointNum%迭代计算f1=h*feval(fun,x(k),y(k,:),varargin{:});f1=f1(:);%得公式中k1f2=h*feval(fun,x(k)+h/2,y(k,:)+f1/2,varargin{:});f2=f2(:);%得公式中k2f3=h*feval(fun,x(k)+h/2,y(k,:)+f2/2,varargin{:});f3=f3(:);%得公式中k3f4=h*feval(fun,x(k)+h,y(k,:)+f3,varargin{:});f4=f4(:);%得公式中k4y(k+1,:)=y(k,:)+(f1+2*(f2+f3)+f4)/6;%得y(n+1)endfunctionex03()%用三种不同方法解微分方程clear,clfx0=0;xt=1;Num=100;h=(xt-x0)/(Num-1);x=x0+[0:Num]*h;a=1;yt=sqrt(1+2*x);%真值解y0=1;%设定函数初值PointNum=4;%设定取点数[x1,ye]=MyEuler(myfun012,x0,xt,y0,PointNum);[x2,yh]=MyEulerPro(myfun012,x0,xt,y0,PointNum);[x3,yr]=MyRunge_Kutta(myfun012,x0,xt,y0,PointNum);plot(x,yt,k,x1,ye,b:,x2,yh,g:,x3,yr,r:)legend(真值,简单欧拉法解,改进欧拉法解,Runge-Kutta法解)holdonplot(x1,ye,bo,x2,yh,b+,x3,yr,r*)ex03.m常微分方程数值解杨瑞琰在许多实际问题中，例如物理中的速率问题，人口的增长问题，放射性衰变问题，经济学中的边际问题等，常常涉及到两个变量之间的变化规律。微分方程是研究上述问题的一种机理分析方法，它在科技、工程、生态、环境、人口以及经济管理等领域中有着十分广泛的应用。01在应用微分方程解决实际问题时，必须经过两个阶段。一是微分方程的建立，建立一个微分方程的实质就是构建函数、自变量以及函数对自变量的