《数值分析与算法》第八讲-常微分方程初值问题

数值分析(8)NumericalAnalysisWenjianYu2第八章常微分方程初值问题

(主要是前3节)常微分方程初值问题WenjianYu3常微分方程基本概念常微分方程

WenjianYu4

常微分方程

WenjianYu5

初值问题:

常微分方程–例子1含电容元件的电路问题通过节点分析法得到微分方程组WenjianYu6问题的解反映了电容充/放电过程C1R2R1R3C2电流/电压关系节点电流方程t常微分方程–例子2双联摆的运动两个摆锤(重物),刚性杆的重量可忽略不考虑摩擦力,运动不会停止,且在初始角

度较大时摆锤的轨迹呈现混沌现象WenjianYu7求解微分方程初值问题,得到

摆锤的运动规律Matlab演示swinger

常微分方程

WenjianYu8

线性齐次常系数微分方程

实际的问题基本上都是稳定的!

(由于历史原因)常微分方程

WenjianYu9

局部稳定简单方法与有关概念WenjianYu10简单的初值问题数值解法

初值问题的数值解法

WenjianYu11

否则为多步法否则为隐格式方法欧拉法

WenjianYu12“左矩形”求积公式

h=0.1h=0.050.11.0000001.0048370.051.0000000.31.0350920.21.0100001.0187310.11.0025000.351.0483370.31.0290001.0408180.151.0073750.41.0634200.41.0561001.0703200.21.0145060.451.0802490.51.0904901.1065310.251.0237810.51.098737

步长h=0.1,和0.05步长小的更准

数值解法的稳定性

WenjianYu13

-10

欧拉法解模型问题的稳定区域

数值解法的稳定性

WenjianYu14-10

欧拉法稳定

数值解法的稳定性

WenjianYu15

00.0250.050.0750.10.1250.151-1.52.25-3.3755.0625-7.5937511.390610.0820850.0067380.0005534.5410-53.7310-63.0610-7这里设的h太大!计算结果如下表:数值解法的局部截断误差

WenjianYu16整体误差

稳定的问题，整体误差小于局部误差之和不稳定的问题呢？一般仅能控制局部误差整体误差？~局部误差数值解法的局部截断误差

WenjianYu17

欧拉法是一阶方法我们讨论的所有方法都至少有1阶准确度数值解法的收敛性：随着h0，误差0向后欧拉法与梯形法从数值积分的角度推导向后欧拉法:梯形法:两者均为单步、隐格式方法,每步计算要求解(非线性)方程例8.6:用向后欧拉法求解WenjianYu18右矩形

梯形

00.0250.050.0750.10.1250.1510.0066630.0019040.00054410.0820850.0067380.0005534.5410-53.7310-63.0610-7向后欧拉法

WenjianYu19

准确解

01稳定区域

无条件稳定(unconditionallystable)!向后欧拉法

WenjianYu20

具有1阶准确度!向后欧拉法与梯形法

WenjianYu21

稳定的条件是:

思考无条件稳定!

具有2阶准确度简单方法与有关概念WenjianYu22Runge-Kutta方法在欧拉法基础上改进再增加一次函数求值:数值积分的中矩形或梯形公式利用欧拉法算半个步长的结果,估算中点处被积函数值先用欧拉法估计区间终点处斜率,再用它与

起始点斜率的平均值算一整步Runge-Kutta方法

中矩形梯形公式(中点公式)(Heun方法/改进的欧拉法)都比欧拉法准确均为2级R-K公式WenjianYu23

Runge-Kutta方法只能估算

,令

WenjianYu24

Runge-Kutta方法其他,用所有前面点的信息

WenjianYu25积分节点………几种显式R-K公式参数的值不按具体数值

积分公式设置,而根据

准确度阶数要求设置2级公式：改进欧拉法、

中点公式经典4级、4阶Runge-Kutta法

(1905)Runge-Kutta方法WenjianYu26

