计算方法引论-第十一章

计算方法引论：

微分方程数值解法常微分方程初值问题的数值解法双曲型方程的差分解法抛物型方程的差分解法橢圆型方程的差分解法有限元方法计算方法引论(第三版)第十一章常微分方程初值问题数值解法几种简单的数值解法R-K方法线性多步法预估—校正公式方程组和高阶方程的数值解法自动选取步长的需要和事后估计刚性方程计算方法引论(第三版)引言问题求(11.1)的解y(x)在存在区间[a,b]中点列xn=xn–1+hn

(n=1,…,N)上的近似值yn. (11.1)hn＞0,步长.为叙述方便，假设hn不变,记为h.假设解y(x)在区间[a,b]上存在、唯一并充分光滑.

f(x,y)也充分光滑.

解法一步法：计算yn+1时只用到前一步的值.

多步法：计算yn+1要用到前k步的值.计算方法引论(第三版)

Euler方法

计算公式

(11.2)

由来几何解释过(xn,yn)沿斜率f(xn,yn)方向行一步

Taylor展开

数值微分

数值积分左矩形公式用于计算方法引论(第三版)

Euler方法算例

计算问题

y′=y-2x/y,0≤x≤1,y(0)=1，h=0.1计算结果

计算方法引论(第三版)Euler方法的误差估计局部截断误差准确解代入

(11.2)式的剩余总体截断误差

,≤M2h2/2

计算方法引论(第三版)Euler方法的绝对稳定性绝对稳定性概念

用一个数值方法求解微分方程y′=λ

y,其中,λ是一个复常数,对给定步长h>0,在计算yn时引入了误差若这个误差在后面的计算中所引起的误差绝对值均不增加,就说这个数值方法对于这个步长h和复数,λ是绝对稳定的.Euler方法绝对稳定性误差方程:绝对稳定条件:绝对稳定区域:复平面上以λh=

-1为中心，1为半径的圆

计算方法引论(第三版)

向后Euler方法

计算公式(11.10)

由数值积分右矩形公式导出误差局部截断误差总体截断误差O(h)隐式方法迭代求解向后Euler方法是A稳定的误差方程:A稳定:绝对稳定区域为整个左半平面,Re(λh)<0

计算方法引论(第三版)

梯形方法

计算公式(11.14)

由数值积分梯形求积公式导出误差局部截断误差总体截断误差O(h2)隐式方法迭代求解梯形方法是A稳定的误差方程:计算方法引论(第三版) 改进Euler方法

计算公式(11.17)

可视为一种预估—校正方法先用Euler方法预估pn+1=yn+hf(xn,yn)再用梯形方法校正yn+1=yn+h/2(f(xn,yn)+f(xn+1,pn+1))局部截断误差O(h3)总体截断误差O(h2)p阶方法

总体截断误差为O(hp)的方法称p阶方法方法的总体截断误差为O(hp)则局部截断误差为O(hp+1)计算方法引论(第三版) 改进Euler方法算例

计算问题

y′=y-2x/y,0≤x≤1,y(0)=1，h=0.1计算结果

计算方法引论(第三版)

Taylor展开法一阶方法到二次项得Euler方法yn+1=yn+hf(xn,yn)

r阶方法到r+1次项得y′(x)=f(x,y(x))计算量很大.可作几个点上f值的线性组合来构造近似公式使其有r+1阶局部截断误差,乃得R-K方法计算方法引论(第三版)几个显式R-K方法二阶方法三阶方法四阶(经典R-K)方法绝对稳定区域误差方程绝对稳定区域

计算方法引论(第三版)

经典R-K方法算例

计算问题

y′=y-2x/y,0≤x≤1,y(0)=1，h=0.2计算结果

计算方法引论(第三版)

Adams外插公式计算公式

(11.28)不能“自开始”可用四阶R-K方法计算y1，y2，y3.由来出发点令F(x)=f(x,y(x))以xn,xn-1,xn-2,xn-3为节点作插值求积公式局部截断误差

计算方法引论(第三版)

Adams内插公式计算公式

(11.30)不能“自开始”.

是隐式方法要迭代求解..由来出发点令F(x)=f(x,y(x))以xn+1,

