n

数値計算の練習
ニュートン・ラフソン法，
ロトカ・ボルテラの方程式
倉橋貴彦
(1) ニュートン・ラフソン法による
非線形連立方程式の計算
x  y 5  0
2
2
2
x
 y2 1  0
9
方程式をf,gとする。
f  x, y   x 2  y 2  5
x2
g  x, y  
 y 2 1
9
各方程式のテイラー展開を
正解
x   4.5  2.121
y
2
 0.707
2
一階微分の項まで示すと・・・
n
n
 n 1
f
f
n
x 
y
f  f 
x
y


n
n

g

g
 g n 1  g n 
x 
y

x
y

行列で示すと・・・
 f n
n 1
n
f 
 f   x
    
n
g
g


g
 
 
 x

n
f 

y  x 

n 
g  y 
y 
x、yが解の場合はf,gは零なので・・・
 f n
n 1
n
0
f
 
   x
    
n
0
 g   g
 x

上式を整理すると・・・
 f n

 x n
 g
 x

n
f 

y  x 

n 
g  y 
y 
n
f 
 x   f n
y 
 
n 
g  y    g 
y 
左辺の行列：
ここに変化量Δx、Δyは・・・
x   x n 1  x n 
    n 1
n

y
y

y
  

ヤコビ行列
よって変化量Δx、Δyは・・・
 f n
x   x
 
n

y


g
 
 x

1
f 
  f n
y 
 
n
g    g 
y 
n
よって、x、yの値の更新式は・・・
x
 
 y
n 1
 x  x 
   
 y  y 
n
具体例で考えると・・・
f  x, y   x 2  y 2  5
x2
g  x, y  
 y 2 1
9
ヤコビ行列は・・・
 f n

 x n
 g
 x

n
f 
  2xn
y  
n  2 n
g   x
9

y 
2 yn 

n
2y 

ヤコビ行列の逆行列は・・・
 f

 x n
 g
 x

n
1
f 

1
y 

n
4 n n
g 
n n
4
x
y

x y

9
y 
n
 2 yn
 2 n
 x
 9
 2 yn 

n
2x 

よって、x、yの更新式は・・・
x
 
 y
n 1
 x  x 
   
 y  y 
n
 f
 x   x
  
n
y


g
 
 x

n
n
1
f 
  f  n
y 
 
n
g    g 
y 
n
n

2
y
x
1
 2 n
  
 y  4 x n y n  4 x n y n  9 x
9
n


n 2
n 2


(
x
)

(
y
) 5 

 2y


n 2



(
x
)

 ( y n ) 2  1
2 x n   
  9

n
例えば、初期値を以下の様に設定して反復計算を行う。
収束判定定数ε=10-6とし、
|dx|<εおよび|dy|<εを満足したら計算を終了する。
x
 
 y
(0)
1.5 
 
0.5
表：反復計算の履歴
iteration
x
y
dx
dy
0
1.5000000
0.5000000
0.7500000
0.2500000
1
2.2500000
0.7500000
-0.1250000
-0.0416667
2
2.1250000
0.7083333
-0.0036765
-0.0012255
3
2.1213235
0.7071078
-0.0000032
-0.0000011
4
2.1213203
0.7071068
正解
x   4.5  2.121
y
2
 0.707
2
Fortranによるプログラム例 (newton.f)
implicit double precision ( a-h, o-z )
c
open(10,file='input.dat')
open(11,file='output.dat')
c
read(20,*) x0
read(20,*) y0
c
eps = 0.000001d0
iter = 0
c
222 continue
c
f0 = x0*x0 + y0*y0 - 5.d0
g0 = (x0*x0/9.d0) + y0*y0 -1.d0
c
x0 = x1
y0 = y1
iter = iter + 1
go to 222
c
111 continue
c
write(11,100) iter, x0, y0, dx, dy
c
100 format(i5,4f12.7)
c
stop
end
c
acoef = 1.d0 / ( 4.d0*x0*y0 - (4.d0/9.d0)*x0*y0 )
c
dx = acoef * ??????????
dy = acoef * ??????????
c
if ( dabs(dx) .lt. eps .and. dabs(dy) .lt. eps ) then
go to 111
else
end if
c
write(11,100) iter, x0, y0, dx, dy
c
x1 = x0 + dx
y1 = y0 + dy
Input.dat
1.50
0.50
初期値の設定による解の収束する値の違い
Case2
(4,3)
Case3
(-3,4)
Case1(1.5,0.5)
注：Case1では，
解の近傍を少し
過ぎてから戻って
収束をしている．
全て解である。初期値の設定に x   4.5  2.121
よっては，解の収束位置が異な
2
y


 0.707
