Transcript 8．3 常微分方程式の数値解放

８．３ 常微分方程式の数値解法
数値解法では，解析的に積分できないような問題
を扱うことが多い
オイラーってどんな人？
数学者・物理学者であり、天体物理学者。
スイスのバーゼル生まれ
現ロシアのサンクトペテルブルクで死んだ。
1707/4/15-1783/9/18
Leonhard Euler
８．３．１ オイラー法
（１）前進差分によるオイラー法
dy
微分方程式
 g ( t , y ),
dt
y (0 )  C1
［例］前進差分近似
y (t   t )  y (t )
t
 g ( t , y ( t ))  y ( t   t )  y ( t )  g ( t , y ( t ))  t
次のような繰返し計算で近似する
y
( t1 , y 1 )
y (0)  C1
(t 0 , y 0 )
y (  t )  y ( 0 )  g ( 0 , y ( 0 ))  t
g (t 0 , y 0 )  t
y ( 2  t )  y (  t )  g (  t , y (  t ))  t

t
0
t
あるいは，
差分近似が微分の近似であることから
（本来は， t は限りなく０に近い値）
dy
dt
 lim
y (t   t )  y (t )
t  0
t

y (t   t )  y (t )
t
y (0)  C1
y (  t )  y ( 0 )  g (  t , y ( 0 ))  t
y ( 2  t )  y (  t )  g ( 2  t , y (  t ))  t

dy
（例）
 ( t  1) ydt
Excelでオイラー法（１）
dt
もちろん，この式は次のようにして厳密解を求めることができるが，誤差評価
のためにこの例を求める。
1
dy  ( t  1) dt 
y

1
y
dy 
 ( t  1) dt
 log y 
t
2
tC
2
［Excelの定義
例］
tとして現時刻をとる定義
（例）
dy
 ( t  1) ydt
