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

常微分方程初值问题的数值解法制作人：赵文波常微分方程初值问题的数值解法制作人：赵文波11引言一般的一阶常微分方程初值问题

y`=f(t,y),a≤t≤b y(a)=η(1.1)定理一 如果f(t,y)在带形区域内R={(t,y

)|a≤t≤b,-∞≤u≤+∞}中连续，且关于y满足Lipschitz条件：存在常数L，使得|f(t,y1)–f(t,y2)|≤L|y1-y2| (1.2)定理二

如果f(t,y)在R={(t,y)|a≤t≤b,-∞≤u≤+∞}中连续，且关于y满足Lipschitz条件，那么初值问题是适定的

(1.3)1引言一般的一阶常微分方程初值问题 22离散变量和离散误差离散化过程：把初值为(1.1)的精确解y(t)在一系列的离散点：t1,t2,……

tN处的近似值y1y2……yN用y(tn)表示在t=tn

处的近似值，n=1,2,……N.离散化方法 1.差商代替导数的方法 2.Taylor级数法误差 1.局部离散误差

2.局部截断误差

3.整体离散误差2离散变量和离散误差离散化过程：把初值为(1.1)的精确解33单步法单步法的一般形式 显式方法：

yn+1=yn+hФ(tn,yn,h),n=0,1,2,…N-1, y0=η,(3.1)

或 隐式方法：

yn+1=yn+hФ(tn,yn,yn+1,h),n=0,1,2,…N-1, y0=η,(3.2)

Ф为增量函数，N是正整数，h=(b-a)/N3单步法单步法的一般形式4Euler法算法10.1

y`=f(t,y),a≤t≤b；y(a)=ηInput端点a,b;区间等分数N，初值ηOuputy(t)在t的N个点处的近似值Step1h←(b-a)/N;t←a;y←η.Step2For

i=1,2,…,NdoStep3-4

Step3y←y+hf(t,y); t←a+ih. Step4Output(t,y).Step5return.Euler法算法10.15改进的Euler法解初值问题的梯形计算公式

yn+1=yn+(h/2)[f(tn,yn)+f(tn+1,yn+1)],n=0,1,2,…N-1 (3.8) y0=η; h=(b-a)/N

局部离散误差为Rn=-(h3/12)y```(ξn) y(0)n+1=yn+hf(tn,yn), yn+1=(1/2)[y(0)n+1+yn+hf(tn+1,y(0)n+1

)] n=0,1,2,…N-1(3.10)改进的Euler法解初值问题的梯形计算公式6二阶Runge-Kutta方法二阶Runge-Kutta方法7m阶Runge-Kutta显式方法m阶Runge-Kutta显式方法8Richardson外推法Richardson外推法94单步法的相容性和稳定性相容性收敛性稳定性4单步法的相容性和稳定性相容性104.1相容性单步法的一般形式

yn+1=yn+hФ(tn,yn,h),n=0,1,2,…N-1, y0=η,(4.1)定义1 Ф(t,y,0)=f(t,y)(4.3) 成立，则称单步法与微分初值问题(1.1)相容，(4.3)为相容条件。定理一假设Ф(t,y,0)关于h是连续的，若单步法与微分初值问题(1.1)相容，则它至少是一阶方法。4.1相容性单步法的一般形式114.2收敛性定义2假设微分方程(1.1)的右端函数f(t,y)在带形区域内R={(t,y)|a≤t≤b,-∞≤u≤+∞}中连续，且关于y满足Lipschitz条件。若对所有的t∈[a,b],limh→0tn=t固定yn=y(t)则称单步法是收敛的。

定理2若Ф(t,y,h)对于a≤t≤b,0≤h≤h0以及一切实数y，关于t,y,h满足Lipschitz条件，则单步法(4.1)收敛的充分必要条件是相容条件成立，即Ф(t,y,0)=f(t,y).4.2收敛性定义2假设微分方程(1.1)的右端函数12定理3在定理2的假设下，单步法(4.1)的局部离散误差满足(4.13),则其整体离散误差误差εn=y(tn)-yn满足估计式 |εn|≤eL(b-a)|ε0|+hpM(eL(b-a)–1)/L L是Ф(t,y,h)关于y满足Lipschitz条件的Lipschitz常数。定理3在定理2的假设下，单步法(4.1)的局部离散误差满134.3稳定性定义3如果存在常数h0及C,使得对任意的初始值y0,y0`,单步法(4.1)的相应的精确解yn,yn`,对所有的0<h≤h0,恒有 |yn-yn`

