matlab常微分方程初值问题的数值解法

Step 1: 先用显式欧拉公式作预测，算出 yn+1 = y n + h f ( x n, y n) 先用显式欧拉公式作预测 显式欧拉公式作预测， Step 2: 再将 yn+1 代入隐式梯形公式的右边作校正，得到 代入隐式梯形公式的右边作校正 隐式梯形公式的右边作校正，
yn+1
h = y n + [ f ( x n, y n) + f ( x n+1 , yn+1 )] 2
y ( x 2 ) ≈ y ( x 0 ) + 2 h f ( x1 , y ( x1 ))
y′( x1 ) ≈ y ( x 2 ) − y ( x0 ) 2h
y n +1 = y n −1 + 2h f (xn , y n ) n = 1, 2, ⋯ 改进欧拉法 /* modified Euler’s method */
改 进 E u le r 公 式
数值分析
数值分析
以上两组公式都使用函数f ( x , y )在某些点上的 值的线性组合来计算y( xn +1 )的近似值yn +1。 的值，为一阶方法。 Euler公式： 公式 Euler公式：每步计算一次f ( x , y )的值，为一阶方法。 改进Euler公式：需计算两次f ( x , y )的值，二阶方法。 的值，二阶方法。 改进Euler公式 Euler公式：
例 ：用尤 拉公 式和 改进 的尤拉 公式 解初 值问 题 2x ' y = y− y y(0) = 1. (0 < x < 1);
数值分析
数值分析
解 ： 取 步 长 h = 0 .1， 尤 拉 公 式 为 ： yn+1 2 xn ). = yn + h( yn − yn
2 xn y p = y n + h ( y n − y ); n 2 x n+1 ); 改 进 的 尤 拉 公 式 为 ： yc = yn + h( y p − yp 1 y n + 1 = ( y p + y c ). 2 计算结果略。
y ( x1 ) ≈ y0 + h f ( x1 , y ( x1 ))
yn+1 = yn + h f ( xn+1 , yn+1 ) (n = 0,1,2, ... )
同时出现在等式的两边，不能直接得到， 由于未知数 yi+1 同时出现在等式的两边，不能直接得到，故 称为隐式 欧拉公式，而前者称为显式 称为隐式 /* implicit */ 欧拉公式，而前者称为显式 /* explicit */ 欧拉公式。 欧拉公式。 一般先用显式计算一个初值，再迭代求解。 一般先用显式计算一个初值， 迭代求解。 求解
数值分析
常微分方程数值解
数值分析
数值分析
第一节 求解初值问题数值方法的基本原理
一、初值问题的数值解
考虑一阶常微分方程的初值问题 /* Initial-Value Problem */:
dy = f ( x, y) dx y(a ) = y0 x ∈ [a , b ]
(10(10-1)
∂f ( j − 2) ∂f ( j − 2) = + f ≡ f ( j −1) ; ∂x ∂y
数值分析
数值分析
Runge-Kutta法的基本思想（2） 法的基本思想（ ） 法的基本思想
如 果 将 E u le r 公 式 与 改 进 E u le r 公 式 写 成 下 列 形 式 ： E u le r 公 式 yn+1 = yn + h K 1 K 1 = f ( xn , yn ) 1 1 yn+1 = yn + h( 2 K 1 + 2 K 2 ) K = f (x , y ) 1 n n K 2 = f ( xn + h, yn + hK 1 )
只要 f (x, y) 在[a, b] × R1 上连续，且关于 y 满足 Lipschitz 条 上连续， 件，即存在与 x, y 无关的常数 L 使 | f ( x, y1 ) − f ( x, y2 ) | ≤ L| y1 − y2 | 数值解 都成立，则上述IVP存 对任意定义在 [a, b] 上的 y1(x) 和 y2(x) 都成立，则上述 存 在唯一解。 在唯一解。 要计算出解函数 y(x) 在一系列节点 a = x0< x1<…< xn= b 处的近似值 y i ≈ y ( x i ) ( i = 1, ... , n ) 为步长，通常采用等距节点 等距节点， 节点间距 hi = xi +1 − xi (i = 0, ... , n −1) 为步长，通常采用等距节点， 常数)。 即取 hi = h (常数 。 常数
数值分析
数值分析
忽略高阶项,取近似值可得到 忽略高阶项 取近似值可得到Euler公式 取近似值可得到 公式
yn+1 = yn + h f ( xn , yn ) (n = 0,1,2, ... )
3. 数值积分法区间
方程y' = f ( x, y )在 区间[ xn , xn+1 ]上 积分 将
∫
数值分析
数值分析
求解(10-1)最基本的方法是单步法 求解(10-1)最基本的方法是单步法 (10 最基本的方法是 单步法:从初值 开始,依次求出 单步法 从初值 y0 开始 依次求出 y1 , y2 , ⋯ ,后一步的值 yn +1 后一步的值 只依靠前一步的 yn 典型的单步法是Euler(欧拉 方法 其计算格式是 欧拉)方法 其计算格式是: 典型的单步法是 欧拉 方法,其计算格式是
法为例说明构造IVP问题数值方法的三种基本途径 以Euler法为例说明构造 法为例说明构造 问题数值方法的三种基本途径 1. 数值微分法 用差商代替微商 数值微分法,用差商代替微商 向前差商近似导数
y( x1 ) − y( x0 ) h y ( x1 ) ≈ y ( x 0 ) + hy ′( x 0 ) = y 0 + h f ( x 0 , y 0 ) 记为 y 1 y′( x0 ) ≈
yn+1 = yn + h f ( xn , yn ) (n = 0,1,2, ... )
例:求解常微分方程初值问题 求解常微分方程初值问题
y' = − y + x + 1 y ( 0) = 1
取步长 h = 0.1, 计算到 x = 0.5
x≥0
解 : f ( x , y ) = − y + x + 1,由Euler 公式
x n +1 xn
y dx = ∫
'
x n +1 xn
f ( x , y )dx ( n = 0,1, ⋯)
用 y n + 1 , y n 代 替 y ( x n + 1 ), y ( x n ), 对 右 端 积 分 采 用 取左端点的矩形公式
∫
则有
x n +1 xn
f ( x , y ) dx ≈ h f ( x n , y n )
数值分析
数值分析
梯形公式 /* trapezoid formula */
— 显、隐式两种算法的平均 隐式两种算法的平均
h y n +1 = y n + [f (xn , y n ) + f (xn+1 , y n +1 )] (n = 0,1, 2 ⋯) 2
中点欧拉公式 /* midpoint formula */ 中心差商近似导数
'
但由于公式中各阶偏导数计算复杂，不实用。
数值分析
数值分析
y ' = f ≡ f (0) ; ∂f (0) ∂f (0) y '' = + f ≡ f (1) ; ∂x ∂y ∂f (1) ∂f (1) y ''' = + f ≡ f (2) ; ∂y ∂x y( j ) 一般地有 j = 2, 3,⋯
数值分析
数值分析
第二节 高精度的单步法
在高精度的单步法中,应用最广泛的是Runge-Kutta(龙格 在高精度的单步法中,应用最广泛的是Runge-Kutta(龙格 Runge 库塔) -库塔)方法 法的基本思想（ ） 一、Runge-Kutta法的基本思想（1） 法的基本思想
若 用 p阶 Taylor 多 项 式 近 似 函 数 y ( x n + 1 )有 ： h2 " h p ( p) y n + 1 ≈ y ( x n + 1 ) = y ( x n ) + hy ( x n ) + y ( xn ) + ⋯ + y ( xn ) 2! P! 其 中 y ' ( x ) = f ( x , y ), y '' ( x ) = f x' ( x , y ) + f y' ( x , y ) f ( x , y ), ⋯ ⋯ 。
yn+1 = yn + h f ( xn , yn ) (n = 0,1,2, .、Euler法的改进及梯形公式 法的改进及梯形公式
隐式欧拉法 /* implicit Euler method */ 向后差商近似导数
y′( x1 ) ≈ y( x1 ) − y( x0 ) h
数值分析
数值分析
二、二阶龙格－库塔方法 二阶龙格－
一 般 地 ， RK方 法 设 近 似 公 式 为
p yn+1 = yn + h∑ ci K i i =1 K 1 = f ( xn , yn ) i −1 K i = f ( x n + a i h , y n + h ∑ bij K j ) ( i = 2, 3 ⋯ , p ) j=1 其 中 a i， bij , c i 都 是 参 数 ， 确 定 它 们 的 原 则 是 使 近 似 公 式
数值分析
数值分析
y n +1 = y n +
h f (xn , y n ) + f ( xn +1 , y n + h f (xn , y n ) ) 2
(n = 0, 1, 2⋯)
注：此法亦称为预测-校正法 /* predictor-corrector 此法亦称为预测predictor预测 */。一方面它有较高精度， method */。一方面它有较高精度，同时可以看 到它是个单步递推格式， 单步递推格式 到它是个单步递推格式，比隐式公式的迭代求 解过程简单 后面将看到，它的稳定性高于 简单。 稳定性高于显 解过程简单。后面将看到，它的稳定性高于显 式欧拉法。 式欧拉法。
线性组合来构造近似公式，构造是要求近似公式在 ( x n , y n )处 的 Taylor 展 开 式 与 解 y ( x )在 x n 处 的 Taylor 展 开 式的前面几项重合，从而使近似公式达到所需要的阶 RK方 数 。 即 避 免 求 偏 导 ， 又 提 高 了 方 法 的 精 度 ， 此 为 RK 方 法的基本思想。
数值分析
数值分析
yn+1 = yn + h( − yn + xn + 1)
代入h = −0.1, 有yn+1 = 0.9 yn + 0.1( xn + 1), 依次算得果如下:
n=0 xn = 0 1 0.1 2 0.2 3 0.3 4 0.4 5 0.5
yn = 1.0 1.0 1.01 1.029 1.0561 1.09049
数值分析
数值分析
Runge-Kutta法的基本思想（3） 法的基本思想（ ） 法的基本思想
p yn+1 = yn + h∑ ci K i i=1 K 1 = f ( xn , yn ) i−1 K i = f ( x n + a i h , y n + h ∑ b ij K j ) ( i = 2 , 3 ⋯ , p ) j=1 于 是 可 考 虑 用 函 数 f ( x , y )在 若 干 点 上 的 函 数 值 的
直接解微分方程可得精确解 : y = f ( x) = x + e − x , 故x5 = 0.5, y(0.5) = 1.106531
由此可见,Euler公式的近似值接近方程的精确值 公式的近似值接近方程的精确值. 由此可见 公式的近似值接近方程的精确值
数值分析
数值分析
二、构造初值问题数值方法的基本途径
yn+1 = yn + h f ( xn , yn ) (n = 0,1,2, ... )
2. Taylor展开法 展开法
将 y ( x n + h ) 在点 x n 作 Taylor 展开 h 2 '' y ( x n + h ) = y ( x n ) + hy ' ( x n ) + y (ξ ) 2!