Excelで定積分・微分方程式

８．３常微分方程式の数値解法
数値解法では，解析的に積分できないような問題
を扱うことが多い
オイラーってどんな人？
数学者・物理学者であり、天体物理学者。
スイスのバーゼル生まれ
現ロシアのサンクトペテルブルクで死んだ。
1707/4/15-1783/9/18
Leonhard Euler
８．３．１オイラー法
（１）前進差分によるオイラー法
微分方程式
dy
 g (t , y ),
dt
y (0)  C1
［例］前進差分近似
y (t  t )  y (t )
 g (t , y (t ))  y (t  t )  y (t )  g (t , y (t )) t
t
次のような繰返し計算で近似する
y
(t1 , y1 )
y (0)  C1
(t0 , y0 )
y (t )  y (0)  g (0, y (0))t
y (2t )  y (t )  g (t , y (t ))t

g (t0 , y0 )t
t
0
t
あるいは，
差分近似が微分の近似であることから
（本来は，t は限りなく０に近い値）
dy
y (t  t )  y (t ) y (t  t )  y (t )
 lim


t

0
dt
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
dt
Excelでオイラー法（１）
もちろん，この式は次のようにして厳密解を求めることができるが，誤差評価
のためにこの例を求める。
1
1
t2
dy  (t  1)dt   dy   (t  1)dt  log y   t  C
y
y
2
［Excelの定義
例］
tとして現時刻をとる定義
（例）
dy
 (t  1) ydt
dt
Excelでオイラー法（２）
グラフを描いてみると
tとして現時刻をとる定義のほうが誤差が少ない
（２）修正オイラー法（修正法１）
データ点における微分係数で増分を計算し，
増分から得られた座標点から更に増分を計算し，
両増分の平均値を最終的な増分とする。
k1  g (t , y )t
y
(t1 , y1 )
k 2  g (t  t , y  k1 )t
k2
1
y (t  t )  y (t )  (k1  k 2 )
2
(t0 , y0 )
0
k1
t
t
Excelによる修正法１
式定義
結果グラフ
（３）修正オイラー法（修正法２）
データ点における微分係数で増分を計算し，
増分の１／２から得られた座標点から増分を計算し，
この増分を最終的な増分とする。
k1  g (t , y )t
y
t
k1
k 2  g (t  , y  )t
2
2
y (t  t )  y (t )  k 2
(t1 
t
k
, y1  1 )
2
2
(t0 , y0 )
0
k1
k2
t
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
k2  g (t    t , y    k1 )t
y(t  t )  y(t )  (c1k2  c2 k2 )
［注］以下，式の展開による導出を示すが，
式の展開自体を覚える必要はない。
どのような流れで導出されているかを理解し，
最終的な結果を確認すること。
式の展開自体が必要なときは，資料を参照すれば良い。
確認
k2 をテーラ展開すると
k 2  g (t    t , y    k1 )t  g (t , y )    t
dg
dg
   k1
 O(t 2 , y 2 )
dt
dy
dg
dg
 g (t , y )    t
   g (t , y )t
 O(t 2 , y 2 )
dt
dy
第３式に代入すると
y (t  t )  y (t )  (c1k 2  c2 k 2 )


dg
dg
 y (t )  c1t  g (t , y )  c2 t  g (t , y )    t
   g (t , y )t
 O(t 2 )
dt
dy


dg
dg
 y (t )  (c1  c2 )t  g (t , y )    c2 t 2
   g (t , y )c2 t 2
 O(t 3 )
dt
dy
テーラ展開の式と比較すると
1
1
(c1  c2 )  1,   c2  ,   c2 
2
2
1
修正法１ c1  c2  ,     1
2
修正法２ c1  0, c2  1,    
1
2
８．３．２ルンゲ・クッタ法（Runge-Kutta method）
（１）２次のルンゲ・クッタ法(その１）
k1  g (t j , x j ), k2  g (t j  b1h, x j  b2hk1 ), x j 1  x j  h(a1k1  a2k2 ) とおくと
k 2  g (t j  b1h, x j  b2 hk1 )  g (t j , x j )  b1h
 g (t j , x j )  b1h
g (t j , x j )
t
 b2 hg(t j , x j )
g (t j , x j )
t
g (t j , x j )
 b2 hk1
g (t j , x j )
x
 O(h 2 )
 O(h 2 )
x


g (t j , x j )
g (t j , x j )

2 


 x j 1  x j  ha1 g (t j , x j )  a2  g (t j , x j )  b1h
 b2 hg(t j , x j )
 O(h ) 

t

x



g (t j , x j )
g (t j , x j ) 

3
 x j  h(a1  a2 ) g (t j , x j )  h a2b1
 a2b2 g (t j , x j )
  O(h )
t
x


2
（１）２次のルンゲ・クッタ法（その２）
d 2 x dg(t , x) g (t , x) g (t , x)
定義から 2 


g (t , x) であるからテーラ展開すると
dt
dt
t
x
h2
x j 1  x j  hg(t j , x j ) 
2
 g (t j , x j ) g (t j , x j )

3

g
(
t
,
x
)

O
(
h
)


j
j
x
 t

係数を比較して
a1  a2  1, a2b1 
1
1
, a2b2 
2
2
1
, b1  b2  1 を選択して
2
k1  g (t j , x j )
これらを満足する値として a1  a2 
k2  g (t j  h, x j  hk1 )
x j 1  x j 
h
k1  k2 
2
この式はホイン（Heun）の２次公式とも呼ばれる。
y j 1  y j  hg(t j , y j ) において
（２）３次のルンゲ・クッタ法
k1  g (t j , y j ), k2  g (t j  b1h, y j  b2 hk1 ), k3  g (t j  c1h, y j  c2 hk1  c3hk2 )
y j 1  y j  h(a1k1  a2 k2  a3k3 )
とおくと，２次の場合と同様，係数比較によって a1 , a2 , a3 , b1 , b2 , c1 , c2 , c3
を決めることができる。
k1  g (t j , x j )
h
h
k 2  g (t j  , x j  k1 )
2
2
k3  g (t j  h, x j  hk1  2hk2 )
x j 1  x j 
h
k1  4k2  k3 
6
この式はクッタ（Kutta）の公式と呼ばれる。
（３）４次のルンゲ・クッタ法
y j 1  y j  hg(t j , y j ) において
k1  g (t j , y j ), k 2  g (t j  b1h, y j  b2 hk1 ),
k3  g (t j  c1h, y j  c2 hk2 ), k 4  g (t j  d1h, y j  d 2 hk2 ),
y j 1  y j  h(a1k1  a2 k 2  a3k3  a4 k 4 )
とおくと，２次の場合と同様，係数比較によって a1 , a2 , a3 , b1 , b2 , c1 , c2 , d1 , d2
を決めることができる。
k1  g (t j , x j )
h
h
k 2  g (t j  , x j  k1 )
2
2
h
h
k3  g (t j  , x j  k 2 )
2
2
k 4  g (t j  h, x j  hk3 )
x j 1  x j 
h
k1  2k2  2k3  k4 
6
この式はルンゲ・クッタの式と呼ばれる。
（４）ルンゲ・クッタ・ギルの公式
（Runge-Kutta-Gill method）
４次のルンゲ・クッタ型の公式の一つとして，次の公式がある。
k1  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 
hk2 
2
2
2 
2



1
1  

k 4  g  t j  h, x j 
hk2  1 
hk3 
2
2 



h
1 
1 


x j 1  x j   k1  21 
k

2
1

k

k
 2 
 3 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 
 f1 (t , x1 ,  , xn ) 
x 
 f (t , x ,  , x ) 
1
n 
x   2 , f   2




 


x
f
(
t
,
x
,

,
x
)
1
n 
 n
 n
①増分ベクトル k1 を求める。
k1  f (t , x )
 k11   f1 (t , x1 ,  , xn ) 
k   f (t , x ,  , x ) 
1
n 
 12    2
   


  

k
f
(
t
,
x
,

,
x
)
1
n 
 1n   n
②増分 k2 ～ k4 を求める。
h
hk
k 2  f (t  , x  1 )
2
2
h
hk
k3  f (t  , x  2 )
2
2
k 4  f (t  h, x  hk3 )
③ x(t  h) を求める。
x (t  h)  x (t ) 
h
k1  2k2  2k3  k4 
6
（２）プログラム例
次の連立１次微分方程式を解くプログラムを例とする。
dy1
dy2
 x  y2 ,
x
dx
dx
y1 (0)  0,
y2 (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
８．３．４高階の常微分方程式
（１）考え方
次の微分方程式を考える。
d nx
d n 1 x
dx
a0 n  a1 n1    an 1  an x  0
dt
dt
dt
各次数の時間微分を別々の変数と考えると
dxn 1
dx1
dx2
 x2 ,
 x3 , ,
 xn ,
dt
dt
dt
dxn
1
  a1 xn  a2 xn 1    an x1 
dt
a0
x  x1 ,
という連立１次微分方程式の形式であるため，
ルンゲ・クッタ法を適用することができる。
［例題］
2
d x
m 2  kx
dt
解析的に解けるが，
厳密解が分かっているので，
誤差検証のために
あえてこの例題を取り上げる。
（連立ルンゲ・クッタ法の変更部分）
（下記，赤線部分の変更）
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
・
・(後略）
・
実行結果
誤差が少ないことを確認する。
（２）ルンゲ・クッタ法を使わない高次微分方程式
（オイラー法による例題）
d 2x
m 2  kx
dt
k / m  , t  s
ここで，s = nhとなる有限の刻み幅hを導入し
，
次のような数式モデルを設定
2
とおくと
d 2x
 x
2
ds
とすることができるので，
dx
 v,
ds
と変形しておく。
dv
 x
ds
xn 1  xn
 vn  xn 1  xn  hvn
(n  1)h  nh
vn 1  vn
  xn  vn 1  vn  hxn
(n  1)h  nh
処理の手順
①
②
③
④
⑤
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
1 2
2
mV  kx ＝一定，
2
2
V
dx dx ds dx
dx k
k

 w  
v
dt ds dt ds
ds m
m
したがって
2
1 
k  1 2 1 2 k 1 2 k 2

  kx  m v   kx  v  x 2
m v 
2 
m  2
2
m 2
2
v2  x2


＝一定
＝一定
一方，数式モデルから


xn21  vn21  xn  hvn   vn  hxn   1  h2 xn2  vn2
2
2

刻み幅が0であれば一定となるが，今刻み幅は有限の値にしているので
保存則が成立していない。
式の改良
現在時刻における速度と位置の値は，お互いに現在時刻における値と前
時刻との平均値を使用して計算するものとする。
1
vn 1  vn   xn 1  xn h
2
1
xn 1  xn  vn 1  vn h
2
現在時刻における値がお互いに入っているので，計算上はループしてし
まうので，お互いに代入して整理すると，

1
2
xn 1 
(
4

h
)  xn  4vn h
2
4h

vn 1 

1
2
(
4

h
)  vn  4 xn h
2
4h
両式を二乗して加えると

確認（エネルギー保存則が満足されている）
xn21  vn21 



2
1
1
2
2
(
4

h
)

x

4
v
h

(
4

h
)  vn  4 xn h
n
n
2 2
2 2
(4  h )
(4  h )

2
(4  h 2 ) 2  xn2  16vn2 h 2  8(4  h 2 )  xn vn h  (4  h 2 ) 2  vn2  16xn2 h 2  8(4  h 2 )  xn vn h

(4  h 2 ) 2
(16h 2  (4  h 2 ) 2 )  xn2  (16h 2  (4  h 2 ) 2 )  vn2 (4  h 2 ) 2  xn2  (4  h 2 ) 2  vn2


2 2
(4  h )
(4  h 2 ) 2


(4  h 2 ) 2 xn2  vn2
2
2


x

v
n
n
(4  h 2 ) 2
考慮した結果
（３）かえる飛び法
次のようにお互いに相手の中間値を使用する方法
xn1  xn  vn1/ 2  h
vn1/ 2  vn1/ 2   xn h
論理的には，中心差分であることに注意しよう。
かえるとび法（元の式は変わらないが誤差が少ない）
（４）テーラ展開による方法
まず，微分方程式から高次の微分を求め，
テーラ展開の式に代入する。
［例］
y (n)
dy
 x2  y
dx
d n y とすると
 n
dx
y (1)  x 2  y,
y ( 2)  2x  y (1) ,
y (3)  2  y ( 2) ,
y ( 4)  y (3)
これを次のようなテーラ展開の式に代入する。
2
3
4

x

x

x
y ( x  x)  y( x)  x  y (1) ( x) 
 y ( 2 ) ( x) 
 y ( 3) ( x ) 
 y ( 4 ) ( x)  
2!
3!
4!
式の変形
次のように近似できる。
yn(1)  xn2  yn , yn( 2)  2 xn  yn(1) , yn(3)  2  yn( 2) , yn( 4)  yn(3)
xn1  xn  x
2
3
4

x

x

x
yn1  yn  x  yn(1) 
 yn( 2) 
 yn(3) 
 yn( 4)
2!
3!
4!
更に次のように変形できる。
y (1)  x 2  y,
y ( 2)  2 x  y (1)  2 x  x 2  y,
y ( 3)  2  y ( 2 )  2  2 x  x 2  y ,
y ( 4 )  y ( 3)  2  2 x  x 2  y
式の整理
すなわち，y
( 5)
y
( 4)
y
( 3)
( n)
(3)
y

y
。同様にして
 x 3 x 4 x 5 
x 2 ( 2)
( 3)
y ( x  x)  y ( x)  x  y ( x) 
 y ( x)  y ( x)



2!
4!
5! 
 3!
(1)
2

x
 y ( x)  x  y (1) ( x) 
 y ( 2 ) ( x)
2!

x 2 x 3 x 4 x 5  (3) 
x 2 
( 3)
 y ( x)1  x 



  y ( x)1  x 

2!
3!
4!
5! 
2! 


2

x
 y ( x)  x  y (1) ( x) 
 y ( 2 ) ( x)
2!

x 2 x 3 x 4 x 5  (3) 
x 2 
( 3)
 y ( x)1  x 



  y ( x)1  x 

2!
3!
4!
5! 
2! 


単純化
一方，定義から
x 2 x 3 x 4 x 5
e  1  x 




2!
3!
4!
5!
x
2
2



x

x
(1)
( 2)
( 3)
x
y( x  x)  y( x)  x  y ( x) 
 y ( x)  y ( x)e  1  x 

2!
2
!


最終的な数式モデル
yn(1)  xn2  yn , yn( 2)  2 xn  yn(1) , yn(3)  2  yn( 2)
xn 1  xn  x,
x 2 ( 2)
x 2 
( 3)  x
yn 1  yn  x  y 
 yn  yn e  1  x 

2!
2! 

(1)
n
テーラ展開による方法
テーラ展開による方法（結果）
［例題２］
d 2x
 x
2
dt
x ( 2)   x, x (3)  1, x ( 4)  x (5)    x ( n )  0
h ( 3)
(1)
(1)
( 2)
 x (t  h)  x (t )  hx (t )  x (t )
2
［結果］

Download Report