Newton-Raphson 法の減速法のコードの書き方(FORTRAN)

Question

前提

fortranで次の非線形方程式 f(x)=tanh(x)+0.2x+0.3= 0 の数値解をNewton-Raphson 法を用いて計算しました。この時、初期値をx=1.3に設定すると解が収束しないので、減速法を用いて解を求めようとしましたが、計算がループしてしまい上手くいきません。

実現したいこと

非線形方程式 f(x)=tanh(x)+0.2x+0.3= 0 の数値解を、初期値をx=1.3に設定し、減速法を用いて求める。

発生している問題・エラーメッセージ

計算が終わらない（ループしている)。

該当のソースコード

PROGRAM suutikeisann
IMPLICIT REAL(A-H,O-Z)
EXTERNAL FUN1,DIVFUN1
DIMENSION X_D(50,3)

DATA EPS/1.E-7/

AA=0.2
BB=0.3
CC=1.0

WRITE(*,*) ' '
WRITE(*,*) '  N          X2           F(X2)            REL_ERROR '

AA=0.2
BB=0.3
CC=1.0
XX=1.3
X1=XX
DO N=1,100
	CC=1.0
	X2=X1-CC*FUN1(X1,AA,BB)/DIVFUN1(X1,AA)
	
	DO WHILE(ABS(FUN1(X2,AA,BB))>=(1-CC/4)*ABS(FUN1(X1,AA,BB)))
		CC=CC*0.5
	ENDDO

	EPSD=ABS(X2-X1)
	EPSX=ABS(EPSD/X1)
	
	WRITE(*,120) N,X2,FUN1(X2,AA,BB),EPSX
120		FORMAT(' ',I3,5E17.7)
		
	IF(EPSX.LT.EPS) EXIT
	X1=X2
	
ENDDO

WRITE(*,*) ' '

END

REAL FUNCTION FUN1(X,AA,BB)

IMPLICIT REAL(A-H,O-Z)

FUN1=tanh(X)+AA*X+BB

RETURN
END

REAL FUNCTION DIVFUN1(X,AA)

IMPLICIT REAL(A-H,O-Z)

DIVFUN1=(1/cosh(X))**2+AA

RETURN
END

試したこと

周りの人にアドバイスを求めたところ、DO WHILEの部分が終わらないので永遠に計算が続くとのことですが、どのようにコードを直せばいいのかわかりません。

Comments

[chirimen]

やろうとしている減速法とやらの内容を説明してください。
たぶんだけど X2 は do while ループの中で再計算しないといけないとかじゃないですか？

do while ループの中で、X1, X2, CC などを表示させて数値の動きを確認するとよいです。

[Yoshi1416]

コメントありがとうございます。減速法というのは解近傍の初期値x0から出発して、減速パラメータCC
(CCはABS(FUN())<(1-CC/4)*ABS(FUN())の条件を満たす)
のついた反復計算X2=X1-CC*FUN()/DIVFUN()を繰り返すことで近似解へ収束させるという計算方法です。

[Yoshi1416]

DO WHILE(ABS(FUN1(X2,AA,BB))>=(1-CC/4)*ABS(FUN1(X1,AA,BB)))の部分を
>=　→　>　へと変えたところ解決しました。ありがとうございました。

[chirimen]

DO WHILE(ABS(FUN1(X2,AA,BB))>=(1-CC/4)*ABS(FUN1(X1,AA,BB)))
のところで、FUN1() の値がループで変化していません。

これだと CC が単純な1次不等式を解くことで求めることができて、
0<CC/4<1 を満たす CC がない場合がでてきます。

[Yoshi1416]

ありがとうございます。修正してみます。

Answer 1

DO WHILE(ABS(FUN1(X2,AA,BB))>=(1-CC/4)*ABS(FUN1(X1,AA,BB)))の部分を>=　→　>　へと変えたところ解決しました。