数值分析第九章节 用Python实现常微分方程初值问题的数值解法 |
您所在的位置：网站首页 › 迭代法解线性方程组的充分必要条件 › 数值分析第九章节 用Python实现常微分方程初值问题的数值解法 |
数值分析第九章节 用Python实现常微分方程初值问题的数值解法
|
参考书籍：数值分析 第五版 李庆杨 王能超 易大义编 第9章 常微分方程初值问题的数值解法 文章声明：如有发现错误，欢迎批评指正 文章目录 欧拉法后退的欧拉方法梯形方法改进欧拉公式补充龙格—库塔方法线性多步法阿当姆斯显示与隐式公式 9.1引言 科学技术中的很多问题可以使用常微分方程的定解描述，包括初值问题以及边值问题两个大类。考虑一阶常微分方程的初值问题， y ′ = f ( x , y ) , x ∈ [ x 0 , b ] , y ( x 0 ) = y 0 . y'=f(x,y),x \in[x_0,b],y(x_0)=y_0. y′=f(x,y),x∈[x0,b],y(x0)=y0.如果存在实数 L > 0 L>0 L>0,使得 ∣ f ( x , y 1 ) − f ( x , y 2 ) ∣ ≤ L ∣ y 1 − y 2 ∣ , ∀ y 1 , y 2 ∈ R , |f(x,y_1)-f(x,y2)|\leq L|y_1-y_2|,\forall y_1,y_2\in R, ∣f(x,y1)−f(x,y2)∣≤L∣y1−y2∣,∀y1,y2∈R,则称f关于y满足利普西茨条件， L L L为 f f f的利普西茨常数。 定理1(唯一性)： 设 f f f在区域 D = { ( x , y ) ∣ a ≤ x ≤ b , y ∈ R } D=\{(x,y)|a\leq x\leq b,y \in R\} D={(x,y)∣a≤x≤b,y∈R}上连续，关于y则满足利普西茨条件，则对任意 x 0 ∈ [ a , b ] , y 0 ∈ R x_0\in[a,b],y_0\in R x0∈[a,b],y0∈R，常微分的初值问题在 x ∈ [ a , b ] x\in[a,b] x∈[a,b]上存在唯一的连续可微解 y ( x ) y(x) y(x)。 定理2(敏感性的分析)： 设 f f f在区域 D D D上连续，且关于y满足利普西茨条件，设初值问题 y ′ ( x ) = f ( x , y ) , y ( x 0 ) = s y'(x)=f(x,y),y(x_0)=s y′(x)=f(x,y),y(x0)=s的解为 y ( x , s ) y(x,s) y(x,s)，则 ∣ y ( x , s 1 1 ) − y ( x , s 2 ) ∣ ≤ e L ∣ x − x 0 ∣ ∣ s 1 − s 2 ∣ |y(x,s1_1)-y(x,s_2)|\leq e^{L|x-x_0|}|s_1-s_2| ∣y(x,s11)−y(x,s2)∣≤eL∣x−x0∣∣s1−s2∣。显然，当 L L L较大时认为坏条件，当 L L L较小时认为好条件。实际问题之中遇到的常微分方程大部分没有解析解，所以主要依靠数值解法。 9.2简单的数值解法 9.2.1欧拉法与后退欧拉法 欧拉法欧拉方法：使用公式 y n + 1 = y n + h f ( x n , y n ) y_{n+1}=y_n+hf(x_n,y_n) yn+1=yn+hf(xn,yn)逐次地推导算出最终结果。 例1：求解初值问题 { y ′ = y − 2 x y , 0 < x < 1 , y ( 0 ) = 1. \left\{\begin{matrix}y'=y-\frac{2x}{y},0y′=−100y,y(0)=1. 可以解得： y ( x ) = e − 100 x y(x)=e^{-100x} y(x)=e−100x lt1=[] x_initial=0;y_initial=1 for i in range(4): y_initial=y_initial+1/40*(-100*y_initial);x_initial=i/40 lt1.append(y_initial) lt2=[] x_initial=0;y_initial=1 for i in range(4): y_initial=y_initial/(1+100*0.025) lt2.append(y_initial) for i in range(4): print("欧拉法：{:>8.4f} 后退的欧拉法：{:>8.4f}".format(lt1[i],lt2[i]))
使用 y n + 1 = y n + h 2 [ f ( x n , y n ) + f ( x n + 1 , y n + 1 ) ] y_{n+1}=y_n+\frac{h}{2}[f(x_n,y_n)+f(x_{n+1},y_{n+1})] yn+1=yn+2h[f(xn,yn)+f(xn+1,yn+1)]为公式作为单步推导的方法叫梯形方法。求解方法等同后退的欧拉法。 使用例4作为演示(就是上面的那一个) lt=[] x_initial=0;y_initial=1 for i in range(4): y_initial=(1-50*0.025)*y_initial/(1+50*0.025) lt.append(y_initial) for i in range(4): print("梯形方法：{:>8.4f}".format(lt[i]))
同样我们这里讨论梯形方法的稳定性： y n + 1 = 1 − h 2 λ 1 + h 2 λ y n + 1 y_{n+1}=\frac{1-\frac{h}{2}\lambda}{1+\frac{h}{2}\lambda}y_{n+1} yn+1=1+2hλ1−2hλyn+1则 ∣ 1 − h 2 λ 1 + h 2 λ ∣ ≤ 1 |\frac{1-\frac{h}{2}\lambda}{1+\frac{h}{2}\lambda}|\leq1 ∣1+2hλ1−2hλ∣≤1，同时这里，我们给出一个定义： 如果一个数值方法的稳定域包含了 { h λ ∣ R e ( h λ ) < 0 ∣ \{h\lambda|Re(h\lambda) |
可以看到
h
h
h为0.025时，欧拉法不稳定但是后退的欧拉法稳定。
h
h
h为0.005时，欧拉法变稳定(这个部分我没有做)。所以得出结论稳定性不仅与方法还与步长有关。现在我们来具体地探讨欧拉方法的稳定性。模型方程
y
′
=
λ
y
y'=\lambda y
y′=λy的欧拉公式为
y
n
+
1
=
(
1
+
h
λ
)
y
n
y_{n+1}=(1+h\lambda)y_n
yn+1=(1+hλ)yn。设在节点
y
n
y_n
yn有一扰动值
ϵ
n
\epsilon_n
ϵn，在
y
n
+
1
y_{n+1}
yn+1产生大小为
ϵ
n
+
1
\epsilon_{n+1}
ϵn+1扰动值。所以有
ϵ
n
+
1
=
(
1
+
h
λ
)
ϵ
n
\epsilon_{n+1}=(1+h\lambda)\epsilon_n
ϵn+1=(1+hλ)ϵn。假设差分方程解是不增长的，根据稳定性的定义有
∣
1
+
h
λ
∣
≤
1
|1+h\lambda|\leq1
∣1+hλ∣≤1。对于后退的欧拉法有
y
n
+
1
=
1
1
−
h
λ
y
n
y_{n+1}=\frac{1}{1-h\lambda}y_n
yn+1=1−hλ1yn则
∣
1
1
−
h
λ
∣
<
1
|\frac{1}{1-h\lambda}|1
∣1−hλ∣>1。
【本文地址】
今日新闻 |
推荐新闻 |