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

常微分方程初值问题的数值解法制作人：赵文波1引言一般的一阶常微分方程初值问题y`=f(t,y),a≤t≤by(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)2离散变量和离散误差离散化过程：把初值为(1.1)的精确解y(t)在一系列的离散点：t1,t2,……tN处的近似值y1y2……yN用y(tn)表示在t=tn处的近似值，n=1,2,……N.离散化方法1.差商代替导数的方法2.Taylor级数法误差1.局部离散误差2.局部截断误差3.整体离散误差3单步法单步法的一般形式显式方法：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)/NEuler法算法10.1y`=f(t,y),a≤t≤b；y(a)=ηInput端点a,b;区间等分数N，初值ηOuputy(t)在t的N个点处的近似值Step1h←(b-a)/N;t←a;y←η.Step2Fori=1,2,…,NdoStep3-4Step3y←y+hf(t,y);t←a+ih.Step4Output(t,y).Step5return.改进的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)二阶Runge-Kutta方法m阶Runge-Kutta显式方法Richardson外推法4单步法的相容性和稳定性相容性收敛性稳定性4.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.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).定理3在定理2的假设下，单步法(4.1)的局部离散误差满足(4.13),则其整体离散误差误差εn=y(tn)-yn满足估计式|εn|≤eL(b-a)|ε0|+hpM(eL(b-a)–1)/LL是Ф(t,y,h)关于y满足Lipschitz条件的Lipschitz常数。4.3稳定性定义3如果存在常数h0及C,使得对任意的初始值y0,y0`,单步法(4.1)的相应的精确解yn,yn`,对所有的0h≤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对给定的微分方程和给定的步长h，若有单步法（显式或隐式）计算yn时有大小为δ的误差，即计算得yn`=yn+δ,而引起其后值ym(mn)的变化小于δ(|ym-ym`||δ|),则说单步法是绝对稳定的。一般限于y`=μy(4.20)考虑数值方法的绝对稳定性，μ为复常数，若对于所有μh∈(α,β),单步法都绝对稳定，称(α,β)为绝对稳定区间。5多步法5.1线性多步法y`=f(t,y),a≤t≤by(a)=η(1.1)线性k步法的一般公式5.2Adams方法显式Adams方法隐式Adams方法5.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)重复迭代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)5.4Hamming方法Milne方法建立线性多步法的待定系数法Hamming方法7线性多步法的相容性、收敛性和稳定性定义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.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.3稳定性定义3若f(t,y)在R={(t,y)|a≤t≤b,-∞≤u≤+∞}中连续，且关于y满足Lipschitz条件。若存在正常数C和h0，使得当0hh0时线性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)定理3线性k步法(5.1)稳定的充分必要条件是ρ(λ)满足特征根条件：ρ(λ)的所有根都在单位圆上，且在在单位圆上的根只能是单重根。定理4若线性k步法(5.1)收敛，则必稳定。定理5若线性k步法(5.1)相容且稳定，则必收敛稳定。7.4绝对稳定性定义4对给定的μ,h,若特征方程(7.24)的所有根的模都1，则称线性k步法(7.23)绝对稳定。若对所有μh∈(α,β),(7.23)都绝对稳定，则说(α,β)为绝对稳定区间。εεεαβρσ