Transcript 数値計算法I - 九州工業大学

数値計算法I
第4回
高次連立方程式へのニュートン法の応用
ガウスの消去法
2006年度
九州工業大学工学部電気工学コース用講義資料
講師：趙孟佑
[email protected]
1
高次連立方程式へのニュートン法の拡張
• 複数の独立変数をもつ連立方程式もニュートン法
で解ける
f (x, y)  0

g(x, y)  0
(1)

2
高次連立方程式へのニュートン法の拡張
• f(x,y),g(x,y)を近似値xi,yiの近傍でテーラー展開す
ると、

f
f
f (xi   x , yi   y )  f (xi , yi )   x   y 


x
y

g(x   , y   )  g(x , y )  g   g  
i
x i
y
i i
x
y

x
y

(2)
3
高次連立方程式へのニュートン法の拡張
• 近似値xi,yiからx, y動けば、解に到達したとすると、
y
ここでf(x,y)=0,g(x,y)=0
y
(xi,yi) 
x
x
f (xi   x , yi   y )  0

g(xi   x , yi   y )  0
(3)
4
高次連立方程式へのニュートン法の拡張
• (2)式の2次以上の項を無視して、

f
f
0  f (xi , yi )   x   y


x
y

0  g(x , y )  g   g 
i i
x
y

x
y

これを満たすx, yは、


 x  


  
y


但し、
fgy  f y g
fx g y  f y gx
f x g  fgx
fx g y  f y gx
(5)
(4)

f
f
fx  , f y 


x
y (6)

g  g , g  g
x
y

x
5y

高次連立方程式へのニュートン法の拡張
• よって、i番目の近似値からi+1番目の近似値に進め
るには、
xi1  xi   x

yi1  yi   y
(7)
x, yは(5)式で与えられる。
収束判定は


 xi1  xi  1