ることに注意。
2
Case5
(-5,-5)
Case4
(3,-2)
ロトカ・ボルテラの方程式（捕食者・非捕食者の関係を表した連立微分方程式）
dx
 x(  y )  f x, y 
dt
dy
  y (  x)  g x, y 
dt
X:捕食者
Y：非捕食者

k x1  tf x n , y n


 txn   y n
修正オイラー法（２次のルンゲ・クッタ法）
による時間進展式
1
k x1  k x 2 
2
1
y n 1  y n  k y1  k y 2
2
x n 1  x n 




k x 2  tf x n 1, y n 1


 txn 1   y n 1






 t x n  tf x n , y n    y n  tg x n , y n

k y1  tg x n , y n


  ty n (  x n )

k y 2  tg x n 1, y n 1


  ty n 1   x n 1



 

  t y n  tg x n , y n    x n  tf x n , y n

Start
フローチャート例
n=0，Δt ，x0=x|t=0 ，y0=y|t=0 ，T
(※T/Δtが整数になるようにΔtを選択する。）
imax = T/Δt
T:全時間
n：ステップ数
imax：最大ステップ数
n*Δt, x0, y0

k x1  tf x n , y n


 txn   y n
n=n+1


n 1
n 1
k x 2  tf x n 1, y n 1
  y 
 t x  tf x , y    y
 tx
kx1, kx2, ky1, ky2の計算
n
1
k x1  k x 2 
2
1
y n 1  y n  k y1  k y 2
2
n

x n 1  x n 


k y1  tg x n , y n



  x 
  t y  tg x , y    x
yn+1
  ty
n 1
n 1
n
n
Yes
No

 tg x n , y n
  ty n (  x n )

n*Δt,
n

k y 2  tg x n 1, y n 1
xn+1,
n
n=imax
end
n
n

 tf x n , y n

Fortranによるプログラム例 (lotoca.f)
implicit double precision (a-h,o-z)
c
open(10,file='input.dat')
open(11,file='output.dat')
c
c----- computational conditions
c
read(10,*) x0
read(10,*) y0
read(10,*) delt
read(10,*) imax
c
read(10,*) alpha
read(10,*) beta
read(10,*) gamma
read(10,*) delta
c
c----- output of initial input value
c
write(11,*) 0.d0, x0, y0
c
do istep = 1,imax
c
c----- value of right hand side term at n step
c
call fx
& ( f0, alpha, beta, x0, y0 )
c
call fy
& ( g0, gamma, delta, x0, y0 )
c
c----- value update by Euler method
c
ax1 = x0 + delt * f0
ay1 = y0 + delt * g0
c
c----- value of right hand side term at n+1 step
c
call fx
& ( f1, alpha, beta, ax1, ay1 )
c
call fy
& ( g1, gamma, delta, ax1, ay1 )
c
c----c
akx1 = delt * f0
akx2 = delt * f1
c
aky1 = delt * g0
aky2 = delt * g1
c
c----- value update by modified Euler method
c
x1 = x0 + 0.5d0 * ( akx1 + akx2 )
y1 = y0 + 0.5d0 * ( aky1 + aky2 )
c
write(11,*) delt*istep, x1, y1
c
c----- step change
c
x0 = x1
y0 = y1
c
end do
c
stop
end
c
c==================================
c
subroutine fx
& ( f0, alpha, beta, x0, y0 )
c
c==================================
c
implicit double precision (a-h,o-z)
c
f0 = x0 * ( alpha - beta * y0 )
c
return
end
c
c==================================
c
subroutine fy
& ( g0, gamma, delta, x0, y0 )
c
c==================================
c
implicit double precision (a-h,o-z)
c
g0 = - y0 * ( gamma - delta * x0 )
c
return
end
計算結果の例
input.dat
計算条件
0.3
0.7
0.2
500
0.2
0.4
0.3
0.3
←x0 = 0.3
←y0 = 0.7
←Δt= 0.2
←time step = 500
←α = 0.2
←β = 0.4
←γ = 0.3
←δ = 0.3
x,yの値の時系列変化
レポート検討課題
(1) ニュートンラフソン法による非線形連立方程式の解を得るプログラムを
作成し，初期値を変えることで，収束する解の値の違いについて検討
をし，結果から分かることについて，式や図，説明を交えて整理せよ。
（※初期値を0,0とすると計算ができない理由についても合わせて検討
し，説明せよ。）
(2) ロトカ・ボルテラの方程式を修正オイラー法により解くプログラムを作成
し，係数α，β，γ，δの値および，初期値x0,y0，また時間刻みΔtの値を変
えることで，得られうる解の違いについて，結果から分かることについ
て式や図，説明を交えて整理せよ。

Download Report