Runge-Kutta方法此时局部截断误差只要对非模型问题也有相同结论！WenjianYu27如

例8.7:用2阶改进欧拉、3阶Ralston、

4阶经典R-K解问题,h=0.1,算到y(2)精确解为Runge-Kutta方法r4对应的r级R-K公式有r阶准确度高于4阶的公式很少单独使用WenjianYu28r>4对应的r级R-K公式达不到r阶准确度二阶Heun三阶Ralston四阶Runge-Kutta准确值10.40.40.40.41.10.4756410.4746260.47463830.47463821.20.5834080.5813640.58138680.58138671.30.7281350.7250340.72506630.72506621.40.9153290.9111370.91117730.91117711.51.1511101.1457851.14583361.14583331.61.4421691.4356641.43572031.43572001.71.7957381.7880041.78806741.78806711.82.2195782.2105612.21063152.21063111.92.7219612.7116062.71168362.71168322.03.3116653.2999163.30000043.3000000

Runge-Kutta方法(此时)一般显式单步法

恰好与欧拉法一样显格式,都不是无条件稳定的局部截断误差判断单步法收敛性的简便方法WenjianYu29

简单方法与有关概念WenjianYu30多步法

多步法WenjianYu31

(线性m步法)

固定步长hTaylor展开法求线性m步法的系数

这些系数应该等于0例8.8:求两步法公式中参数值多步法WenjianYu32

满足它们才可能收敛(相容性)

,

有二阶准确度

多步法WenjianYu33

同例8.8的结果!

多步法公式中包含p个待定参数，至少可达到p-1阶准确度

多步法WenjianYu34

Vandermonde阵T,非奇异

插值节点函数值

常用的多步法公式Adams公式的推导：用插值多项式近似被积函数例8.10:推导m=4对应的显式Adams公式多步法WenjianYu35

可证明满足最高准确度阶数

类似地算其他系数，得

显式四阶Adams-

Bashforth公式(单项式函数代入法)Adams公式几种显式公式几种隐式公式多步法WenjianYu36

阶数稳定阈值误差常数11

-21/223/2-1/2

-15/12323/12-16/125/12

-6/113/8455/24-59/2437/24-9/24-3/10251/720欧拉法

阶数稳定阈值误差常数11

--1/221/21/2

--1/1235/128/12-1/12

-6-1/2449/2419/24-5/241/24-3-19/720向后欧拉法梯形法并非无条件稳定!

多步法WenjianYu37

用Matlab求解初值问题WenjianYu38用Matlab解ODE-IVP

Matlab中的ODE-IVP求解器

WenjianYu39[T,

Y,

TE,

YE,

IE]=solver(odefun,

tspan,

y0,

options)

求解单个ODE火焰燃烧问题当点燃一根火柴时,火焰迅速增大直到一个临界体积,然后维持这一体积不变,此时火焰内部燃烧耗费的氧气和其表面现存的氧气达到了一种平衡.火焰(近似为球)半径y满足ODE>>f=@(t,y)y^2-y^3;>>ode23(f,[0,2.0e4],1e-4)WenjianYu40设初始半径=0.0001例8.15尝试用ode23s求解它求解ODE方程组

WenjianYu41例8.16functionydot=myode2(t,y);ydot=[y(2);6*t];%列向量

双联摆问题的求解WenjianYu42隐格式非线性常微分方程初值

问题,初值为[,,0,0]T

odeset可以设置质量矩阵质量矩阵swinger_solve.m很大时混沌现象:确定的但不可预测(非线性科学)

小结WenjianYu43函数名内部算法说明ode23显格式，单步法，采用BS23算法的自动变步长R-K方法对于精度要求不高的情况，效率好于ode45ode45显格式，单步法，包含一个4阶和5阶公式的自动变步长R-K方法一般情况下，首先尝试使用它来求解ode113显格式，多步法，采用变阶数的Adam