|≤C|y0-y0`|,nh≤b-a,则说单步法(4.1)是稳定。定理4若Ф(t,y,h)对于a≤t≤b,0≤h≤h0以及一切实数y，关于t,y,h满足Lipschitz条件，则单步法(4.1)是稳定的。4.3稳定性定义3如果存在常数h0及C,使得对任意的初始值14定义4对给定的微分方程和给定的步长h，若有单步法（显式或隐式）计算yn时有大小为δ的误差，即计算得yn`=yn+δ,而引起其后值ym(m>n)的变化小于δ(|ym-ym`

|<|δ|),则说单步法是绝对稳定的。 一般限于y`=μy(4.20)

考虑数值方法的绝对稳定性，μ为复常数，若对于所有μh∈(α,β),单步法都绝对稳定，称(α,β)为绝对稳定区间。定义4对给定的微分方程和给定的步长h，若有单步法（显式或155多步法5多步法165.1线性多步法

y`=f(t,y),a≤t≤b y(a)=η(1.1)线性k步法的一般公式

5.1线性多步法 y`=f(t,y),a≤t≤175.2Adams方法显式Adams方法隐式Adams方法5.2Adams方法显式Adams方法185.3预测-校正方法y(0)n+1=yn+hf(tn,yn)(5.20)

yn+1=(1/2)[y(0)n+1+yn+hf(tn+1,y(0)n+1

)](5.21)

n=0,1,2,…N-1;y0=η(5.20)起预测yn+1

的作用，(5.21)起校正作用。 记f(i)n=f(tn,y(i)n).用P表预测过程，C表校正过程，E表计算f的过程 P:y(0)n+1=yn+hfn,

E:f(0)n+1=f(tn+1,y(0)n+1) C:yn+1=yn+(h/2)(fn+f(0)n+1)

E:fn+1=f(tn+1,yn+1)5.3预测-校正方法y(0)n+1=yn+hf(t19 重复迭代P,E,Ct次可提高精度。 通常，把Adams隐式和

显式方法联合使用，构成预测-校正方法。 预测公式： y(0)n+1=yn+h[βk0fn+βk1fn-1+…+βkkfn-k) 校正公式：

y(i+1)n+1=yn+h[β*k0fn+β*k1fn-1+…+β*kkfn-k) 重复迭代P,E,Ct次可提高精度。205.4Hamming方法Milne方法建立线性多步法的待定系数法Hamming方法5.4Hamming方法Milne方法217线性多步法的相容性、收敛性和稳定性定义1若求解初值问题(1.1)的线性k步法(5.1)至少是一阶方法，则称他们是相容的。 记ρ(λ)=

αkλk+αk-1λk-1+…+α1λ

+α0，

σ(λ)=βkλk+βk-1λk-1+…+β1λ

+β0。 他们由线性k步法(5.1)完全确定。反之，若给定了ρ(λ)和σ(λ)，则他们唯一确定一个线性k步法。我们称ρ(λ)为线性k步法(5.1)的特征多项式。定理1线性k步法(5.1)相容的充分必要条件是

ρ(1)=0，ρ`(1)=σ(1)，7线性多步法的相容性、收敛性和稳定性定义1若求解初值问题(227.2收敛性定义2假设f(t,y)在R={(t,y)|a≤t≤b,-∞≤u≤+∞}中连续，且关于y满足Lipschitz条件。若对任意的t∈[a,b],但t→0,而a+nh=tn=t固定时，(5.1)的解yn收敛于问题(1.1)的解y(t)，则说线性k步法(5.1)是收敛的。定理2若线性k步法(5.1)是收敛，则必相容。7.2收敛性定义2假设f(t,y)在R={(t,y)|a237.3稳定性定义3若f(t,y)在R={(t,y)|a≤t≤b,-∞≤u≤+∞}中连续，且关于y满足Lipschitz条件。若存在正常数C和h0，使得当0<h<h0

时线性k步法(5.1)的任意两个yn和yn`满足不等式 maxnh

≤(b-a)|yn-yn`|≤CM0

其中M0=max0≤i≤k-1|yi-yi`|那么说线性k步法(5.1)是稳定的。引理若k为非负整数，数列

{εn}满足递推不等式 |εn|≤β+αh(ε0+ε1+…εn-1),n=k,k+1,…nh≤(b-a),

其中0≤

α，β,M0=max0≤i≤k-1|εi|, 则|εn|≤eα(b-a)(β+αhkM0),n=k,k+1,…nh≤(b-a)7.3稳定性定义3若f(t,y)在R={(t,y)|a24定理3线性k步法(5.1)稳定的充分必要条件是ρ(λ)满足特征根条件：ρ(λ)的所有根都在单位圆上，且在在单位圆上的根只能是单重根。定理4若线性k步法(5.1)收敛，则必稳定。定理5若线性k步法(5.1)相容且稳定，则必收敛稳定。定理3线性k步法(5.1)稳定的充分必要条件是ρ(λ)满257.4绝对稳定性定义4对给定的μ,h,若特征方程(7.24)的所有根的模都<1，则称线性k步法(7.23)绝对稳定。若对所有μh∈(α,β),(7.23)都绝对稳定，则说(α,β)为绝对稳定区间。εε

ε

αβρσ7.4绝对稳定性定义4对给定的μ,h,若特征方程(7.26常微分方程初值问题的数值解法制作人：赵文波常微分方程初值问题的数值解法制作人：赵文波271引言一般的一阶常微分方程初值问题

y`=f(t,y),a≤t≤b y(a)=η(1.1)定理一 如果f(t,y)在带形区域内R={(t,y

)|a≤t≤b,-∞≤u≤+∞}中连续，且关于y满足Lipschitz条件：存在常数L，使得|f(t,y1)–f(t,y2)|≤L|y1-y2| (1.2)定理二

如果f(t,y)在R={(t,y)|a≤t≤b,-∞≤u≤+∞}中连续，且关于y满足Lipschitz条件，那么初值问题是适定的

(1.3)1引言一般的一阶常微分方程初值问题 282离散变量和离散误差离散化过程：把初值为(1.1)的精确解y(t)在一系列的离散点：t1,t2,……

tN处的近似值y1y2……yN用y(tn)表示在t=tn

处的近似值，n=1,2,……N.离散化方法 1.差商代替导数的方法 2.Taylor级数法误差 1.局部离散误差

2.局部截断误差

3.整体离散误差2离散变量和离散误差离散化过程：把初值为(1.1)的精确解293单步法单步法的一般形式 显式方法：

yn+1=yn+hФ(tn,yn,h),n=0,1,2,…N-1, y0=η,(3.1)

或 隐式方法：

yn+1=yn+hФ(tn,yn,yn+1,h),n=0,1,2,…N-1, y0=η,(3.2)

Ф为增量函数，N是正整数，h=(b-a)/N3单步法单步法的一般形式30Euler法算法10.1

y`=f(t,y),a≤t≤b；y(a)=ηInput端点a,b;区间等分数N，初值ηOuputy(t)在t的N个点处的近似值Step1h←(b-a)/N;t←a;y←η.Step2For

i=1,2,…,NdoStep3-4

Step3y←y+hf(t,y); t←a+ih. Step4Output(t,y).Step5return.Euler法算法10.131改进的Euler法解初值问题的梯形计算公式

yn+1=yn+(h/2)[f(tn,yn)+f(tn+1,yn+1)],n=0,1,2,…N-1 (3.8) y0=η; h=(b-a)/N

局部离散误差为Rn=-(h3/12)y```(ξn) y(0)n+1=yn+hf(tn,yn), yn+1=(1/2)[y(0)n+1+yn+hf(tn+1,y(0)n+1

)] n=0,1,2,…N-1(3.10)改进的Euler法解初值问题的梯形计算公式32二阶Runge-Kutta方法二阶Runge-Kutta方法33m阶Runge-Kutta显式方法m阶Runge-Kutta显式方法34Richardson外推法Richardson外推法354单步法的相容性和稳定性相容性收敛性稳定性4单步法的相容性和稳定性相容性364.1相容性单步法的一般形式

yn+1=yn+hФ(tn,yn,h),n=0,1,2,…N-1, y0=η,(4.1)定义1 Ф(t,y,0)=f(t,y)(4.3) 成立，则称单步法与微分初值问题(1.1)相容，(4.3)为相容条件。定理一假设Ф(t,y,0)关于h是连续的，若单步法与微分初值问题(1.1)相容，则它至少是一阶方法。4.1相容性单步法的一般形式374.2收敛性定义2假设微分方程(1.1)的右端函数f(t,y)在带形区域内R={(t,y)|a≤t≤b,-∞≤u≤+∞}中连续，且关于y满足Lipschitz条件。若对所有的t∈[a,b],limh→0tn=t固定yn=y(t)则称单步法是收敛的。

定理2若Ф(t,y,h)对于a≤t≤b,0≤h≤h0以及一切实数y，关于t,y,h满足Lipschitz条件，则单步法(4.1)收敛的充分必要条件是相容条件成立，即Ф(t,y,0)=f(t,y).4.2收敛性定义2假设微分方程(1.1)的右端函数38定理3在定理2的假设下，单步法(4.1)的局部离散误差满足(4.13),则其整体离散误差误差εn=y(tn)-yn满足估计式 |εn|≤eL(b-a)|ε0|+hpM(eL(b-a)–1)/L L是Ф(t,y,h)关于y满足Lipschitz条件的Lipschitz常数。定理3在定理2的假设下，单步法(4.1)的局部离散误差满394.3稳定性定义3如果存在常数h0及C,使得对任意的初始值y0,y0`,单步法(4.1)的相应的精确解yn,yn`,对所有的0<h≤h0,恒有 |yn-yn`

|≤C|y0-y0`|,nh≤b-a,则说单步法(4.1)是稳定。定理4若Ф(t,y,h)对于a≤t≤b,0≤h≤h0以及一切实数y，关于t,y,h满足Lipschitz条件，则单步法(4.1)是稳定的。4.3稳定性定义3如果存在常数h0及C,使得对任意的初始值40定义4对给定的微分方程和给定的步长h，若有单步法（显式或隐式）计算yn时有大小为δ的误差，即计算得yn`=yn+δ,而引起其后值ym(m>n)的变化小于δ(|ym-ym`

|<|δ|),则说单步法是绝对稳定的。 一般限于y`=μy(4.20)

考虑数值方法的绝对稳定性，μ为复常数，若对于所有μh∈(α,β),单步法都绝对稳定，称(α,β)为绝对稳定区间。定义4对给定的微分方程和给定的步长h，若有单步法（显式或415多步法5多步法425.1线性多步法

y`=f(t,y),a≤t≤b y(a)=η(1.1)线性k步法的一般公式

5.1线性多步法 y`=f(t,y),a≤t≤435.2Adams方法显式Adams方法隐式Adams方法5.2Adams方法显式Adams方法445.3预测-校正方法y(0)n+1=yn+hf(tn,yn)(5.20)

yn+1=(1/2)[y(0)n+1+yn+hf(tn+1,y(0)n+1

)](5.21)

n=0,1,2,…N-1;y0=η(5.20)起预测yn+1

的作用，(5.21)起校正作用。 记f(i)n=f(tn,y(i)n).用P表预测过程，C表校正过程，E表计算f的过程 P:y(0)n+1=yn+hfn,

E:f(0)n+1=f(tn+1,y(0)n+1) C:yn+1=yn+(h/2)(fn+f(0)n+1)

E:fn+1=f(tn+1,yn+1)5.3预测-校正方法y(0)n+1=yn+hf(t45 重复迭代P,E,Ct次可提高精度。 通常，把Adams隐式和

显式方法联合使用，构成预测-校正方法。 预测公式： y(0)n+1=yn+h[βk0fn+βk1fn-1+…+βkkfn-k) 校正公式：

y(i+1)n+1=yn+h[β*k0fn+β*k1fn-1+…+β*kkfn-k) 重复迭代P,E,Ct次可提高精度。465.4Hamming方法Milne方法建立线性多步法的待定系数法Hamming方法5.4Hamming方法Milne方法477线性多步法的相容性、收敛性和稳定性定义1若求解初值问题(1.1)的线性k步法(5.1)至少是一阶方法，则称他们是相容的。 记ρ(λ)=

αkλk+αk-1λk-1+…+α1λ

+α0，

σ(λ)=βkλk+βk-1λk-1+…+β1λ

+β0。 他们由线性k步法(5.1)完全确定。反之，若给定了ρ(λ)和σ(λ)，则他们唯一确定一个线性k步法。我们称ρ(λ)为线性k步法(5.1)的特征多项式。定理1线性k步法(5.1)相容的充分必要条件是

ρ(1)=0，ρ`(1)=σ(1)，7线性多步法的相容性、收敛性和稳定性定义1若求解初值问题(487.