yi1  yi  2
の2つが同時に成り立つ時。1と2は十分に小さな数
（自分で決める）
6
高次連立方程式(例）
f (x, y)  x 3  3xy  y 3  0

2
2
 g(x, y)  x  y 1  0
• f(x,y),g(x,y)をx,yで偏微分すると、


f x  3x 2  3y, f y  3x  3y 2


g x  2x, g y  2y

これより、(5)と(7)を使って計算することができる。

7
高次連立方程式(例）
回数
1
2
3
4
5
x
1.00000000
1.00000000
0.95185185
0.95022108
0.95021879
y
0.00000000
0.33333333
0.31111111
0.31158110
0.31158345
f(x,y)
1.00000000
0.03703704
0.00411609
0.00001010
0.00000000
g(x,y)
0.00000000
0.11111111
0.00281207
0.00000288
0.00000000
8
例題のプログラム
implicit real*8 (a-h,o-z)
c 初期値の入力
write(6,*)' input x0'
read(5,*)x0
write(6,*)' input y0'
read(5,*)y0
c 繰り返し最大数
itermax=100
c 最小誤差
errlimx=1.d-7
errlimy=1.d-7
iter=0
xp=x0
yp=y0
c 繰り返しの開始
10 continue
iter=iter+1
c
deltaxとdeltayを計算
deltax=-(f(xp,yp)*gy(xp,yp)-fy(xp,yp)*g(xp,yp))/
# (fx(xp,yp)*gy(xp,yp)-fy(xp,yp)*gx(xp,yp))
deltay=-(fx(xp,yp)*g(xp,yp)-f(xp,yp)*gx(xp,yp))/
# (fx(xp,yp)*gy(xp,yp)-fy(xp,yp)*gx(xp,yp))
c 新しいxとyを計算
xn=xp+deltax
yn=yp+deltay
c
結果の出力
write(6,100)iter,xp,yp,f(xp,yp),g(xp,yp)
c 収束判定
if(dabs(deltax).lt.errlimx.and.
# dabs(deltay).lt.errlimy)goto 30
if(iter.gt.itermax)goto 30
c xp、yp値の更新
xp=xn
yp=yn
goto 10
30 continue
100 format(i5,4f19.8)
stop
end
c f(x,y)の計算
function f(x,y)
implicit real*8 (a-h,o-z)
f=x*x*x-3*x*y+y*y*y
return
end
c fのx微分の計算
function fx(x,y)
implicit real*8 (a-h,o-z)
fx=3*x*x-3*y
return
end
c fのy微分の計算
function fy(x,y)
implicit real*8 (a-h,o-z)
fy=-3*x+3*y*y
return
end
c g(x,y)の計算
function g(x,y)
implicit real*8 (a-h,o-z)
g=x*x+y*y-1
return
end
c gのx微分の計算
function gx(x,y)
implicit real*8 (a-h,o-z)
gx=2*x
return
end
c gのy微分の計算
function gy(x,y)
implicit real*8 (a-h,o-z)
gy=2*y
return
end
9
連立一次方程式
• 連立1次方程式は常に行列の形式で表される。
a11x1  a12 x2  a13 x3  a14 x4  b1

a21x1  a22 x2  a23 x3  a24 x4  b2

a31x1  a32 x2  a33 x3  a34 x4  b3

a41x1  a42 x2  a43 x3  a44 x4  b4

a11

a21
a31

a41
a12
a13
a22
a23
a32
a33
a42
a43
a14 x1  b1 
   
a24 x2  b2 

a34 x3  b3 
   
a44 x4  b4 
連立一次方程式の計算=行列の計算
10
n次元連立1次方程式の解法
n元連立1次方程式AX=b
 a11
a12

a22
 a21


ai2
 ai1


an1,1 an1,2

an,2
 an,1
a1 j
a1,n1
a2 j
a2,n1
aii
ai,n1
an1, j
an1,n1
an, j
an,n1
a1n  x1   b1 
   
a2n  x2   b2 
   
   
ain  xi   bi 
   
   
an1,n xn1 bn1
an,n 
 xn 
 
 bn 


式はn個
未知数はx1,x2,--,xi,---,xn-1,xnのn個
11
前進消去
上三角行列への変換
 a11
a12

a
22
 0


ai2
 ai1


an1,1 an1,2

an,2
 an,1
a1 j
a
2j
a1,n1
a
2,n1
aii
ai,n1
an1, j
an1,n1
an, j
an,n1
 x   b 
1
1
   
 x2   b2 
   
   
ain  xi   bi 
   
   
an1,n xn1 bn1
an,n 
 xn 
 
 bn 


a1n
a
2n
ここで1行目を使って、2行目の第1列を消す。
a21
a21
a21
a
, a
, b2  b2  b1
22  a22  a12
2 j  a2 j  a1 j
a11
a11
a11
12
前進消去
上三角行列への変換
a11 a12

a
22
 0


a
i2
 0


 0 an1,2

an,2
 0
a1 j
a
2j
a1,n1
a
2,n1
a
ii
a
i,n1
an1, j
an, j
an1,n1
an,n1
 x   b 
1
1
   
 x2   b2 
   
   
a
in  xi   bi 
   
   
an1,n xn1 bn1
 
an,n 
 xn 
 
 bn 


a1n
a
2n
他の行も同様に、1行目を使って、第1列の要素を消す。
ai1
ai1
ai1
a
, a
, bi bi  b1
i2  ai2  a12
ij  aij  a1 j
a11
a11
a11
13
前進消去
上三角行列への変換
a11 a12

22
 0 a


0
 0


0
 0

0
 0
a1 j
a
2j
a1,n1
a
2,n1
a
ii

a
i,n1
 j
an1,
 j
an,

an1,n1

an,n1
 x   b 
1
1
   
 x2   b2 
   
   
  xi   bi 
a
in
   
   
 xn1 b
an1,n
 
n1
 
an,n
 xn 
 
 b


n 
a1n
a
2n
新たな2行目を使って、3行目以降の第2列の要素を消す。
a
a
a
32
i2
i2












a

a

a
,
a

a

a
,
b

b

b
33
33
23
ij
ij
2j
i
i
2
a
a
a
22
22
22
14
前進消去
上三角行列への変換
a11 a12

22
 0 a


0
 0


0
 0

0
 0
a1 j
a
2j
a1,n1
a
2,n1

a
ii

a
i,n1
0

an1,n1
0
0
 x   b 
1
1
   
 x2   b2 
   
   
  xi   bi
a
in
   
   
 xn1 b
an1,n
 
n1
 
an,n

 xn 
 
 b


n
a1n
a
2n
同様にして、対角要素から下を全て0にする。
aim
aim
aij  aij  amj
,bi  bi  bi
amm
amm
15
後退代入
a11 a12
a1 j
a1,n1
a1n  x1   b1 

   
a
a
a
22
2j
2,n1
2n  x2   b2 
 0 a

   

   


  xi   bi
0
a
a
a
ii
i,n1
in
 0

   

   

 xn1 b
0
0
an1,n1
an1,n
 
n1
 0

 
0
0
0
an,n

 xn 
 
 b

 0

n
一番下のn行目は、係数が一個しかないので、
bn
xn 

an,n
16
後退代入
xnが判ると、後は上に順繰りに計算していく。
ni
 j xi j
bi  ai,i
xi 
j1
aii
前進消去と後退代入で求める方法を
ガウスの消去法と呼ぶ
行列に要素が一杯入っている時（密行列）に有効
17
連立方程式の解法例
2x1  4x2  2x3  6

 x1  3x2  2x3  5
 2x  2x  x  5
2
3
 1
これを行列形式で書くと

2 4 2x1  6

   
1 3 2x2  5

2 2 1

x3 
 
5

18
連立方程式の解法例
2 4 2x1  6

   
1 3 2x2  5

2 2 1

x3 
 
5



2 4 2 x1  6 

    (2)-(1)*1/2
0 1 1 x2  2 

0 2 1

x3 
 
1
 (3)-(1)*2/2
2 4 2x1  6

   
0 1 1x2  2

0 0 1

x3 
 
3

(3)’-(2)’*-2/1
19
連立方程式の解法例
2 4 2x1  6

   
0 1 1x2  2

0 0 1

x3 
 
3

これより、

x3  3
2 1 x3
x2 
 1
1
6  2  x3  4  x2
x1 
2
2
が求まる。
20
ピボット(pivot)選び
2x1  4x2  2x3  6

 x1  2x2  3x3  9
4x  6x  2x  8
2
3
 1

行列に直して、
2 4 2x1  6

   
1 2 3x2  9

4 6 2

x3 
 
8

真値はx1=2,x2=-1,x3=3

21
ピボット(pivot)選び
2 4 2x1  6

   
1 2 3x2  9

4 6 2

x3 
 
8


2 4 2 x1  6 

   
0 0 2 x2  6 

0 2 2

x3 
 
4

(2)-(1)*1/2
(3)-(1)*4/2

22
ピボット(pivot)選び
2 4 2 x1  6 

   
0 0 2 x2  6 

0 2 2

x3 
 
4

第3行の-2を0にするには、
(3)’-(2)’*-2/0

とするが、そのままでは0で割っているので
発散する。
23
ピボット(pivot)選び

2 4 2 x1  6 

   
0 0 2 x2  6 

0 2 2

x3 
 
4

は 2x1  4x2  2x3  6

2x3  6

 2x  2x  4
2
3

2番目の式と3番目を入れ替えても答えは同じ。
2x1  4x2  2x3  6 2 4 2 x  6 
1







2x

2x

4

2
3
0 2 2x2  4
 
2x3  6


0 0 2 

x3 
 
6 

24
ピボット(pivot)選び
2 4 2 x1  6 

   
0 0 2 x2  6 

0 2 2

x3 
 
4


2 4 2 x1  6 

   
0 2 2x2  4

0 0 2 

x3 
 
6 

は
に置き換え可能
解が求まる。

25
Pivoting
a11 a12

22
 0 a


0
 0


0
 0

0
 0
a1 j
a
2j
a1,n1
a
2,n1
0
a
ii

a
i,n1
0
 j
an1,
 j
an,

an1,n1

an,n1
0
 x   b 
 1   1 
 x2   b2 
   
   
  xi   bi 
a
in
   
   
 xn1 b
an1,n
 
n1
 
an,n
 xn 
 
 b


n 
a1n
a
2n
j列の係数を0にするときに、a”ii,----,a”n-1,j,a”n,j
の内で最大の行を探し、それがi行目に来るよ
うに行を入れ替える。
26
スケーリング
1000 940 380 x1   340 

  

0.002 0.01 0.005x2  0.001

0.6
0.2 
 1

x3 
 
 0.8 

係数の桁が違い過ぎて、誤差の元

各式を各行の最大値で割ってから計算する。
 1 0.94 0.38x1  0.34

   
1
0.5 x2  0.1 
0.2

0.6
0.2 
 1

x3 
 
0.8 

27
消去法プログラム例（例題3.5)
implicit real*8 (a-h,o-z)
parameter(nmax=4)
dimension a(nmax,nmax),x(nmax),b(nmax)
dimension atmp(nmax)
c
a(1,1)=1
a(1,2)=2
a(1,3)=6
a(1,4)=1
a(2,1)=2
a(2,2)=4
a(2,3)=2
a(2,4)=2
a(3,1)=1
a(3,2)=5
a(3,3)=7
a(3,4)=1
a(4,1)=2
a(4,2)=4
a(4,3)=1
a(4,4)=3
b(1)=7
b(2)=4
b(3)=2
b(4)=2
c
n=4
c
c* scaling
c
do i=1,n
colmax=abs(a(i,1))
do j=2,n
if(abs(a(i,j)).gt.colmax)colmax=abs(a(i,j))
end do
do j=1,n
a(i,j)=a(i,j)/colmax
end do
b(i)=b(i)/colmax
end do
c
c* forward deleting
c
do j=1,n-1
c
c* pivotting
aamax=abs(a(j,j))
imax=j
do i=j+1,n
if(dabs(a(i,j)).gt.aamax)then
aamax=dabs(a(i,j))
imax=i
end if
end do
c
28
消去法プログラム例（例題3.5)
c* exchange
c
do jj=1,n
atmp(jj)=a(j,jj)
end do
btmp=b(j)
c
do jj=1,n
a(j,jj)=a(imax,jj)
end do
b(j)=b(imax)
c
do jj=1,n
a(imax,jj)=atmp(jj)
end do
b(imax)=btmp
c
do i=j+1,n
do jj=j+1,n
a(i,jj)=a(i,jj)-a(j,jj)*a(i,j)/a(j,j)
end do
b(i)=b(i)-b(j)*a(i,j)/a(j,j)
end do
end do
c
c* back substitution
c
x(n)=b(n)/a(n,n)
c
do i=n-1,1,-1
sum=0
do ii=i+1,n
sum=sum+a(i,ii)*x(ii)
end do
x(i)=(b(i)-sum)/a(i,i)
end do
c
do i=1,n
write(6,100)i,x(i)
end do
100 format(i3,f19.4)
c
stop
end
29