我们现在简要讨论如何求解初值问题。有关更多信息，请参见 Bradie (Chap. 7)[ 1] [ 2]
为了有效地计算计算机上的微分方程解，方便地在有限数量的指定点进行计算，然后在这些点之间插值。对于许多计算来说，使用网格很方便，网格上的点彼此相隔等距。
在本节的其余部分 h {\displaystyle h} h {\displaystyle h} t n {\displaystyle t_{n}} n th {\displaystyle n^{\text{th}}} y n {\displaystyle y^{n}} n th {\displaystyle n^{\text{th}}} y n ≡ y ( t = t n ) {\displaystyle y^{n}\equiv y(t=t^{n})} h {\displaystyle h} t {\displaystyle t} h = t n + 1 − t n {\displaystyle h=t^{n+1}-t^{n}} f ( t , y ) {\displaystyle f(t,y)} y ( t ) {\displaystyle y(t)}
d y d t = f ( t , y ) y ( t 0 ) = y 0 . {\displaystyle {\frac {dy}{dt}}=f(t,y)\qquad y(t^{0})=y^{0}.}
微分方程可以用有限差分来逼近，
y n + 1 − y n h = f ( t n , y n ) . {\displaystyle {\frac {y^{n+1}-y^{n}}{h}}=f(t^{n},y^{n}).}
现在我们要做的就是用代数方法求解 y n + 1 {\displaystyle y^{n+1}}
y n + 1 = y n + h f ( t n , y n ) (Forward Euler/Explicit method) {\displaystyle y^{n+1}=y^{n}+hf(t^{n},y^{n})\qquad {\text{(Forward Euler/Explicit method)}}}
如果我们想要计算 d y d t {\displaystyle {\frac {dy}{dt}}} t 0 {\displaystyle t^{0}} 2 t n + 1 {\displaystyle t^{n+1}} y ( t 0 ) {\displaystyle y(t^{0})} y n + 1 {\displaystyle y^{n+1}}
让我们以公式 1 y ( t 0 ) = 1 {\displaystyle y(t^{0})=1} t 0 = 0 {\displaystyle t^{0}=0} h = 1 {\displaystyle h=1} h = .1 {\displaystyle h=.1}
公式 1 毫无疑问，与 h = 1 {\displaystyle h=1} h = .1 {\displaystyle h=.1} h = 0.1 {\displaystyle h=0.1} y ( t ) {\displaystyle y(t)} h = .1 {\displaystyle h=.1} [ 脚注 1]
y ( t n + h ) = y ( t n ) + h d y d t | t n + O ( h 2 ) , {\displaystyle y(t^{n}+h)=y(t^{n})+h\left.{\frac {dy}{dt}}\right|_{t^{n}}+{\text{O}}(h^{2}),}
其中，比 O ( h 2 ) {\displaystyle (h^{2})} 2
y n + h d y d t | t n + O ( h 2 ) = y n + h f ( t n , y n ) {\displaystyle y^{n}+h\left.{\frac {dy}{dt}}\right|_{t^{n}}+{\text{O}}(h^{2})=y^{n}+hf(t^{n},y^{n})}
在消去项并除以 h {\displaystyle h}
d y d t | t n + O ( h ) = f ( t n , y n ) , {\displaystyle \left.{\frac {dy}{dt}}\right|_{t^{n}}+{\text{O}}(h)=f(t^{n},y^{n}),}
由此可以清楚地看出，该方法的精度随步长线性变化，因此它是 一阶精度的。
可以通过使用向后差分商来逼近导数，从而获得向前欧拉法的变体。使用等式 1
y n − y n − 1 h ≈ f ( t n , y n ) {\displaystyle {\frac {y^{n}-y^{n-1}}{h}}\approx f(t^{n},y^{n})}
y n = y n − 1 + h f ( t n , y n ) . {\displaystyle y^{n}=y^{n-1}+hf(t^{n},y^{n}).}
将索引从 n {\displaystyle n} n + 1 {\displaystyle n+1}
y n + 1 = y n + h f ( t n + 1 , y n + 1 ) (Backwards Euler/Implicit method) {\displaystyle y^{n+1}=y^{n}+hf(t^{n+1},y^{n+1})\qquad {\text{(Backwards Euler/Implicit method)}}}
注意 y n + 1 {\displaystyle y^{n+1}} y n + 1 {\displaystyle y^{n+1}} y n + 1 {\displaystyle y^{n+1}} f {\displaystyle f} f {\displaystyle f} y 2 {\displaystyle y^{2}}
通过对向前欧拉方法和向后欧拉方法进行平均，我们可以找到 Crank-Nicolson 方法
y n + 1 − y n h = 1 2 f ( t n + 1 , y n + 1 ) + 1 2 f ( t n , y n ) {\displaystyle {\frac {y^{n+1}-y^{n}}{h}}={\frac {1}{2}}f(t^{n+1},y^{n+1})+{\frac {1}{2}}f(t^{n},y^{n})}
重新排列，我们得到：
y n + 1 = y n + h 2 [ f ( t n + 1 , y n + 1 ) + f ( t n , y n ) ] (Crank-Nicolson) {\displaystyle y^{n+1}=y^{n}+{\frac {h}{2}}\left[f(t^{n+1},y^{n+1})+f(t^{n},y^{n})\right]\qquad {\text{(Crank-Nicolson)}}}
再次注意 y n + 1 {\displaystyle y^{n+1}} y n + 1 {\displaystyle y^{n+1}} y n + 1 {\displaystyle y^{n+1}}
让我们看一下以下常微分方程
d y d t = − λ y ( t ) {\displaystyle {\frac {dy}{dt}}=-\lambda y(t)}
其中 λ {\displaystyle \lambda } y ( t 0 ) = 1 {\displaystyle y(t^{0})=1} t 0 = 0 {\displaystyle t^{0}=0}
y n + 1 = y n − λ h y n (Forward Euler) {\displaystyle y^{n+1}=y^{n}-\lambda hy^{n}\qquad {\text{(Forward Euler)}}}
y n + 1 = y n ( 1 + λ h ) (Backward Euler) {\displaystyle y^{n+1}={\frac {y^{n}}{(1+\lambda h)}}\qquad {\text{(Backward Euler)}}}
y n + 1 = y n ( 2 − λ h 2 + λ h ) (Crank-Nicolson) {\displaystyle y^{n+1}=y^{n}\left({\frac {2-\lambda h}{2+\lambda h}}\right)\qquad {\text{(Crank-Nicolson)}}}
而精确解由下式给出
y ( t ) = e − λ t (Exact solution) {\displaystyle y(t)=e^{-\lambda t}\qquad {\text{(Exact solution)}}}
方程 2 上图显示了两种方法如何收敛到解，但前向欧拉解对于所选时间步长是非稳定的。下面的清单是一个 Matlab 程序，你可以用它来改变 λ {\displaystyle \lambda }
Matlab 程序
代码下载
% A program to demonstrate instability of timestepping methods
% when the timestep is inappropriately choosen.
%Differential equation: y'(t)=-y(t) y(t_0)=y_0
%Initial Condition, y(t_0)=1 where t_0=0
clear all ; clc ; clf ;
%Grid
h = .1 ;
tmax = 4 ;
Npoints = tmax / h ;
lambda = .1 ;
%Initial Data
y0 = 1 ;
t_0 = 0 ;
t ( 1 )= t_0 ;
y_be ( 1 )= y0 ;
y_fe ( 1 )= y0 ;
y_imr ( 1 )= y0 ;
for n = 1 : Npoints
%Forward Euler
y_fe ( n + 1 )= y_fe ( n ) - lambda * h * y_fe ( n );
%Backwards Euler
y_be ( n + 1 )= y_be ( n ) / ( 1 + lambda * h );
%Crank Nicolson
y_imr ( n + 1 )= y_imr ( n ) * ( 2 - lambda * h ) / ( 2 + lambda * h )
t ( n + 1 )= t ( n ) + h ;
end
%Exact Solution
tt =[ 0 : .001 : tmax ];
exact = exp ( - lambda * tt );
%Plot
figure ( 1 ); clf ; plot ( tt , exact , 'r-' , t , y_fe , 'b:' , t , y_be , 'g--' , t , y_imr , 'k-.' );
xlabel time ; ylabel y ;
legend ( 'Exact' , 'Forward Euler' , 'Backward Euler' , ...
'Crank Nicolson' );
一个 Python 程序，用来演示不同时间步进方法的非稳定性。
代码下载
#!/usr/bin/env python
"""
A program to demonstrate instability of timestepping methods#
when the timestep is inappropriately choosen.################
"""
from math import exp
import matplotlib.pyplot as plt
import numpy
#Differential equation: y'(t)=-l*y(t) y(t_0)=y_0
#Initial Condition, y(t_0)=1 where t_0=0
# Definition of the Grid
h = 0.1 # Time step size
t0 = 0 # Initial value
tmax = 4 # Value to be computed y(tmax)
Npoints = int (( tmax - t0 ) / h ) # Number of points in the Grid
t = [ t0 ]
# Initial data
l = 0.1
y0 = 1 # Initial condition y(t0)=y0
y_be = [ y0 ] # Variables holding the value given by the Backward Euler Iteration
y_fe = [ y0 ] # Variables holding the value given by the Forward Euler Iteration
y_imr = [ y0 ] # Variables holding the value given by the Midpoint Rule Iteration
for i in xrange ( Npoints ):
y_fe . append ( y_fe [ - 1 ] * ( 1 - l * h ))
y_be . append ( y_be [ - 1 ] / ( 1 + l * h ))
y_imr . append ( y_imr [ - 1 ] * ( 2 - l * h ) / ( 2 + l * h ))
t . append ( t [ - 1 ] + h )
print "Exact Value: y( %d )= %f " % ( tmax , exp ( - 4 ))
print "Backward Euler Value: y( %d )= %f " % ( tmax , y_be [ - 1 ])
print "Forward Euler Value: y( %d )= %f " % ( tmax , y_fe [ - 1 ])
print "Midpoint Rule Value: y( %d )= %f " % ( tmax , y_imr [ - 1 ])
# Exact Solution
tt = numpy . arange ( 0 , tmax , 0.001 )
exact = numpy . exp ( - l * tt )
# Plot
plt . figure ()
plt . plot ( tt , exact , 'r-' , t , y_fe , 'b:' , t , y_be , 'g--' , t , y_imr , 'k-.' );
plt . xlabel ( 'time' )
plt . ylabel ( 'y' )
plt . legend (( 'Exact' , 'Forward Euler' , 'Backward Euler' ,
'Implicit Midpoint Rule' ))
plt . show ()
y ′ = − λ y {\displaystyle y'=-\lambda y} [ edit | edit source ] 我们现在尝试理解这些观察结果，以便我们有一些指导方针来设计稳定的数值方法。如果数值解可以被不依赖步长的函数所限制，那么具有有界解的初值问题的数值解是稳定的 。通常使用两种方法来理解稳定性。第一种方法是线性化稳定性，它涉及计算线性系统的特征值，以查看小的扰动是增长还是衰减。第二种方法是计算与微分方程相关的能量类数量，并检查它是否保持有界。
我们假设 λ ≥ 0 {\displaystyle \lambda \geq 0}
y n + 1 − y n h = − λ y n {\displaystyle {\frac {y^{n+1}-y^{n}}{h}}=-\lambda y^{n}}
y n + 1 = ( 1 − λ h ) y n {\displaystyle y^{n+1}=(1-\lambda h)y^{n}}
⇒ | y n + 1 | ≥ | ( 1 − λ h ) | | y n | if | ( 1 − λ h ) | > 1 {\displaystyle \Rightarrow |y^{n+1}|\geq |(1-\lambda h)||y^{n}|\quad {\text{ if }}|(1-\lambda h)|>1}
⇒ | y n + 1 | ≤ | ( 1 − λ h ) | | y n | if | ( 1 − λ h ) | < 1. {\displaystyle \Rightarrow |y^{n+1}|\leq |(1-\lambda h)||y^{n}|\quad {\text{ if }}|(1-\lambda h)|<1.}
我们可以对后向欧拉方法进行类似的计算，得到
y n + 1 − y n h = − λ y n + 1 {\displaystyle {\frac {y^{n+1}-y^{n}}{h}}=-\lambda y^{n+1}}
y n + 1 = y n 1 + λ h {\displaystyle y^{n+1}={\frac {y^{n}}{1+\lambda h}}}
⇒ | y n + 1 | ≤ | y n 1 + λ h | ≤ | y n | since | 1 1 + λ h | < 1. {\displaystyle \Rightarrow |y^{n+1}|\leq \left|{\frac {y^{n}}{1+\lambda h}}\right|\leq |y^{n}|\quad {\text{ since }}\left|{\frac {1}{1+\lambda h}}\right|<1.}
因此，后向欧拉方法是无条件稳定的，而前向欧拉方法不是。我们留下了Crank-Nicolson方法的分析作为练习。
第二种方法，常用于显示偏微分方程的稳定性，是寻找能量式的量并证明该量可以限制解，并防止解变得过大或过小。通常，该量被选为非负的，那么只需要推导出一个上限。我们概述了如何对常微分方程进行此操作，以便我们可以在观察偏微分方程时使用相同的思想。回想一下，前向欧拉算法由下式给出
y n + 1 − y n h = − λ y n . {\displaystyle {\frac {y^{n+1}-y^{n}}{h}}=-\lambda y^{n}.}
用 y n + 1 {\displaystyle y^{n+1}}
( y n + 1 ) 2 = ( 1 − h λ ) y n y n + 1 . {\displaystyle (y^{n+1})^{2}=(1-h\lambda )y^{n}y^{n+1}.}
现在要获得关于 | y n + 1 | {\displaystyle |y^{n+1}|} | y n | {\displaystyle |y^{n}|}
( a − b ) 2 ≥ 0 ⇒ a 2 + b 2 ≥ 2 a b ⇒ ( y n + 1 ) 2 + ( y n ) 2 2 ≥ y n y n + 1 . {\displaystyle (a-b)^{2}\geq 0\Rightarrow a^{2}+b^{2}\geq 2ab\Rightarrow {\frac {(y^{n+1})^{2}+(y^{n})^{2}}{2}}\geq y^{n}y^{n+1}.}
因此，如果
( 1 − h λ ) > 0 {\displaystyle (1-h\lambda )>0}
是
( y n + 1 ) 2 ≤ ( 1 − h λ ) ( y n + 1 ) 2 + ( y n ) 2 2 {\displaystyle (y^{n+1})^{2}\leq (1-h\lambda ){\frac {(y^{n+1})^{2}+(y^{n})^{2}}{2}}}
( y n + 1 ) 2 1 + h λ 2 ≤ 1 − h λ 2 ( y n ) 2 {\displaystyle (y^{n+1})^{2}{\frac {1+h\lambda }{2}}\leq {\frac {1-h\lambda }{2}}(y^{n})^{2}}
( y n + 1 ) 2 ≤ 1 − h λ 1 + h λ ( y n ) 2 , {\displaystyle (y^{n+1})^{2}\leq {\frac {1-h\lambda }{1+h\lambda }}(y^{n})^{2},}
因此，如果 1 − h λ > 0 {\displaystyle 1-h\lambda >0} 0 < 1 − h λ 1 + h λ < 1 {\displaystyle 0<{\frac {1-h\lambda }{1+h\lambda }}<1} λ {\displaystyle \lambda } h {\displaystyle h}
稳定性不是数值方法逼近初值问题解的唯一要求。我们还希望证明，随着时间步长的减小，数值逼近会变得更好。对于前向欧拉方法，我们有
y n + 1 − y n h = − λ y n {\displaystyle {\frac {y^{n+1}-y^{n}}{h}}=-\lambda y^{n}}
现在如果
y n = y ( t ) {\displaystyle y^{n}=y(t)}
y n + 1 = y ( t + h ) {\displaystyle y^{n+1}=y(t+h)}
那么[ 脚注 2]
y ( t + h ) = y ( t ) + h d y d t + O ( h 2 ) {\displaystyle y(t+h)=y(t)+h{\frac {\mathrm {d} y}{\mathrm {d} t}}+O(h^{2})}
所以
y n + 1 − y n h + λ y n {\displaystyle {\frac {y^{n+1}-y^{n}}{h}}+\lambda y^{n}}
= y ( t + h ) − y ( t ) h + λ y ( t ) {\displaystyle ={\frac {y(t+h)-y(t)}{h}}+\lambda y(t)}
= d y d t + O ( h ) + λ y ( t ) {\displaystyle ={\frac {\mathrm {d} y}{\mathrm {d} t}}+O(h)+\lambda y(t)}
= O ( h ) . {\displaystyle =O(h).}
我们可以做类似的计算来证明Crank-Nicolson方法是二阶的。在这种情况下，我们使用围绕 y ( t + h / 2 ) {\displaystyle y(t+h/2)}
y n + 1 − y n h = − λ y n + 1 + y n 2 {\displaystyle {\frac {y^{n+1}-y^{n}}{h}}=-\lambda {\frac {y^{n+1}+y^{n}}{2}}}
所以
y n + 1 = y ( t + h ) = y ( t + h / 2 ) + ( h / 2 ) d y d t + ( h / 2 ) 2 1 2 d 2 y d t 2 + O ( h 3 ) {\displaystyle y^{n+1}=y(t+h)=y(t+h/2)+(h/2){\frac {\mathrm {d} y}{\mathrm {d} t}}+(h/2)^{2}{\frac {1}{2}}{\frac {\mathrm {d} ^{2}y}{\mathrm {d} t^{2}}}+O(h^{3})}
y n = y ( t ) = y ( t + h / 2 ) − ( h / 2 ) d y d t + ( h / 2 ) 2 1 2 d 2 y d t 2 + O ( h 3 ) {\displaystyle y^{n}=y(t)=y(t+h/2)-(h/2){\frac {\mathrm {d} y}{\mathrm {d} t}}+(h/2)^{2}{\frac {1}{2}}{\frac {\mathrm {d} ^{2}y}{\mathrm {d} t^{2}}}+O(h^{3})}
因此
y n + 1 − y n h + λ y n + 1 + y n 2 {\displaystyle {\frac {y^{n+1}-y^{n}}{h}}+\lambda {\frac {y^{n+1}+y^{n}}{2}}}
= d y d t + O ( h 2 ) + λ [ y ( t + h / 2 ) + O ( h 2 ) ] {\displaystyle ={\frac {\mathrm {d} y}{\mathrm {d} t}}+O(h^{2})+\lambda \left[y(t+h/2)+O(h^{2})\right]}
= O ( h 2 ) . {\displaystyle =O(h^{2}).}
因此，这是一个二阶方法。
确定隐式中点规则对于常微分方程 d y d t = − λ y {\displaystyle {\frac {\mathrm {d} y}{\mathrm {d} t}}=-\lambda y} λ {\displaystyle \lambda } h {\displaystyle h} λ {\displaystyle \lambda } h {\displaystyle h}
证明后向欧拉方法是一阶方法。
↑ 泰勒展开式的推导可以在大多数微积分书籍中找到。 ↑ 我们将使用大O符号来表示，如果 f O ( h ) {\displaystyle f~O(h)} h → 0 {\displaystyle h\rightarrow 0} | f h | < C {\displaystyle \left|{\frac {f}{h}}\right|<C} C {\displaystyle C}
↑ Bradie (2006) ↑ Boyce and DiPrima (2010)
Boyce, W.E.; DiPrima, R.C. (2010). Elementary Differential Equations and Boundary Value Problems . Wiley.
Bradie, B. (2006). A Friendly Introduction to Numerical Analysis . Pearson.
![{\displaystyle y^{n+1}=y^{n}+{\frac {h}{2}}\left[f(t^{n+1},y^{n+1})+f(t^{n},y^{n})\right]\qquad {\text{(Crank-Nicolson)}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/55edaf19a5fc35c43ddb0ed5d0039f1db92c543e)
![{\displaystyle ={\frac {\mathrm {d} y}{\mathrm {d} t}}+O(h^{2})+\lambda \left[y(t+h/2)+O(h^{2})\right]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/6cb2f3771bc16131e61c5f6570e243d0321ef0eb)
