数值分析(常微分方程数值解)

解一阶常微分方程欧拉法局部截断误差与p阶精度Range-Kutta公式一阶常微分方程组和二阶方程线性多步法简介《数值分析》23*常微分方程(ODE)的定解问题主要有初值问题(IVP)和边值问题(BVP)两大类。

急性传染病数学模型(SIR模型)常微分方程的初值问题是描述系统发展演变的重要工具。洛伦兹吸引子

lorenzgui关于爱情的动力系统

Loveraffairesanddifferentialequations*一阶常微分方程初值问题:离散化(discretization)

取定离散点x0<x1<x2<···<xN,步长h

其中y=y(x)是未知函数,右端函数f(x,y)是已知函数,初值y0是已知数据。求未知函数y(x)在离散点处的近似值y1,y2,y3,·····,yN考察点xn列出微分方程Euler法与修正的Euler法其中yn是y(xn)的近似。足够多的点来近似连续对象*求解常微分方程初值问题的Euler方法

取定步长:h,记

xn=x0+nh,(n=1,2,···,N)称计算格式:yn+1=yn+hf(xn,yn)为显式Euler公式。对应的求初值问题数值解的方法称为Euler方法。例1.

用Euler法求初值问题的数值解。解:记

f(x,y)=y－xy2,xn=nh(n=0,1,2,···,N)

由Euler公式得:yn+1=yn+h(yn－xnyn2)(n=0,1,···,N)*取步长

h=2/10,2/20,2/30,2/40,用Euler法求解的数值实验结果如下.N10203040h0.20.10.06670.05误差0.10590.05210.03420.0256解析解:o——数值解----——准确解

*y'=f(x,y)左矩形公式用数值积分方法离散化常微分方程显式Euler法求解常微分方程初值问题的数值方法可以由很多,各种方法都可以看成使用不同的数值积分方法计算*梯形公式:

右矩阵公式:

隐式Euler法*预-校方法算法如下k1=f(xn,yn),

k2=f(xn+1,yn+hk1)

,预报-校正方法(修正的Euler法):预报校正*预-校方法(h=0.2)误差最大值0.0123n10203040h0.20.10.06670.05误差20.01230.00260.00115.9612e-004误差10.10590.05210.03420.0256欧拉方法(h=0.2)误差最大值0.1059*

基本概念单步法:即为了求得后一点xn+1上数值解yn+1,只要知道一点xn的数值解yn就可以了。多步法:即计算yn+1时,要用到前面多点的数值解yn,yn-1,···,yn-r等信息。显式格式:数值解yn+1可以用xn和yn

解析表出。这类格式称为显式格式。反之称为隐式格式。对隐式格式而言,计算数值解yn+1形式上需要解方程(组)。显式方法易于计算,但通常数值稳定性较差。而隐式方法通常稳定性较好,但计算上显然很不方便。将显式方法和隐式方法联合使用是重要的思想方法。*

Tn+1=y(xn+1)-yn+1称为局部截断误差。设

yn=y(xn)即由Taylor公式Euler公式:yn+1=yn+hf(xn,yn)的局部截断误差Euler公式的局部截断误差记为

O(h2),则称Euler公式具有1阶精度。局部截断误差仅考虑xn到xn+1的局部情况,并假定xn之前的计算没有误差*若局部截断误差为O(hp

+1)则称显式单步法具有

p阶精度。例2.

证明梯形方法具有2阶精度从x0开始计算,如果考虑每一步产生的误差直到xn,则误差y(xn)-yn称为数值方法在节点xn处的整体误差。粗略的讲,整体误差把许多步的局部截断误差加在一起,而步数正比于1/h(即步长倒数),因此整体误差O(hp)。*

(Runge-Kutta)龙格-库塔法改进的Euler方法这样理解:它用xn和xn+1两个点的斜率值K1和K2取算术平均作为平均斜率,而xn+1处的斜率值则利用已知信息xn来预报。RK方法是一大类的方法,其基本思想是采用如下形式:是否可以推广改进的Euler方法？*取点个数s称为龙格-库塔法的级数。一般的RK公式如下给出:

RK方法是单步法(为了求得后一点xn+1上数值解yn+1,只要知道一点xn的数值解yn)。我记得我的朋友JohnvonNeumann曾经说过,用四个参数我可以拟合出一头大象,而用五个参数我可以让它的鼻子摆动-费米下面以2级方法为例子具体介绍龙格-库塔法。构造的基本思想是选择适当的系数使得方法的局部截断误差阶数尽可能高。**三阶Range-Kutta公式一般形式yn+1=yn+h[k1+4k2+k3]/6k1=f(xn,yn),k2=f(xn+0.5h,yn+0.5hk1)k3=f(xn+h,yn–hk1+2hk2)四阶Range-Kutta公式一般形式yn+1=yn+h[k1+2k2+2k3+k4]/6k1=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)*例4数值实验:几种不同求数值解公式的误差比较n10203040h0.20.10.06670.05RK46.862e-0053.747e-0067.071e-0072.186e-007RK30.00121.529e-0044.517e-0051.906e-005RK20.01230.00260.00115.9612e-004Euler0.10590.05210.03420.0256*常微分方程初值问题的MATLAB数值求解命令:[x,y]=ode23('f',[a,b],y0)f=inline('y-x.*y.^2');[x,y]=ode23(f,[0,2],1)符号求解命令:dsolve('eqn1',...)symsxydsolve('Dy=y-x*y^2','y(0)=1','x')解析解:参考:TheMatlabODEsuite,SIAMJScientificComputing*一阶常微分方程组欧拉公式:(n=0,1,2,·······,N)一阶常微分方程组向量形式记*修改的欧拉公式:(n=0,1,·······,N)经典龙格-库塔公式:*洛伦兹模型取=8/3，=10，=28。x(0)=0，y(0)=0，z(0)=0.01。t∈[0，80]，微分方程右端函数: