● 講義資料

● 講義資料

▼ 講義予定

• 常微分方程式の初期値問題の数値解法 – Runge-Kutta Methodの安定性 – Symplectic解法

● 前回の講義のまとめ

★ Runge-Kutta法

• 次の計算スキームをRunge-Kutta 法と呼び,s をその段数と呼ぶ.

ki =f(x0+h

s

X

j=1

aijkj), i= 1, . . . , s,

x1=x0+h

s

X

i=1

biki,

• 特に,i≤j に対しては, aij= 0が成り立つものを陽的 Runge-Kutta法と呼ぶ. すなわち, ki =f(x0+h

i−1

X

j=1

aijkj), i= 1, . . . , s,

x1=x0+h

s

X

i=1

biki, である. 具体的に書けば,











k1=f(x0),

k2=f(x0+ha21k1), x1=x0+h(b1k1+b2k2)

s= 2,















k1=f(x0),

k2=f(x0+ha21k1),

k3=f(x0+ha31k1+ha32k2), x1=x0+h(b1k1+b2k2+b3k3)

s= 3,

などとなる.

• Runge-Kutta法の代表的なものは,以下の通りである.

改良オイラー法 （２段２次）

b1= 0, b2= 1, a21= 1 2,











k1=f(x0), k2=f(x0+h

2k1), x1=x0+hk2

ホインの２次公式 （２段２次）

b1=b2=1

2, a21= 1,











k1=f(x0), k2=f(x0+hk1), x1=x0+h

2(k1+k2) ホインの３次公式 （３段３次）

b1= 1/4, b2= 0, b3= 3/4, a21= 1/3, a31= 0, a32= 2/3























k1=f(x0), k2=f(x0+h

3k1), k3=f(x0+2h

3 k2), x1=x0+h

4(k1+ 3k3) クッタの３次公式 （３段３次）

b1= 1/6, b2= 2/3, b3= 1/6, a21= 1/2, a31=−1, a32= 2























k1=f(x0), k2=f(x0+h

2k1), k3=f(x0−hk1+ 2hk2), x1=x0+h

6(k1+ 4k2+k3) Runge-Kutta法 （４段４次公式）

a21= 1/2, a31= 0, a32= 1/2, a41= 0, a42= 0, a43= 1, b1= 1/6, b2= 1/3, b3= 1/3, b4= 1/6































k1=f(x0), k2=f(x0+h

2k1), k3=f(x0+h

2k2), k4=f(x0+h

2k3), x1=x0+h

6(k1+ 2k2+ 2k3+k4)

通常「Runge-Kutta法」と言うときには,上記の４段４次のRunge-Kutta法を指す.

• s段のRunge-Kutta型公式がp次となる条件を求めるためには, 以下の方法が考えられる.

（多くの書籍では,この方法を組み合わせ論的に整理した証明を用いている）

▼ テイラー展開を書く

x(t+h) =x(t) +hx^′(t) +h²

2 x^′′(t) +h³

3!x⁽³⁾(t) +h⁴

4!x⁽⁴⁾(t) +O(h⁵), ここで,x[0, T]−→Rⁿ とすれば,t∈[0, T]を固定したとき,

f(x(t)) =f ∈Rⁿ,

f^′(x(t)) =f^′∈Hom(Rⁿ,Rⁿ), であることに注意して,

x^′(t) =f, x^′′(t) = d

dtf(x(t)) =f^′(x(t))x^′(t) =f^′f, となる. よって,

x(t+h) =x(t) +hx^′(t) +h²

2 x^′′(t) +O(h³)

=x(t) +hf+h²

2 f^′f+O(h³) が成り立つ.

▼ ki の展開式を計算する ki をhの関数として展開すると,以下が成り立つ.

ki(h) =f(x0+hX

aijkj(h)), ki^′(h) =f^′(x0+hX

aijkj(h)) hX

aijkj^′(h) +X

aijkj(h) , これらのh= 0での値は,

ki(0) =f, k_i^′(0) =f^′X

aijkj

=X aij

f^′f となる. これを書き直せば,

ki(0) =f, k^′_i(0) =X

aij

f^′f

▼ テイラー展開に代入する これらのki(0),k^′_i(0)を使って, x1=x0+hX

biki(h)

=x0+hX

biki(0)

+h²X

bik^′_i(0)

+O(h³)

=x0+hX bi

f+h²X biaij

f^′f +O(h³) が成り立つ.

▼ テイラー展開に代入する この計算を比較すれば,

hの項 X

bi= 1, h² の項 X

biaij =1 2, をみたせば,局所離散化誤差をh³ とすることができる.

• より高次の条件式は以下の通りである.

３次の条件式







Xbiaijaik = 1 3, Xbiaijajℓ =1

6

４次の条件式























Xbiaijaikaiℓ=1 4, Xbiaijajkaiℓ= 1

8, Xbiaijajkajℓ= 1 12, Xbiaijajkakℓ = 1

24

• これらの条件式を各sについて具体的に書けば以下のようになる.

▼ ２段２次公式

b1+b2= 1, b2a21= 1

2

▼ ３段３次公式

b1+b2+b3= 1,

b2a21+b3(a31+a32) = 1 2, b2a²₂₁+b3(a31+a32)²= 1 3, b3a32a21=1

6

▼ ４段４次公式

b1+b2+b3+b4= 1,

a21b2+a31b3+a32b3+a41b4+a42b4+a43b4=1 2, a²₂₁b2+ (a31+a32)²b3+ (a41+a42+a43)²b4=1

3, a21a32b3+ (a21a42+a31a43+a32a43)b4= 1

6, a³₂₁b2+ (a31+a32)³b3+ (a41+a42+a43)³b4=1

4,

a32a21(a31+a32)b3+b4(a42a21+a43(a31+a32))(a41+a42+a43) = 1 8, a²₂₁a32b3+ (a²₂₁a42+ (a31+a32)²a43)b4= 1

12, a21a32a43b4= 1

24

これらの条件を満たせば,s段p次公式,すなわち,

|x(tn)−xn| ≤Ch^p となる公式となる.

• なお,１次の条件式P

bi= 1 が成り立つことは,計算スキームが意味を持つための必要条件 であるので,適切性の条件と呼ばれる. 逆に,Pbi = 1さえ成り立っていれば, aij の数値の 入力ミスがあっても, より低い次数の数値解をつくることができてしまうことに注意が必要 である.

• 陽的Runge-Kutta型公式については,以下の事実が知られている.

– s段でs+ 1次以上の公式をつくることはできない.

– s≥5の時,s段s次公式をつくることはできない.

• なお, Runge-Kutta法の係数{aij},{bi},{ci},ただし,ci=P^s

j=1aij を c1 a11 · · · a1s

... ... ... cs as1 · · · ass

b1 · · · bs

と表記することが多く,これをButcher 行列と呼ぶ.

• Butcher行列を使って,代表的な陽的ルンゲ・クッタ型公式を書くと以下のようになる.

1段 1次公式

前進オイラー法 後退オイラー法 0

1

1 1 1 2段 2次公式

改良オイラー法 ホインの２次公式 0

1/2 1/2

0 1

0

1 1

1/2 1/2 3段 3次公式

ホインの３次公式 クッタの３次公式 0

1/3 1/3

2/3 0 2/3

1/4 0 3/4

0 1/2 1/2

1 −1 2

1/6 2/3 1/6 4段 4次公式

ルンゲ・クッタの公式 ルンゲ・クッタ・ジルの公式 0

1/2 1/2

1/2 0 1/2

1 0 0 1

1/6 1/3 1/3 1/6

0

1/2 1/2 1/2 √¹

2−¹₂ 1−√¹ 2

1 0 −^√1₂ 1 + ^√1₂ 1/6 ¹₃(1−√¹

2) ¹₃(1 + √¹

2) 1/6

● 講義資料

▼ 単振動の数値計算例

以下の図は, 種々の方法で x^′′ = −x の初期値問題の数値解を計算した例である. x(0) = 1.0, x^′(0) = 0.0,h= 0.01t∈[0,16]である.

それぞれの左上図は位相の厳密解からの差であり,左下図は全エネルギーの初期値での値E0= 0.5 からの差である.

0 2 4 6 8 10 12 14 16 +0.0e+00 +1.0e-02 +2.0e-02 +3.0e-02 +4.0e-02 +5.0e-02 +6.0e-02 +7.0e-02 +8.0e-02 +9.0e-02 -6.0e-04

-5.0e-04 -4.0e-04 -3.0e-04 -2.0e-04 -1.0e-04 +0.0e+00

x’’ = -x (Euler)

-1.0 0.0 +1.0

-1.0 0.0 +1.0 x’’ = -x (Euler)

前進オイラー法

0 2 4 6 8 10 12 14 16 -8.0e-02 -7.0e-02 -6.0e-02 -5.0e-02 -4.0e-02 -3.0e-02 -2.0e-02 -1.0e-02 +0.0e+00 -6.0e-04

-5.0e-04 -4.0e-04 -3.0e-04 -2.0e-04 -1.0e-04 +0.0e+00

x’’ = -x (B_Euler)

-1.0 0.0 +1.0

-1.0 0.0 +1.0 x’’ = -x (B_Euler)

後退オイラー法

0 2 4 6 8 10 12 14 16 +0.0e+00 +5.0e-07 +1.0e-06 +1.5e-06 +2.0e-06 +2.5e-06 +0.0e+00

+5.0e-05 +1.0e-04 +1.5e-04 +2.0e-04 +2.5e-04 +3.0e-04

x’’ = -x (I_Euler)

-1.0 0.0 +1.0

-1.0 0.0 +1.0 x’’ = -x (I_Euler)

改良オイラー法

0 2 4 6 8 10 12 14 16 +0.0e+00 +5.0e-07 +1.0e-06 +1.5e-06 +2.0e-06 +2.5e-06 +0.0e+00

+5.0e-05 +1.0e-04 +1.5e-04 +2.0e-04 +2.5e-04 +3.0e-04

x’’ = -x (Heun_2)

-1.0 0.0 +1.0

-1.0 0.0 +1.0 x’’ = -x (Heun_2)

ホインの２次公式

0 2 4 6 8 10 12 14 16 -7.0e-07 -6.0e-07 -5.0e-07 -4.0e-07 -3.0e-07 -2.0e-07 -1.0e-07 +0.0e+00 +0.0e+00

+1.0e-09 +2.0e-09 +3.0e-09 +4.0e-09 +5.0e-09 +6.0e-09

x’’ = -x (Heun_3)

-1.0 0.0 +1.0

-1.0 0.0 +1.0 x’’ = -x (Heun_3)

ホインの３次公式

0 2 4 6 8 10 12 14 16 -1.2e-11 -1.0e-11 -8.0e-12 -6.0e-12 -4.0e-12 -2.0e-12 +0.0e+00 -1.4e-09

-1.2e-09 -1.0e-09 -8.0e-10 -6.0e-10 -4.0e-10 -2.0e-10 +0.0e+00

x’’ = -x (RK4)

-1.0 0.0 +1.0

-1.0 0.0 +1.0 x’’ = -x (RK4)

Runge-Kutta法

0 2 4 6 8 10 12 14 16 -3.0e-03 -2.0e-03 -1.0e-03 +0.0e+00 +1.0e-03 +2.0e-03 +3.0e-03 -5.0e-03

-4.0e-03 -3.0e-03 -2.0e-03 -1.0e-03 +0.0e+00 +1.0e-03

x’’ = -x (S_Euler)

-1.0 0.0 +1.0

-1.0 0.0 +1.0 x’’ = -x (S_Euler)

Symplectic Euler法

0 2 4 6 8 10 12 14 16 -1.4e-05 -1.2e-05 -1.0e-05 -8.0e-06 -6.0e-06 -4.0e-06 -2.0e-06 +0.0e+00 -1.0e-05

+0.0e+00 +1.0e-05 +2.0e-05 +3.0e-05 +4.0e-05 +5.0e-05 +6.0e-05 +7.0e-05

x’’ = -x (SV)

-1.0 0.0 +1.0

-1.0 0.0 +1.0 x’’ = -x (SV)

St¨olmer-Verlet法

0 2 4 6 8 10 12 14 16 -1.5e-08 -1.0e-08 -5.0e-09 +0.0e+00 +5.0e-09 +1.0e-08 +1.5e-08 +0.0e+00

+5.0e-09 +1.0e-08 +1.5e-08 +2.0e-08 +2.5e-08 +3.0e-08

x’’ = -x (Ruth)

-1.0 0.0 +1.0

-1.0 0.0 +1.0 x’’ = -x (Ruth)

Ruth 法

0 2 4 6 8 10 12 14 16 -1.6e-11 -1.4e-11 -1.2e-11 -1.0e-11 -8.0e-12 -6.0e-12 -4.0e-12 -2.0e-12 +0.0e+00 +2.0e-12 -8.0e-12

-6.0e-12 -4.0e-12 -2.0e-12 +0.0e+00 +2.0e-12 +4.0e-12 +6.0e-12 +8.0e-12 +1.0e-11 +1.2e-11 +1.4e-11

x’’ = -x (SS)

-1.0 0.0 +1.0

-1.0 0.0 +1.0 x’’ = -x (SS)

Sanz-Serna法

0 2 4 6 8 10 12 14 16 -4.0e-16 -2.0e-16 +0.0e+00 +2.0e-16 +4.0e-16 +6.0e-16 +8.0e-16 +1.0e-15 +1.2e-15 +1.4e-15 -1.4e-04

-1.2e-04 -1.0e-04 -8.0e-05 -6.0e-05 -4.0e-05 -2.0e-05 +0.0e+00

x’’ = -x (GL1)

-1.0 0.0 +1.0

-1.0 0.0 +1.0 x’’ = -x (GL1)

陰的中点法（１段２次Gauss-Legrandre法）

0 2 4 6 8 10 12 14 16 -5.0e-16 +0.0e+00 +5.0e-16 +1.0e-15 +1.5e-15 +2.0e-15 +2.5e-15 +3.0e-15 -2.5e-10

-2.0e-10 -1.5e-10 -1.0e-10 -5.0e-11 +0.0e+00

x’’ = -x (GL2)

-1.0 0.0 +1.0

-1.0 0.0 +1.0 x’’ = -x (GL2)

２段４次Gauss-Legrandre法

0 2 4 6 8 10 12 14 16 -2.5e-15 -2.0e-15 -1.5e-15 -1.0e-15 -5.0e-16 +0.0e+00 +5.0e-16 +1.0e-15 +1.5e-15 -3.5e-15

-3.0e-15 -2.5e-15 -2.0e-15 -1.5e-15 -1.0e-15 -5.0e-16 +0.0e+00 +5.0e-16 +1.0e-15

x’’ = -x (GL3)

-1.0 0.0 +1.0

-1.0 0.0 +1.0 x’’ = -x (GL3)

３段６次Gauss-Legrandre法

▼ 単振り子の数値計算例

以下の図は,種々の方法でx^′′=−sin(x)の初期値問題の数値解を計算した例である. x(0) = 2.0, x^′(0) = 0.0,h= 0.01t∈[0,16]である.

それぞれの左図は全エネルギーの初期値での値E0= 0.5からの差である.

+0.0e+00 +2.0e-02 +4.0e-02 +6.0e-02 +8.0e-02 +1.0e-01 +1.2e-01

0 2 4 6 8 10 12 14 16 Pendulum (Enery, Euler)

-3.0 -2.0 -1.0 +0.0 +1.0 +2.0 +3.0

-3.0 -2.0 -1.0 +0.0 +1.0 +2.0 +3.0 Pendulum (Euler)

前進オイラー法

-1.2e-01 -1.0e-01 -8.0e-02 -6.0e-02 -4.0e-02 -2.0e-02 +0.0e+00

0 2 4 6 8 10 12 14 16 Pendulum (Enery, B_Euler)

-3.0 -2.0 -1.0 +0.0 +1.0 +2.0 +3.0

-3.0 -2.0 -1.0 +0.0 +1.0 +2.0 +3.0 Pendulum (B_Euler)

後退オイラー法

-2.0e-06 +0.0e+00 +2.0e-06 +4.0e-06 +6.0e-06 +8.0e-06 +1.0e-05 +1.2e-05

0 2 4 6 8 10 12 14 16 Pendulum (Enery, I_Euler)

-3.0 -2.0 -1.0 +0.0 +1.0 +2.0 +3.0

-3.0 -2.0 -1.0 +0.0 +1.0 +2.0 +3.0 Pendulum (I_Euler)

改良オイラー法

-1.8e-05 -1.6e-05 -1.4e-05 -1.2e-05 -1.0e-05 -8.0e-06 -6.0e-06 -4.0e-06 -2.0e-06 +0.0e+00 +2.0e-06 +4.0e-06

0 2 4 6 8 10 12 14 16 Pendulum (Enery, Heun_2)

-3.0 -2.0 -1.0 +0.0 +1.0 +2.0 +3.0

-3.0 -2.0 -1.0 +0.0 +1.0 +2.0 +3.0 Pendulum (Heun_2)

ホインの２次公式

-6.0e-07 -5.0e-07 -4.0e-07 -3.0e-07 -2.0e-07 -1.0e-07 +0.0e+00 +1.0e-07

0 2 4 6 8 10 12 14 16 Pendulum (Enery, Heun_3)

-3.0 -2.0 -1.0 +0.0 +1.0 +2.0 +3.0

-3.0 -2.0 -1.0 +0.0 +1.0 +2.0 +3.0 Pendulum (Heun_3)

ホインの３次公式

-1.0e-10 -8.0e-11 -6.0e-11 -4.0e-11 -2.0e-11 +0.0e+00 +2.0e-11 +4.0e-11 +6.0e-11

0 2 4 6 8 10 12 14 16 Pendulum (Enery, RK4)

-3.0 -2.0 -1.0 +0.0 +1.0 +2.0 +3.0

-3.0 -2.0 -1.0 +0.0 +1.0 +2.0 +3.0 Pendulum (RK4)

Runge-Kutta法

-6.0e-03 -4.0e-03 -2.0e-03 +0.0e+00 +2.0e-03 +4.0e-03 +6.0e-03

0 2 4 6 8 10 12 14 16 Pendulum (Enery, S_Euler)

-3.0 -2.0 -1.0 +0.0 +1.0 +2.0 +3.0

-3.0 -2.0 -1.0 +0.0 +1.0 +2.0 +3.0 Pendulum (S_Euler)

Symplectic Euler法

-3.0e-05 -2.5e-05 -2.0e-05 -1.5e-05 -1.0e-05 -5.0e-06 +0.0e+00 +5.0e-06

0 2 4 6 8 10 12 14 16 Pendulum (Enery, SV)

-3.0 -2.0 -1.0 +0.0 +1.0 +2.0 +3.0

-3.0 -2.0 -1.0 +0.0 +1.0 +2.0 +3.0 Pendulum (SV)

St¨olmer-Verlet法

-3.0e-08 -2.0e-08 -1.0e-08 +0.0e+00 +1.0e-08 +2.0e-08 +3.0e-08

0 2 4 6 8 10 12 14 16 Pendulum (Enery, Ruth)

-3.0 -2.0 -1.0 +0.0 +1.0 +2.0 +3.0

-3.0 -2.0 -1.0 +0.0 +1.0 +2.0 +3.0 Pendulum (Ruth)

Ruth 法

-3.0e-11 -2.5e-11 -2.0e-11 -1.5e-11 -1.0e-11 -5.0e-12 +0.0e+00 +5.0e-12 +1.0e-11 +1.5e-11 +2.0e-11

0 2 4 6 8 10 12 14 16 Pendulum (Enery, SS)

-3.0 -2.0 -1.0 +0.0 +1.0 +2.0 +3.0

-3.0 -2.0 -1.0 +0.0 +1.0 +2.0 +3.0 Pendulum (SS)

Sanz-Serna法

+0.0e+00 +1.0e-06 +2.0e-06 +3.0e-06 +4.0e-06 +5.0e-06 +6.0e-06 +7.0e-06 +8.0e-06 +9.0e-06

0 2 4 6 8 10 12 14 16 Pendulum (Enery, GL1)

-3.0 -2.0 -1.0 +0.0 +1.0 +2.0 +3.0

-3.0 -2.0 -1.0 +0.0 +1.0 +2.0 +3.0 Pendulum (GL1)

陰的中点法（１段２次Gauss-Legrandre法）

-5.0e-12 +0.0e+00 +5.0e-12 +1.0e-11 +1.5e-11 +2.0e-11 +2.5e-11

0 2 4 6 8 10 12 14 16 Pendulum (Enery, GL2)

-3.0 -2.0 -1.0 +0.0 +1.0 +2.0 +3.0

-3.0 -2.0 -1.0 +0.0 +1.0 +2.0 +3.0 Pendulum (GL2)

２段４次Gauss-Legrandre法

-5.0e-15 -4.0e-15 -3.0e-15 -2.0e-15 -1.0e-15 +0.0e+00 +1.0e-15 +2.0e-15 +3.0e-15

0 2 4 6 8 10 12 14 16 Pendulum (Enery, GL3)

-3.0 -2.0 -1.0 +0.0 +1.0 +2.0 +3.0

-3.0 -2.0 -1.0 +0.0 +1.0 +2.0 +3.0 Pendulum (GL3)

３段６次Gauss-Legrandre法

▼ 相流による写像の計算例

以下の図は,x^′′=−xおよびx^′′=−sin(x)の相流で,相平面上の図形がどのように変化したか を示している図である. h= 0.05,t= 5,t= 10,t= 15での図を示している.

それぞれ,青破線の図がt= 0での図形であり,厳密解による相流の場合には,図形の面積が保存 されているはずである.

-1.0 0.0 +1.0

-1.0 0.0 +1.0 x’’ = -x (Euler)

-1.0 0.0 +1.0

-1.0 0.0 +1.0 Pendulum (Euler)

前進オイラー法

-1.0 0.0 +1.0

-1.0 0.0 +1.0 x’’ = -x (B_Euler)

-1.0 0.0 +1.0

-1.0 0.0 +1.0 Pendulum (B_Euler)

後退オイラー法

-1.0 0.0 +1.0

-1.0 0.0 +1.0 x’’ = -x (I_Euler)

-1.0 0.0 +1.0

-1.0 0.0 +1.0 Pendulum (I_Euler)

改良オイラー法

-1.0 0.0 +1.0

-1.0 0.0 +1.0 x’’ = -x (Heun_2)

-1.0 0.0 +1.0

-1.0 0.0 +1.0 Pendulum (Heun_2)

ホインの２次公式

-1.0 0.0 +1.0

-1.0 0.0 +1.0 x’’ = -x (Heun_3)

-1.0 0.0 +1.0

-1.0 0.0 +1.0 Pendulum (Heun_3)

ホインの３次公式

-1.0 0.0 +1.0

-1.0 0.0 +1.0 x’’ = -x (RK4)

-1.0 0.0 +1.0

-1.0 0.0 +1.0 Pendulum (RK4)

Runge-Kutta法

-1.0 0.0 +1.0

-1.0 0.0 +1.0 x’’ = -x (S_Euler)

-1.0 0.0 +1.0

-1.0 0.0 +1.0 Pendulum (S_Euler)

Symplectic Euler法

-1.0 0.0 +1.0

-1.0 0.0 +1.0 x’’ = -x (SV)

-1.0 0.0 +1.0

-1.0 0.0 +1.0 Pendulum (SV)

St¨olmer-Verlet法

-1.0 0.0 +1.0

-1.0 0.0 +1.0 x’’ = -x (Ruth)

-1.0 0.0 +1.0

-1.0 0.0 +1.0 Pendulum (Ruth)

Ruth 法

-1.0 0.0 +1.0

-1.0 0.0 +1.0 x’’ = -x (SS)

-1.0 0.0 +1.0

-1.0 0.0 +1.0 Pendulum (SS)

Sanz-Serna法

-1.0 0.0 +1.0

-1.0 0.0 +1.0 x’’ = -x (GL1)

-1.0 0.0 +1.0

-1.0 0.0 +1.0 Pendulum (GL1)

陰的中点法（１段２次Gauss-Legrandre法）

-1.0 0.0 +1.0

-1.0 0.0 +1.0 x’’ = -x (GL2)

-1.0 0.0 +1.0

-1.0 0.0 +1.0 Pendulum (GL2)

２段４次Gauss-Legrandre法

-1.0 0.0 +1.0

-1.0 0.0 +1.0 x’’ = -x (GL3)

-1.0 0.0 +1.0

-1.0 0.0 +1.0 Pendulum (GL3)

３段６次Gauss-Legrandre法

● 実習内容

（以下で「★」の数は推奨の程度を示します. 多いほど推奨の度合いが大きい. また,「†」は 難易度またはマニアックな程度を表します. 多いほど難しいかマニアックかです.）

1. （★★★）６回目資料の「単振動」と「単振り子」の常微分方程式の初期値問題の数値解を, Symplectic Euler法で計算しなさい.

2. （★★††）６回目資料の「単振動」と「単振り子」の常微分方程式の初期値問題の数値解 を,陰的中点法で計算しなさい.