dt
Excelでオイラー法（２）
グラフを描いてみると
tとして現時刻をとる定義のほうが誤差が少ない
（２）修正オイラー法（修正法１）
データ点における微分係数で増分を計算し，
増分から得られた座標点から更に増分を計算し，
両増分の平均値を最終的な増分とする。
y
k1  g (t , y )  t
( t1 , y 1 )
k 2  g (t   t , y  k1 )  t
y (t   t )  y (t ) 
1
2
k2
( k1  k 2 )
(t 0 , y 0 )
k1
t
0
t
Excelによる修正法１
式定義
結果グラフ
（３）修正オイラー法（修正法２）
データ点における微分係数で増分を計算し，
増分の１／２から得られた座標点から増分を計算し，
この増分を最終的な増分とする。
y
k1  g (t , y )  t
k 2  g (t 
t
2
,y
k1
( t1 
t
2
, y1 
k1
)
2
)t
2
k2
(t 0 , y 0 )
y (t   t )  y (t )  k 2
k1
t
0
t
Excelによる修正法２
式定義
結果グラフ
VBAによる定義
①データ宣言，データ設定および結果設定用シートへのX値設定
Private DX As Double
Private Y0 As Double
Private XST As Double
Private Y(11) As Double
Private N As Integer
Sub データ設定()
N = 11: Y0 = 1: DX = 0.1: XST = 0
End Sub
Sub X値設定()
With Worksheets("Sheet4")
X = XST
For i = 1 To N
.Cells(i + 1, 1) = X: X = X + DX
Next
End With
End Sub
VBAによる定義
②結果設定とボタンのClickイベントハンドラ
Sub 結果設定(ID)
With Worksheets("Sheet4")
For i = 1 To N
.Cells(i + 1, ID) = Y(i)
Next
End With
End Sub
Sub ボタン1_Click()
データ設定
X値設定
オイラー法前進差分 Y, Y0, XST, DX, N
結果設定 2
修正オイラー法1 Y, Y0, XST, DX, N
結果設定 3
修正オイラー法2 Y, Y0, XST, DX, N
結果設定 4
厳密解 Y, Y0, XST, DX, N
結果設定 5
End Sub
VBAによる定義
③関数値，オイラー法，修正法１
Function g(X, Y) As Double
g = (X + 1) * Y
End Function
Sub オイラー法前進差分(Y, Y0, XST, DX, N)
Y(1) = Y0: X = XST
For i = 2 To N
X = X + DX
Y(i) = Y(i - 1) + g(X, Y(i - 1)) * DX
Next
End Sub
Sub 修正オイラー法1(Y, Y0, XST, DX, N)
Y(1) = Y0: X = XST
For i = 2 To N
K1 = DX * g(X, Y(i - 1))
K2 = DX * g(X + DX, Y(i - 1) + K1)
X = X + DX
Y(i) = Y(i - 1) + (K1 + K2) / 2
Next
End Sub
VBAによる定義
④修正法２，厳密解設定
Sub 修正オイラー法2(Y, Y0, XST, DX, N)
Y(1) = Y0: X = XST
For i = 2 To N
K1 = DX * g(X, Y(i - 1))
K2 = DX * g(X + DX / 2, Y(i - 1) + K1 / 2)
Y(i) = Y(i - 1) + K2
X = X + DX
Next
End Sub
Sub 厳密解(Y, Y0, XST, DX, N)
X = XST
For i = 1 To N
Y(i) = Exp(X * X / 2 + X)
X = X + DX
Next
End Sub
VBAによる定義による結果
結果（Excelによる結果より誤差が少ない）
（４）修正オイラー法（一般化）
以下の係数に適当な値を代入する。修正法１，２はその代表的な例。
k1  g (t , y )  t
k 2  g (t     t , y    k1 )  t
y ( t   t )  y ( t )  ( c1 k 2  c 2 k 2 )
［注］以下，式の展開による導出を示すが，
式の展開自体を覚える必要はない。
どのような流れで導出されているかを理解し，
最終的な結果を確認すること。
式の展開自体が必要なときは，資料を参照すれば良い。
確認
k 2 をテーラ展開すると
k 2  g (t     t , y    k1 )  t  g (t , y )     t
 g (t , y )     t
dg
   g (t , y )  t
dt
dg
dg
   k1
dt
dg
 O (t , y )
2
2
dy
 O (t , y )
2
2
dy
第３式に代入すると
y ( t   t )  y ( t )  ( c1 k 2  c 2 k 2 )

dg
dg
2 
 y ( t )  c1  t  g ( t , y )  c 2  t  g ( t , y )     t
   g (t , y )  t
 O (t )
dt
dy


 y ( t )  ( c1  c 2 )  t  g ( t , y )    c 2  t
テーラ展開の式と比較すると
修正法１ c1  c 2 
1
2
,
2
dg
dt
   g (t , y ) c 2  t
( c1  c 2 )  1,
   1
  c2 
1
2
dg
,
2
修正法２ c1  0 ,
 O (t )
3
dy
  c2 
c 2  1,
1
2
  
