ニュートン法による非線形連立方程式の解法 [ネコ騙し数学]
ニュートン法による非線形連立方程式の解法
f(x,y)、g(x,y)を2回微分可能な関数とする。
の近似解法について考えることにする。
(x,y)を連立方程式の解、(x⁰,y⁰)を(1)の解の推測値とする。
f(x,y)、g(x,y)は仮定より微分可能なので、次のようにテーラー展開することが可能である。
したがって、(x⁰,y⁰)が(x,y)の近傍の点とすると、
と近似することができる。
(x,y)は連立方程式(1)の解なので、f(x,y)=0、g(x,y)=0であるから、(4)、(5)式より、(1)式は
と書き換えることが可能で、この連立方程式を解くことによって、(1)の解の推測値(x⁰,y⁰)をもとに解(x,y)の推測値を求めることができる。
なおここで、
である。
(6)式を行列を用いて書き換えると、
ここで、
とおくと、(8)の行列式|J|≠0のとき、(7)式より
と解くことができる。
(9)式で得られた連立方程式の解(x,y)をあらたに(x⁰,y⁰)にし、(8)、(9)を再計算し、(x,y)を求める。
この操作を繰り返すほど、(1)の近づくことが予想される。
この方法は、そのまま、より一般の
に拡張することが可能で、この場合は次のようになる。
(12)から(非線形)連立方程式の近似解を求めることができるが、(12)を利用する場合、ヤコビ行列Jの逆行列J⁻¹を求める必要がある。しかし、連立方程式(11)を解くよりヤコビ行列Jの逆行列を求めることの方が大変なので――連立方程式を直接解く場合と比較すると、Jの逆行列を求める計算は計算量が約n倍!!――なので、連立方程式
をガウス消去法などで解き、
と計算するとよい。
このように(非線形)連立方程式の近似解を求める方法をニュートン(・ラプソン)法という。
参考までに、非線形連立方程式
をニュートン法で解くプログラムを示す。
! ニュートン・ラプソン法で非線形連立方程式f(x,y)=0,g(x,y)=0を解く
parameter (n=2)
real a(n,n+1)
limit = 10 ! 繰り返し計算の上限
eps=1.e-6 ! 収束判定に使用する何か
a=0. ! 行列の初期化
write(*,*) 'input x0,y0'
read(*,*) x,y ! 計算開始時の初期値設定
write(*,*)
do iter=1,limit
call partial(f,fx,fy,g,gx,gy,x,y)
! 連立方程式の係数の設定
a(1,1)=fx; a(1,2)=fy; a(1,3)=-f
a(2,1)=gx; a(2,2)=gy; a(2,3)=-g
call gauss(a,n,iflag) ! ガウス消去法で連立方程式を解く
if (iflag.eq.0) then ! エラー判定
write(*,*) 'エラー!! 連立方程式を解けない'
stop
end if
dx=a(1,3); dy=a(2,3) ! 修正量 Δx,Δyをセット
x=x+dx; y=y+dy ! 修正
err=amax1(abs(dx),abs(dy))
write(*,100) 'iter=',iter, 'x=',x,'y=',y
if (err.lt.eps) exit ! 収束したかの判定
end do
write(*,*)
if (iter.le.limit) then
write(*,*) '***solution***'
write(*,*) 'x=',x,'y=',y
else
write(*,*) '収束しない'
end if
100 format(a,i3,2x,a,f10.6,2x,a,f10.6)
end
subroutine partial(f,fx,fy,g,gx,gy,x,y) ! 偏微分などの計算
f=x*x+y*y-4 ! f(x,y)=x²+y²−4=0
g=x+2*y-3 ! g(x,y)=x+2y-3=0
fx=2.*x !∂f/∂x
fy=2.*y !∂f/∂y
gx=1. !∂g/∂x
gy=2. !∂g/∂y
end
subroutine Gauss(a,n,iflag) ! ガウス消去法
real a(n,n+1)
eps = 1.e-10
iflag = 1
! 前進消去
do k=1, n-1
! ピボット選択
i_max = k
do i=k+1, n
if (abs(a(i,k)).gt.abs(a(i_max,k))) i_max=i
end do
if (i_max.ne.k) then
do j=k,n
w=a(k, j)
a(k, j) = a(i_max,j)
a(i_max,j) = w
end do
end if
if (abs(a(k,k)).lt.eps) then
iflag = 0
return
end if
! ピボット選択終わり
do i=k+1, n
t=a(i,k)/a(k,k)
do j=k+1, n+1
a(i,j)=a(i,j)-t*a(k,j)
end do
end do
end do
if (abs(a(n,n)).lt.eps) iflag = 0
! 後退代入
do i=n,1, -1
d = a(i,n+1)
do j= n, i+1, -1
d=d- a(i,j)*a(j,n+1)
end do
a(i,n+1)=d/a(i,i)
end do
end