PROGRAM EXERCISE26 IMPLICIT NONE ! Broyden法でf(x,y)=0, g(x,y)=0の解を求めるプログラム REAL*8 :: x_bro(1:2), f_bro(1:2), f_bro2(1:2), deltax_bro(1:2) INTEGER :: iter, i,j REAL*8 :: eps=1.0d-14 REAL*8 :: Binv(1:2,1:2), denom, nom1(1:2), nom2(1:2), deltaf_bro(1:2) ! 反復の回数カウンタ iter = 1 ! xの初期値 x_bro(1) = 0.0d0 x_bro(2) = 0.0d0 f_bro = f(x_bro(1:2)) PRINT *, "initial x =", x_bro PRINT *, "initial f(x)=", f_bro ! Jacobianの逆行列の初期値(単位行列に設定) Binv(1:2,1:2) = 0.0d0 Binv(1,1) = 1.0d0 Binv(2,2) = 1.0d0 ! 反復のループ DO ! Newton法で次の点を計算 ! \Delta x^{(k)} deltax_bro(1:2) = -MATMUL( Binv(1:2,1:2), f_bro(1:2)) ! x^{(k+1)} x_bro(1:2) = x_bro(1:2) + deltax_bro(1:2) ! f(x^{(k+1)}) f_bro2 = f(x_bro(1:2)) ! \Delta f deltaf_bro(1:2) = f_bro2(1:2) - f_bro(1:2) ! f_broにf(x^{(k+1)})を代入して次の反復計算で用いる f_bro = f_bro2 ! B^{-1}行列の分母部分 denom = DOT_PRODUCT( deltax_bro(1:2), MATMUL( Binv(1:2,1:2), deltaf_bro(1:2))) ! B^{-1}行列の分子部分(左側) nom1(1:2) = deltax_bro(1:2) - MATMUL( Binv(1:2,1:2), deltaf_bro(1:2)) ! B^{-1}行列の分子部分(右側) nom2(1:2) = MATMUL(deltax_bro(1:2), Binv(1:2,1:2)) ! B^{(k+1)}^{-1}を計算 DO i = 1, 2 DO j = 1, 2 Binv(i,j) = Binv(i,j) + nom1(i)*nom2(j)/denom END DO END DO ! 反復の途中経過を出力 PRINT *, "iter = ", iter PRINT *, "x = ", x_bro(1), x_bro(2) PRINT *, "f(x) = ", f_bro(1), f_bro(2) ! f_broが十分ゼロに近づいたらDOループから抜ける IF( ABS(SQRT(f_bro(1)**2 + f_bro(2)**2)) .LT. eps) EXIT ! 反復カウンタを1増やす iter = iter + 1 ! もし問題が起きた場合1000回で反復を強制的に終了する IF( iter .EQ. 1000) STOP "iteration did not converge." END DO CONTAINS ! f(x,y)の値を返す関数副プログラム FUNCTION f(x) REAL*8, INTENT(IN) :: x(1:2) REAL*8 :: f(1:2) f(1) = x(1)**2 + x(2)**2 - 10.0d0 f(2) = x(1) - x(2)**3 RETURN END FUNCTION f END PROGRAM EXERCISE26