数值分析2常微分方程解

常微分方程的数值解法电子科技大学常微分方程的数值解引言简单的数值方法Runge-Kutta方法一阶常微分方程组和高阶方程在高等数学中我们见过以下常微分方程：6.1引言(1)，(2)式称为初值问题，(3)式称为边值问题。在实际应用中还经常需要求解常微分方程组：本章主要研究问题(1)的数值解法，对(2)~(4)只作简单介绍。（其中L为Lipschitz常数）则初值问题（1）存在唯一的连续解。

考虑一阶常微分方程初值问题其中，y=y(x)是未知函数，y(x0)=y0

是初值条件，而f(x,y)是给定的二元函数.由常微分方程理论知，若f(x)在x[a,b]连续且f满足对y的Lipschitz条件：常微分方程的数值解法有单步法和多步法之分：单步法：在计算yn＋1时只用到前一点yn

的值；多步法：计算yn＋1时不仅利用yn，还要利用yn-1,yn-2,…….，一般k步法要用到yn,yn-1,yn-2,…….,yn-k+1。求问题（1）的数值解，就是要寻找解函数在一系列离散节点x1<x2<……<xn<xn+1

上的近似值y1,y

2,…,yn

。为了计算方便，可取xn=x0+nh，(n=0,1,2,…)，h称为步长。6.2简单的数值方法一、欧拉（Euler）方法在x=x0处，用差商代替导数：由得同理，在x=xn

处，用差商代替导数：由得若记则上式可记为此即为求解初值问题的Euler方法，又称显式Euler方法。Pn+1yOxx0x1x2xnP0P1P2Pny=f(x)xn+1Euler方法的几何意义：（Euler折线法）例:用Euler方法求解常微分方程初值问题并将数值解和该问题的解析解比较。解：Euler方法的具体格式：

xn

y(xn) yn

yn-y(xn) 0.0 0 0 00.2 0.1923 0.2000 0.00770.4 0.3448 0.3840 0.03920.6 0.4412 0.5170 0.07580.8 0.4878 0.5824 0.0946 1.0 0.5000 0.5924 0.09241.2 0.4918 0.5705 0.07871.4 0.4730 0.5354 0.0624取h=0.2,xn=nh,(n=0,1,2…,15),f(x,y)=y/x–2y2

计算中取f(0,0)=1.计算结果如下：xn

y(xn) yn

yn-y(xn)1.6 0.4494 0.4972 0.0478 1.8 0.4245 0.4605 0.03592.0 0.4000 0.4268 0.02682.2 0.3767 0.3966 0.0199 2.4 0.3550 0.3698 0.01472.6 0.3351 0.3459 0.01082.8 0.3167 0.3246 0.0079 3.0 0.3000 0.3057 0.0057由表中数据可以看到，微分方程初值问题的数值解和解析解的误差一般在小数点后第二位或第三位小数上，这说明Euler方法的精度是比较差的。O:数值解—:准确解h=0.2;y(1)=0.2;x=h:h:3;forn=1:14xn=x(n);yn=y(n);y(n+1)=yn+h*(yn/xn-2*yn*yn);endx0=0:h:3;y0=x0./(1+x0.^2);plot(x0,y0,x,y,x,y,'o')若直接对y`=f(x,y)在[xn,xn+1]积分，利用数值积分中的左矩形公式：此即为Euler公式。设y(xn)=yn，则得若用右矩形公式：得上式称后退的Euler方法，又称隐式Euler方法。可用迭代法求解：初值：迭代：k=0,1,……因故当hL<1时，迭代法收敛。二、梯形方法由利用梯形求积公式：得上式称梯形方法，是一种隐式方法。用迭代法求解：初值：迭代：k=0,1,……因故当hL/2<1时，迭代法收敛。由以上分析可以看出，隐式方法的计算比显式方法复杂，需要用迭代法求解非线性方程才能得出计算结果。可采用将显式Euler格式与梯形格式结合使用的方法来避免求解非线性方程。记再用梯形格式计算：－－预测－－校正上面两式统称预测－校正法，又称改进的Euler方法。三、单步法的局部截断误差和精度单步法的一般形式为：（

与f有关）显式单步法形式为：整体截断误差：从x0开始，考虑每一步产生的误差，直到xn，则有误差称为数值方法在节点xn处的整体截断误差。但en不易分析和计算，故只考虑从xn到xn+1的局部情况。定义：设y(x)是初值问题(1)的精确解，则称为显式单步法在节点xn+1处的局部截断误差。若存在最大整数p使局部截断误差满足则称显式单步法具有p阶精度或称p阶方法。注：将Tn+1表达式各项在xn处作Taylor展开，可得具体表达式。Euler方法的局部截断误差：故Tn+1=O(h2)，p=1，（设yn=y(xn)）其中称局部截断误差主项。即Euler方法具1阶精度。（设yn=y(xn)）故Tn+1=O(h3)，p=2，梯形方法的局部截断误差：局部截断误差主项为：梯形方法具2阶精度。6.3Runge－Kutta方法一、Runge-Kutta方法的基本思想由Taylor展式Tn+1=O(hp+1)，若提高p，可提高精度。但因……高阶导数计算复杂，故可从另外角度考虑。分析Euler公式及改进的Euler公式：局部截断误差：O(h2)局部截断误差：O(h3)可用f(x,y)在某些点处值的线性组合得yn+1，增加计算f(x,y)的次数可提高阶数。设法计算f(x,y)在某些点上的函数值，然后对这些函数值作线性组合，构造近似计算公式，再把近似公式和解的泰勒展开式相比较，使前面的若干项吻合，从而获得达到一定精度的数值计算公式。Runge-Kutta方法的基本思想：设ci,

λi,

μij为待定常数。上面第一个式子的右端在(xn,yn)作泰勒展开后,按h的幂次作升序排列：再与初值问题的精确解y(x)在点x=xn处的泰勒展开式相比较，使其有尽可能多的项重合。例如，要求就得到p个方程,从而定出参数ci，

i，

ij，再代入K1，K2，…，Kr的表达式，就可得到计算微分方程初值问题的数值计算公式：若

Tn+1=O(hp+1)，则称其为p阶r级R－K方法。上式称为r级Runge－Kutta方法的计算公式。当r=1时，就是Euler方法。要使Runge-Kutta公式具有更高的阶p，就要增加r的值。下面我们只就r=2推导R-K方法。二、二阶Runge-Kutta方法其中c1,c2,

2,

21待定。上式的局部截断误差为：又由利用二元函数的Taylor展开，得代入Tn+1的表达式，得即c1=1-a，

2=

21=1/(2a)要使上式p=2阶，则需方程组解不唯一，可令c2=a

0

，则

满足上述条件的公式都为2阶R-K公式。称中点公式，相当于数值积分的中矩形公式：如取a=½，则c1=

c2=½，

2=

21=1，即为改进Euler公式。若取a=1，则c1=0，c2=1，

2=

21

=½，得例:蛇形曲线的初值问题令f(x,y)=y/x–2y2,取f(0,0)=1,h=0.2,xn=hn,(n=1,2,…,15)2阶龙格-库塔公式计算格式：

k1=yn/xn–2yn

2,k2=(yn+hk1)/(xn+h)–2(yn+hk1)2yn+1=yn+0.5h[k1+k2]x0=0;y0=0;h=.2;x=.2:h:3;k1=1;k2=(y0+h*k1)/x(1)-2*(y0+h*k1)^2;y(1)=y0+.5*h*(k1+k2);forn=1:14k1=y(n)/x(n)-2*y(n)^2;k2=(y(n)+h*k1)/x(n+1)-2*(y(n)+h*k1)^2;y(n+1)=y(n)+0.5*h*(k1+k2);endy1=x./(1+x.^2);subplot(221)plot(x,y1,x,y1,'b*')subplot(222)plot(x,y,x,y,'o')subplot(223)plot(x,y,x,y,'o',x,y1,x,y1,'b*')

三、三阶与四阶Runge-Kutta方法当r=3时，R-K公式表示为其中为8个待定常数。上式的局部截断误差为类似二阶的推导过程，将K2,K3按二元函数展开，使Tn+1=O(h4)，得方程有8个未知数，解不唯一。满足该条件的公式统称为三阶R-K公式。其中一个常用公式为：当r=4时，利用相同的推导过程，经过较复杂的计算，可以得出四阶R-K公式的成立条件。下列经典公式是其中常用的一个：k1=f(xn,yn),k2=f(xn+0.5h,yn+0.5hk1)k3=f(xn+0.5h,yn+0.5hk2),k4=f(xn+h,yn+hk3)

yn+1=yn+h[k1+2k2+2k3+k4]/6(n=0,1,2,·····,N)四阶龙格-库塔公式数值求解:狐狸和野兔问题——常微分方程组在一个生物圈中,只有野兔和狐狸两种动物,记y1

为野兔数量,y2

为狐狸数量.设有足够的食物供野兔生存,而狐狸只以野兔为食物.

模型如下自变量x∈(0,15),初值条件为:y1(0)=20,y2(0)=20t0=0;t1=15; %设置区域y0=[2020]; %给定初值[t,y]=ode23(‘lot1’,t0,t1,y0); %求解plot(t,y) %绘图functionz=lot1(x,y)z(1,:)=y(1)-0.01*y(1).*y(2);z(2,:)=-y(2)+0.02*y(1).*y(2);求解常微分方程:(1)定义函数;(2)调用ode236.5一阶常微分方程组和高阶微分方程的数值解法简介一、一阶常微分方程组的数值解法：下