!###########################################################
! page.78 表4.3(例題3)
! 2次元ニュートン・ラフソン法のプログラム
!
! e(z)=z**3+6*z**2+21*z+32=0
! f(x,y)=x**3-3*x*y**2+6*x**2-6*y**2+21*x+32=0
! g(x,y)=3*x**2*y-y**3+12*x*y+21*y=0
! z=x+y*i
! e(z)=f(x,y)+g(x,y)*i
!
! df/dx=3*x**2-3*y**2+12*x+21
! df/dy=-6*x*y-12*y
! dg/dx=6*x*y+12*y
! dg/dy=3*x**2-3*y**2+12*x+21
!
!###########################################################
PROGRAM prog15
IMPLICIT REAL(8)(a-h,o-z)
IMPLICIT INTEGER(i-n)

INTEGER,PARAMETER :: KMAX=10
!-------------
x=-1.d0
y=-2.d0
delta=1.d-10
!-------------
fxy=funcf(x,y)
gxy=funcg(x,y)
k=0
WRITE(*,"('k=',i2,', x=',1PE12.5,', y=',1PE12.5,', f=',1PE12.5,', g=',1PE12.5)") k,x,y,fxy,gxy


DO k=1,KMAX

CALL calcJ(x,y,dfdx,dfdy,dgdx,dgdy,fj)

xk=fnewtonx(x,y)
yk=fnewtony(x,y)
fxy=funcf(xk,yk)
gxy=funcg(xk,yk)
WRITE(*,"('k=',i2,', x=',1PE12.5,', y=',1PE12.5,', f=',1PE12.5,', g=',1PE12.5)") k,xk,yk,fxy,gxy

IF(fxy**2+gxy**2<delta**2)GOTO 100

x=xk
y=yk
    
END DO
WRITE(*,*) '解が求まりませんでした'
STOP

100 CONTINUE
WRITE(*,*)
WRITE(*,"('解は z=',1PE12.5,' ',SP,1PE12.5,'*i')") xk,yk
STOP
!###########################################################
CONTAINS

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
FUNCTION funcf(xx,yy)
   REAL(8),INTENT(IN) :: xx,yy
   REAL(8) :: funcf
   funcf=xx**3-3.d0*xx*yy**2+6.d0*xx**2-6.d0*yy**2+21.d0*xx+32.d0
RETURN
END FUNCTION funcf
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
FUNCTION funcg(xx,yy)
   REAL(8),INTENT(IN) :: xx,yy
   REAL(8) :: funcg
   funcg=3.d0*xx**2*yy-yy**3+12.d0*xx*yy+21.d0*yy
RETURN
END FUNCTION funcg
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
FUNCTION fnewtonx(xx,yy)
   REAL(8),INTENT(IN) :: xx,yy
   REAL(8) :: fnewtonx

fnewtonx=xx-(fxy*dgdy-gxy*dfdy)/fj

RETURN
END FUNCTION fnewtonx
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
FUNCTION fnewtony(xx,yy)
   REAL(8),INTENT(IN) :: xx,yy
   REAL(8) :: fnewtony

fnewtony=yy-(gxy*dfdx-fxy*dgdx)/fj

RETURN
END FUNCTION fnewtony
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  

END PROGRAM prog15

! df/dx=3*x**2-3*y**2+12*x+21
! df/dy=-6*x*y-12*y
! dg/dx= 6*x*y+12*y
! dg/dy=3*x**2-3*y**2+12*x+21
SUBROUTINE calcJ(x,y,dfdx,dfdy,dgdx,dgdy,fj)

IMPLICIT REAL(8)(a-h,o-z)
IMPLICIT INTEGER(i-n)

REAL(8),INTENT(IN) :: x,y
REAL(8),INTENT(OUT) :: dfdx,dfdy,dgdx,dgdy,fj

dfdx=3.d0*x**2-3.d0*y**2+12.d0*x+21.d0
dfdy=-6.d0*x*y-12.d0*y

dgdx= 6.d0*x*y+12.d0*y
dgdy=3.d0*x**2-3.d0*y**2+12.d0*x+21.d0

fj=dfdx*dgdy-dfdy*dgdx

RETURN
END SUBROUTINE calcJ
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!