我們現在簡要討論如何求解初值問題。有關更多資訊，請參見 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)