xn,xn-1,xn为节点作插值求积公式局部截断误差

计算方法引论(第三版)Adams方法的稳定性Adams外插法的绝对稳定性误差方程:通解：特征方程：绝对稳定条件:绝对稳定区域计算方法引论(第三版)Adams方法的稳定性(续)Adams内插法的绝对稳定性误差方程:通解：特征方程：绝对稳定条件:绝对稳定区域计算方法引论(第三版)多步法稳定性的注记多步法的绝对稳定性本章从某步误差对后面计算的影响考虑提出绝对稳定性概念并讨论了各方法的绝对稳定性及绝对稳定区域.一般多步法特征方程的根结构复杂,绝对稳定性也可定义为特征方程的根按模小于1(ρn→0,n→∞).稳定性与收敛性尙有一些稳定性概念以处理不同问题.类似于十二、三章有一种稳定性保证(Dahlquist定理):

相容的差分方法(局部截断误差为:o(h))收敛的充要条件是其稳定,即特征方程h=0时的根在单位圆,圆周上只有单根.因此本章各方法收敛计算方法引论(第三版)预估-校正公式显式方法计算初值隐式方法进行迭代

Ｅuler-梯形Adams只迭代一次预估

校正

预估

校正

计算方法引论(第三版) 算例:R-K与AdamsP-C法计算问题

y′=y-2x/y,0≤x≤1,y(0)=1，h=0.1计算结果

计算方法引论(第三版)误差估计Adams方法局部截断误差的估计由于得到估计据此可改善结果,乃有Adams预估—校正法的修正方案.

计算方法引论(第三版)修正的Adams预估-校正方法Adams预估—校正方法的一种修正方案原误差的主要部分得以补偿,是五阶方法.y1，y2，y3另用单步法计算.并令

p3=c3=0.计算方法引论(第三版)常微分方程组解法常微分方程组n=2Euler方法向后Euler方法梯形公式k=0,1,…Adams外插公式计算方法引论(第三版)常微分方程组解法(续)常微分方程组n=2Adams预估—校正方法k=0,1,…R-K方法计算方法引论(第三版)高阶常微分方程解法化成方程组的形式求解二阶方程化为一阶方程组函数的特定形式计算公式可以化简例如R-K公式

计算方法引论(第三版)事后估计误差事后估计误差依据是误差的先验估计方法一:利用同一公式不同步长的误差来估计计算结果的误差.数值积分中逐次分半法和Romberg求积用过.方法二:利用同阶的不同公式的误差来估计计算结果的误差.本章四阶Adams外插和内插公式存在不同阶的二公式,计算结果的差可估计其中一个的误差.改善结果事后误差估计可改善结果.Romberg求积公式Adams预估—校正法的修正方案事后误差估计可调整步长.以一步r阶公式为例介绍一种方法.计算方法引论(第三版)自动选取步长事后估计误差局部截断误差为cnhr+1从xn–1出发两步(步长为h)计算得从xn–1出发一步(步长为2h)计算得误差可估计误差并改善结果调整步长一种策略:如果

成立步长满足要求.如过小,只1/2r容许误差,步长加倍.否则步长减小一半自动变步长的R-K方法计算例子,效果显著.尚有其它策略.多步法也可设计策略调整步长.不过调整代价较大.计算方法引论(第三版)

自动选取步长R-K方法算例

计算问题

其中

计算方法引论(第三版)

自动选取步长R-K方法算例

计算结果计算方法引论(第三版)

自动选取步长R-K方法算例

说明问题的解是以T=6.192169为周期的轨迹开始计算时给出步长h=0.4程序自动调节为0.0735892最小步长是第79步的0.000106687同样精度要求用固定步长计算大约要计算62000步计算方法引论(第三版)Stiff问题Stiff现象

方程通解Reλ<0,|Reλ|愈大解衰减越快方程组二解衰减快慢差千倍系统解衰减速度局部地决定于Jacobi矩阵的特征值大小方程组求稳态解(暂态解能忽略后)困难截断误差由慢变分量确定(方程组中y,对应小特征值).稳定性的要求决定于快变化分量(方程组中z,对应大特征值).必须依稳定性要求采用小步长,步数多,过程长.求稳态解,