非線形方程式の数値解法〜ニュートン法（2）〜

非線形方程式の数値解法〜ニュートン法（２）〜
荻田　武史
目的
・2元連立非線形方程式の解を数値計算で求める。 ! f1 (x1, x2 ) = 0
"
# f2 (x1, x2 ) = 0
・本講義では、前回に引き続き、ニュートン法
（Newton's method）によってx1 , x2 の近似解を求
める方法を学習し、そのアルゴリズムを理解した
上で、プログラミングを行い、様々な方程式の問題
に取り組む。
※前回講義資料　ニュートン法を理解しておくこと。
2
ニュートン法（２元連立非線形方程式）
! f1 (x1, x2 ) = 0
•  連立非線形方程式 "
# f2 (x1, x2 ) = 0
これらについて、テイラー展開を1次の項で打ち切り
連立一次方程式で近似する。以下の連立一次方程式の解を（x(k+1), y(k+1)）とする。
(k ) (k )
(k ) (k )
#
∂f
(x
,
x
)
∂f
(x
(k ) (k )
(k )
(k )
1
1
2
1
1, x 2 )
f
(x
,
x
)
+
(x
−
x
)
+
(x
−
x
% 1 1 2
1
1
2
2 )=0
%
∂x1
∂x2
$
(k ) (k )
(k ) (k )
∂f
(x
,
x
)
∂f
(x
% f (x (k, )x (k ) ) + 2 1 2 (x − x (k ) ) + 2 1, x2 ) (x − x (k ) ) = 0
2
1
2
1
1
2
2
%
∂x
∂x
&
1
2
3
ニュートン法（２元連立非線形方程式）続き
x1,x2の変位をΔx1, Δx2と置くと、
(k ) (k )
(k ) (k )
#
∂f
(x
,
x
)
∂f
(x
(k ) (k )
1
1
2
1
1, x 2 )
Δx1 +
Δx2 = 0
% f1 (x1, x2 ) +
%
∂x1
∂x2
$
(k ) (k )
(k ) (k )
% f (x (k, )x (k ) ) + ∂f2 (x1, x2 ) Δx + ∂f2 (x1, x2 ) Δx = 0
2
1
2
1
2
%
∂x
∂x
&
1
2
これを行列表記にすると
"
(k ) (k )
∂f
(x
1
1, x 2 )
$
$
∂x1
$
(k ) (k )
∂f
(x
$
2
1, x 2 )
$
∂x1
#
%
∂f1 (x1(k, )x2(k ) ) '
'" Δx1
∂x2
$
'
∂f2 (x1(k, )x2(k ) ) '$# Δx2
'
∂x2
&
% " − f (x (k, )x (k ) )
1
1
2
'=$
' $ − f (x (k, )x (k ) )
2
1
2
& #
4
%
'
'
&
ニュートン法（２元連立非線形方程式）続き
変位について求めると
" Δx
1
$
$ Δx2
#
"
(k ) (k )
∂f
(x
1
1, x 2 )
$
%
$
∂x1
' = −$
(k ) (k )
'
∂f
(x
$
2
1, x 2 )
&
$
∂x1
#
−1
%
∂f1 (x , x ) '
'
∂x2
'
∂f2 (x1(k, )x2(k ) ) '
'
∂x2
&
(k ) (k )
1
2
"
(k ) (k )
f
(x
, x2 )
1
1
$
$ f (x (k, )x (k ) )
# 2 1 2
%
'
'
&
これをニュートン法の漸化式
x1(k+1) = x1(k ) + Δx1, x2(k+1) = x2(k ) + Δx2 に適用。
! (k+1)
# x1
# x (k+1)
" 2
$ ! (k )
& = # x1
& # x (k )
% " 2
!
(k ) (k )
∂f
(x
1
1, x 2 )
#
$
∂x1
&−#
& # ∂f (x (k, )x (k ) )
% # 2 1 2
#
∂x1
"
−1
$
∂f1 (x1(k, )x2(k ) ) &
&
∂x2
&
(k ) (k )
∂f2 (x1, x2 ) &
&
∂x2
%
!
(k ) (k )
f
(x
# 1 1, x2 )
# f (x (k, )x (k ) )
" 2 1 2
5
$
&
&
%
ヤコビ行列
以下のような行列をヤコビ行列という。
" ∂f (x, y)
$ 1
∂x
$
F '(x, y) = $
∂f2 (x, y)
$
$
∂x
#
∂f1 (x, y) %
'
∂y
'
∂f2 (x, y) '
'
'
∂y
&
2
"
$
f
(x,
y)
=
2x
−
y
+ log x = 0
（例） # 1
2
$
% f2 (x, y) = x − xy − x +1 = 0
! f (x, y)
1
F(x, y) = #
# f2 (x, y)
"
$
&
&
%
(x > 0)
"
%
1
$ 2+
−2y '
x
とおくとF '(x, y) = $
'
$ 2x − y −1 −x '
#
&
6
【実習】2元連立非線形方程式対応
•  関数の定義方法
（例） ! f (x , x )
1
1
2
F =#
# f2 (x1, x2 )
"
$ ! 2x − x 2 + log x
1
2
1
&=#
& # x 2 − x x − x +1
1 2
1
% " 1
$ !
$
&=# 0 &
& " 0 %
%
(x1 > 0)
>> F = inline('[2*x1 -‐ x2^2 + log(x1); x1^2 -‐ x1*x2 -‐ x1 + 1]') （注意） Fは2次元ベクトルで、F = [f1(x1,x2); f2(x1,x2)] >> Fd = inline('[2 + 1/x1, -‐2*x2; 2*x1 -‐ x2 – 1, -‐x1]') （注意） Fdは2x2のヤコビ行列 7
【実習】2次元連立非線形方程式対応
funcHon [x,k] = mynewton2(F,Fd,x0) tol = 1e-‐6; x = x0; % 初期値 for k = 1:50 fv = F(x(1),x(2)); Fm = Fd(x(1),x(2)); % 収束条件 if (norm(fv) <= tol) break; end d = -‐Fm\fv; x = x + d; end ← 連立一次方程式を解いて変位 d を計算
8
演習問題
•  作成したプログラムで以下の問題を解いてみ
よう（初期値x0は自分で適当に与える）。
(1)
(2)
(3)
" f (x, y) = sin(xy) − 0.51x − 0.32y = 0
#
$ g(x, y) = y cos(xy) − 0.51 = 0
" f (x, y) = −5x + 2sin x + cos y = 0
#
$ g(x, y) = 4 cos x + 2sin y − 5y = 0
"
$ f (x, y) = x(x − 3)2 − y 2 +1 = 0
#
x
−x
2
$
% g(x, y) = y(e − e ) − 4y +1 = 0
9

Download Report