1
2
８．３．２ ルンゲ・クッタ法（Runge-Kutta method）
（１）２次のルンゲ・クッタ法(その１）
k 1  g ( t j , x j ),
k 2  g ( t j  b1 h , x j  b 2 hk 1 ),
k 2  g ( t j  b1 h , x j  b 2 hk 1 )  g ( t j , x j )  b1 h
 g ( t j , x j )  b1 h
 g (t j , x j )
t
 b 2 hg ( t j , x j )
x j  1  x j  h ( a1 k 1  a 2 k 2 ) とおくと
 g (t j , x j )
t
 g (t j , x j )
x
 b 2 hk 1
 g (t j , x j )
x
 O (h )
2
 O (h )
2


 g (t j , x j )
 g (t j , x j )

2 
 x j  1  x j  h  a 1 g ( t j , x j )  a 2  g ( t j , x j )  b1 h
 b 2 hg ( t j , x j )
 O ( h )  
t
x



 g (t j , x j )
 g (t j , x j ) 

3
 x j  h ( a 1  a 2 ) g ( t j , x j )  h  a 2 b1
 a 2 b2 g (t j , x j )
  O (h )
t
x


2
（１）２次のルンゲ・クッタ法（その２）
2
定義から
d x
dt
x j 1
2

dg ( t , x )
dt

 g (t , x )
t

 g (t , x )
x
g (t , x )
であるからテーラ展開すると
2

h   g (t j , x j )  g (t j , x j )
3
 x j  hg ( t j , x j ) 

g (t j , x j )   O ( h )

2 
t
x

係数を比較して
a1  a 2  1,
これらを満足する値として a1  a 2 
1
2
,
a 2 b1 
1
2
,
a 2 b2 
1
2
b1  b 2  1 を選択して
k1  g (t j , x j )
k 2  g ( t j  h , x j  hk 1 )
x j 1  x j 
h
2
k1  k 2 
この式はホイン（Heun）の２次公式とも呼ばれる。
y j  1  y j  hg ( t j , y j )
（２）３次のルンゲ・クッタ法
において
k 1  g ( t j , y j ), k 2  g ( t j  b1 h , y j  b 2 hk 1 ), k 3  g ( t j  c1 h , y j  c 2 hk 1  c 3 hk 2 )
y j 1  y j  h ( a1 k 1  a 2 k 2  a 3 k 3 )
とおくと，２次の場合と同様，係数比較によって a1 , a 2 , a 3 , b1 , b2 , c1 , c 2 , c 3
を決めることができる。
k1  g (t j , x j )
k 2  g (t j 
h
2
, xj 
h
2
k1 )
k 3  g ( t j  h , x j  hk 1  2 hk 2 )
x j 1  x j 
h
6
k1  4 k 2  k 3 
この式はクッタ（Kutta）の公式と呼ばれる。
y j  1  y j  hg ( t j , y j )
k 1  g ( t j , y j ),
（３）４次のルンゲ・クッタ法
において
k 2  g ( t j  b1 h , y j  b 2 hk 1 ),
k 3  g ( t j  c1 h , y j  c 2 hk 2 ),
k 4  g ( t j  d 1 h , y j  d 2 hk 2 ),
y j 1  y j  h ( a1 k 1  a 2 k 2  a 3 k 3  a 4 k 4 )
とおくと，２次の場合と同様，係数比較によって a1 , a 2 , a 3 , b1 , b2 , c1 , c 2 , d 1 , d 2
を決めることができる。
k1  g (t j , x j )
k 2  g (t j 
h
k 3  g (t j 
h
2
2
, xj 
h
, xj 
h
2
2
k1 )
k2 )
k 4  g ( t j  h , x j  hk 3 )
x j 1  x j 
h
6
k1  2 k 2  2 k 3  k 4 
この式はルンゲ・クッタの式と呼ばれる。
（４）ルンゲ・クッタ・ギルの公式
（Runge-Kutta-Gill method）
４次のルンゲ・クッタ型の公式の一つとして，次の公式がある。
k 1  g t j , x j 
h
h 

k 2  g  t j  , x j  k1 
2
2 



h
1 
1 
1


k3  g  t j  , x j   
 hk 1   1 
 hk 2 
2
2
2
2





1
1 

k 4  g  t j  h , x j 
hk 2   1 
 hk 3 
2
2



x j 1

h
1 
1 



 x j   k1  2  1 
k 2  21 
 k 3  k 4 
6
2
2



ルンゲとクッタ
• カール・ルンゲ(Runge,Carl David Tolme) ドイツ
の数学者で物理学者。 1856/08/30-1927/01/03
• ウイルヘルム・クッタ(Wilhelm Kutta)ドイツの数
学者でかつ物理学者。1867~1944。渦による揚
力の発生クッタ-ジェコフスキーの揚力理論で有
名。
VBAによる定義
①データ宣言，データ設定および結果設定用シートへのX値設定
Private DX As Double
Private Y0 As Double
Private XST As Double
Private Y(11) As Double
Private N As Integer
Sub データ設定()
N = 11: Y0 = 1: DX = 0.1: XST = 0
End Sub
Sub X値設定()
With Worksheets("Sheet4")
X = XST
For i = 1 To N
.Cells(i + 1, 1) = X: X = X + DX
Next
End With
End Sub
VBAによる定義
②結果設定とボタンのClickイベントハンドラ
Sub 結果設定(ID)
With Worksheets("Sheet4")
For i = 1 To N
.Cells(i + 1, ID) = Y(i)
Next
End With
End Sub
Sub ボタン1_Click()
データ設定
X値設定
ルンゲクッタ Y, Y0, XST, DX, N
結果設定 2
ルンゲクッタギル Y, Y0, XST, DX, N
結果設定 3
厳密解 Y, Y0, XST, DX, N
結果設定 4
End Sub
VBAによる定義
③関数値，４次のルンゲ・クッタ法
Function g(X, Y) As Double
g = (X + 1) * Y
End Function
Sub ルンゲクッタ(Y, Y0, XST, DX, N)
Y(1) = Y0: X = XST
For i = 2 To N
k1 = g(X, Y(i - 1))
k2 = g(X + DX / 2, Y(i - 1) + DX * k1 / 2)
k3 = g(X + DX / 2, Y(i - 1) + DX * k2 / 2)
k4 = g(X + DX, Y(i - 1) + DX * k3)
Y(i) = Y(i - 1) + (k1 + 2 * k2 + 2 * k3 + k4) * DX / 6
X = X + DX
Next
End Sub
VBAによる定義
④ルンゲ・クッタ・ギル法
Sub ルンゲクッタギル(Y, Y0, XST, DX, N)
Y(1) = Y0: X = XST
A1 = -(1 / 2 - 1 / Sqr(2))
A2 = (1 - 1 / Sqr(2))
B1 = -1 / Sqr(2)
B2 = 1 + 1 / Sqr(2)
C1 = 2 * (1 - 1 / Sqr(2))
C2 = 2 * (1 + 1 / Sqr(2))
For i = 2 To N
k1 = g(X, Y(i - 1))
k2 = g(X + DX / 2, Y(i - 1) + DX * k1 / 2)
k3 = g(X + DX / 2, Y(i - 1) + A1 * DX * k1 + A2 * DX * k2)
k4 = g(X + DX, Y(i - 1) + B1 * DX * k2 + B2 * DX * k3)
Y(i) = Y(i - 1) + (k1 + C1 * k2 + C2 * k3 + k4) * DX / 6
X = X + DX
Next
End Sub
VBAによる定義
⑤厳密解の設定
Sub 厳密解(Y, Y0, XST, DX, N)
X = XST
For i = 1 To N
Y(i) = Exp(X * X / 2 + X)
X = X + DX
Next
End Sub
VBAによる定義
④修正法２，厳密解設定
Sub 修正オイラー法2(Y, Y0, XST, DX, N)
Y(1) = Y0: X = XST
For i = 2 To N
K1 = DX * g(X, Y(i - 1))
K2 = DX * g(X + DX / 2, Y(i - 1) + K1 / 2)
Y(i) = Y(i - 1) + K2
X = X + DX
Next
End Sub
Sub 厳密解(Y, Y0, XST, DX, N)
X = XST
For i = 1 To N
Y(i) = Exp(X * X / 2 + X)
X = X + DX
Next
End Sub
ルンゲクッタ法による結果
非常に良い結果が得られていることに注意
８．３．３ 連立微分方程式
（１）手順
次の連立１次微分方程式を考える。
dx
 f (t , x )
dt
 x1 
 
x2

,
x 
  
 
 xn 
 f 1 ( t , x1 ,  , x n ) 


f 2 ( t , x1 ,  , x n )

f  





f
(
t
,
x
,

,
x
)
1
n 
 n
①増分ベクトル k 1 を求める。
k1  f (t , x )
 k 11   f 1 ( t , x1 ,  , x n ) 
  

k 12
f 2 ( t , x1 ,  , x n )
 

   


  

k
f
(
t
,
x
,

,
x
)
1
n 
 1n   n
②増分 k 2 ～ k 4 を求める。
k 2  f (t 
h
k 3  f (t 
h
,x
2
h k1
)
2
,x
2
hk 2
)
2
k 4  f (t  h , x  h k 3 )
③ x ( t  h ) を求める。
x (t  h )  x (t ) 
h
6
k 1  2 k 2  2 k 3  k 4 
（２）プログラム例
次の連立１次微分方程式を解くプログラムを例とする。
dy 1
dx
 x  y2 ,
y1 ( 0 )  0 ,
dy 2
 x
dx
y 2 (0)  0
①データ宣言，データ設定および結果設定用シートへのX値設定
Private DX As Double
Private Y0() As Double
Private XST As Double
Private Y(11, 2) As Double
Private N As Integer
Private M As Integer
Sub データ設定()
N = 11: DX = 0.1: XST = 0
M = 2: ReDim Y0(M)
For k = 1 To M
Y0(k) = 0
Next
End Sub
Sub X値設定()
With Worksheets("Sheet2")
X = XST
For i = 1 To N
.Cells(i + 1, 1) = X: X = X + DX
Next
End With
End Sub
②結果設定およびボタンのClickイベントハンドラ
Sub 結果設定()
With Worksheets("Sheet2")
For i = 1 To N
For k = 1 To M
.Cells(i + 1, k + 1) = Y(i, k)
Next
Next
End With
End Sub
Sub Sheet2_ボタン1_Click()
データ設定
X値設定
連立ルンゲクッタ Y, Y0, XST, DX, N, M
結果設定
End Sub
③関数設定および処理本体（その１）
Function g(ID, X, Y) As Double
If ID = 1 Then
g = X + Y(2)
Else
g = X
End If
End Function
Sub 連立ルンゲクッタ(Y, Y0, XST, DX, N, M)
Dim K1() As Double: ReDim K1(M)
Dim K2() As Double: ReDim K2(M)
Dim K3() As Double: ReDim K3(M)
Dim K4() As Double: ReDim K4(M)
Dim YY() As Double: ReDim YY(M)
For k = 1 To M
Y(1, k) = Y0(k)
Next
⑤処理本体（その３）
X = XST
For i = 2
For k =
YY(k)
Next
For k =
K1(k)
Next
For k =
YY(k)
Next
For k =
K2(k)
Next
For k =
YY(k)
Next
For k =
K3(k)
Next
To N
1 To M
= Y(i - 1, k)
1 To M
= g(k, X, YY)
1 To M
= Y(i - 1, k) + DX * K1(k) / 2
1 To M
= g(k, X + DX / 2, YY)
1 To M
= Y(i - 1, k) + DX * K2(k) / 2
1 To M
= g(k, X + DX / 2, YY)
⑤処理本体（その４）
For k = 1
YY(k) =
Next
For k = 1
K4(k) =
Next
For k = 1
Y(i, k)
Next
X = X + DX
Next
End Sub
To M
Y(i - 1, k) + DX * K3(k)
To M
g(k, X + DX, YY)
To M
= Y(i - 1, k) + _
(K1(k) + 2 * K2(k) + 2 * K3(k) + K4(k)) * DX / 6
８．３．４ 高階の常微分方程式
（１）考え方
次の微分方程式を考える。
n
a0
d x
dt
n
 a1
d
n 1
x
n 1
dt
   a n 1
dx
dt
 an x  0
各次数の時間微分を別々の変数と考えると
x  x1 ,
dx 1
dx n
1
dt

dt
a0
 x2 ,
dx 2
dt
 x3 ,
,
dx n 1
dt
 a 1 x n  a 2 x n  1    a n x1 
という連立１次微分方程式の形式であるため，
ルンゲ・クッタ法を適用することができる。
 xn ,
［例題］
2
m
d x
dt
2
  kx
解析的に解けるが，
厳密解が分かっているので，
誤差検証のために
あえてこの例題を取り上げる。
（連立ルンゲ・クッタ法の変更部分）
（下記，赤線部分の変更）
Private DX As Double
Private Y0() As Double
Private XST As Double
Private Y(60, 2) As Double
Private N As Integer
Private M As Integer
Sub データ設定()
N = 60: DX = 0.1: XST = 0
M = 2: ReDim Y0(M)
Y0(1) = 0
Y0(2) = 4
End Sub
・
・(中略）
・
Function g(ID, X, Y) As
Double
If ID = 1 Then
g = Y(2)
Else
g = -Y(1)
End If
End Function
・
・(後略）
・
実行結果
誤差が少ないことを確認する。
（２）ルンゲ・クッタ法を使わない高次微分方程式
（オイラー法による例題）
2
m
d x
dt
  kx
2
k /m  ,
2
t  s
2
とおくと
d x
ds
2
 x
 v,
ds
と変形しておく。
次のような数式モデルを設定
x n 1  x n
( n  1) h  nh
v n 1  v n
とすることができるので，
dx
ここで，s = nhとなる有限の刻み幅hを導入し
，
dv
ds
 x
( n  1) h  nh
 v n  x n  1  x n  hv n
  x n  v n  1  v n   hx n
処理の手順
①
②
③
④
⑤
DT = 刻み幅
V0 = 初速度，X0 = 0，T = 0とする．
T時刻の速度，位置をV0，X0とする．
T = T + DTとする．
T ≦ 最終時刻の間，以下を繰り返す．
・X = X0 ＋ V0 × DT
・V = V0 － X0 × DT
・T時刻の速度，位置をV，Xとする．
・X0 = X，V0=V
・T = T + DT
単純な方法による式定義
（厳密解との比較を含む）
単純な方法による結果
（誤差が蓄積している）
この理由（エネルギー保存則）
1
mV
2

2
1
kx
2
＝一定，
2
V 
dx

dt
dx ds

ds dt
dx
dx
w 
ds
ds

k
v
m
したがって
2
1 
m v 
2 
v x
2
k 
1
1
k
1
k 2
2
  kx 2  mv 2   kx 2 
v x
m 
2
2
m 2
2
2


＝一定
＝一定
一方，数式モデルから

x n  1  v n  1   x n  hv n   v n  hx n   1  h
2
2
2
2
2
 x
2
n
 vn
2

刻み幅が0であれば一定となるが，今刻み幅は有限の値にしているので
保存則が成立していない。
k
m
式の改良
現在時刻における速度と位置の値は，お互いに現在時刻における値と前
時刻との平均値を使用して計算するものとする。
x n 1  x n 
1
2
v n 1  v n h
v n 1  v n  
1
2
 x n 1  x n h
現在時刻における値がお互いに入っているので，計算上はループしてし
まうので，お互いに代入して整理すると，
x n 1 
1
4h
2
( 4  h
2
)  xn  4vn h

v n 1 
1
4h
両式を二乗して加えると
2
( 4  h
2
)  vn  4 xn h

確認（エネルギー保存則が満足されている）
x n 1  v n 1 
2
1
2
2
2
2

2
1

(4  h )
2
2
2
2
2
2
2
( 4  h
2
)  vn  4 xn h

2
2
2
2
2
(4  h )
2

(4  h ) xn  vn
2
2
(4  h )
2
2
2
 x
2
n
2
 vn
2
2
2
2
2
2
2
2
2
(16 h  ( 4  h ) )  x n  (16 h  ( 4  h ) )  v n
2

)  xn  4vn h
(4  h )
2

2
2
( 4  h )  x n  16 v n h  8 ( 4  h )  x n v n h  ( 4  h )  v n  16 x n h  8 ( 4  h )  x n v n h
2

(4  h )
( 4  h
2
(4  h )  xn  (4  h )  vn
2

2
2
2
(4  h )
2
2
2
2
考慮した結果
（３）かえる飛び法
次のようにお互いに相手の中間値を使用する方法
x n 1  x n  v n 1 / 2  h
v n  1 / 2  v n 1 / 2   x n h
論理的には，中心差分であることに注意しよう。
かえるとび法（元の式は変わらないが誤差が少ない）
（４）テーラ展開による方法
まず，微分方程式から高次の微分を求め，
テーラ展開の式に代入する。
dy
［例］
 x  y
2
dx
n
y
(n)

d y
dx
y
(1 )
とすると
n
 x  y,
2
y
(2)
 2x  y
(1 )
,
y
(3)
 2 y
(2)
,
y
(4)
 y
(3)
これを次のようなテーラ展開の式に代入する。
y( x  x)  y( x)  x  y
(1 )
( x) 
x
2
2!
y
(2)
( x) 
x
3!
3
y
(3)
( x) 
x
4
4!
y
(4)
( x)  
式の変形
次のように近似できる。
yn  xn  yn , yn
(1 )
2
(2)
 2 xn  yn , yn
(1 )
(3)
 2  yn , yn
(2)
(4)
 yn
(3)
x n 1  x n   x
y n 1  y n   x  y n 
(1 )
x
2
2!
 yn 
(2)
x
3
 yn 
(3)
3!
x
4
4!
 yn
(4)
更に次のように変形できる。
y
(1 )
 x  y,
y
(3)
 2 y
2
(2)
y
(2)
 2x  y
(1 )
 2 x  x  y,
 2  2 x  x  y,
2
2
y
(4)
 y
(3)
 2  2x  x  y
2
式の整理
すなわち，y
(5)
 y
(4)
y( x  x)  y( x)  x  y
(1 )
 y
( x) 
(3)
x
2
2!
 y( x)  x  y
(1 )
(x) 
x
。同様にして y
y
(2)
( x)  y
(3)
(n)
 y
(3)
 x3 x 4 x5 
( x)



4!
5!
 3!

2
y
(2)
( x)
2!
 y
(3)
2
3
4
5
2



x
x
x
x
x 
(3)
( x ) 1   x 



   y ( x ) 1   x 

2
!
3
!
4
!
5
!
2
!




 y( x)  x  y
(1 )
(x) 
x
2
y
(2)
( x)
2!
 y
(3)
2
3
4
5
2



x
x
x
x
x 
(3)
( x ) 1   x 



   y ( x ) 1   x 

2
!
3
!
4
!
5
!
2
!




単純化
一方，定義から
e
x
 1  x 
x
2

x
2!
3!
3

x
4

4!
x
5

5!
2



x
(1 )
(2)
(3)
x
y( x  x)  y( x)  x  y ( x) 
 y ( x )  y ( x )e  1   x 

2!
2! 

x
2
最終的な数式モデル
yn  xn  yn , yn
(1 )
2
(2)
 2 xn  yn , yn
(1 )
(3)
 2  yn
(2)
x n 1  x n   x ,
y n 1  y n   x  y
(1 )
n

x
2
2!
y
(2)
n
 y
(3)
n
2
 x
x 
e  1   x 

2! 

テーラ展開による方法
テーラ展開による方法（結果）
［例題２］
2
d x
dt
x
(2)
x
  x,
(1 )
2
 x
x
(3)
(t  h )  x
  1,
(1 )
x
(4)
( t )  hx
(2)
 x
(5)
(t ) 
h
2
  x
x
(3)
(t )
(n)
0
［結果］