Programming Example 12
Red characters are codes I wrote by myself today. 2018/01/05
But these can be changed like below using "matmul".
---------------------------------------------------------------------
do while ( (err>TOL).and.(k<k_max) )
k= k + 1
H=Hessian(x) !Get Hessian matrix as H using Hessian Function
Minv=Minverse(H) !Get inverse of Hessian matrix using Minverse Function
x=x-matmul(Minv,R) ! Update x
R = Residual(x) !Set Residual matrix as R.
err=sqrt(dot_product(R,R))
write(*,'(I10, $)') k
do i=1, m
write(*,'(F10.3, $)') x(i)
end do
write(*,'(E20.3)') err
end do
---------------------------------------------------------------------
This is the main code.
------------------------------------------------------------------------
Program nonlinear_system
! Global variables
integer, parameter:: m=3 ! Maximum number of variables
real(8):: S(m), H(m,m) ! Unknown variable vector S
real(8), parameter:: TOL=1D-6 ! Numerical tolerance
Call Guess(S)
Call Solve(S)
contains
Function Residual(x)
real(8):: Residual(m), x(m)
Residual(1)=-x(1)+2*x(2)+x(3)-2
Residual(2)=2*( (x(1)+3)**2) + 5*x(2)**2- (x(3)+1)**2
Residual(3)=(x(1)-5)*(x(2)+2) + 3
end Function Residual
Function Hessian(x)
real(8):: Hessian(m,m), x(m), x_(m), v(m), v_(m)
integer:: i,j
do i = 1, m
do j = 1, m
x_ = x; x_(j) = x(j) + TOL ! Perturb x
v_ = Residual(x_)
v = Residual(x)
Hessian(i,j) = ( v_(i) - v(i) ) / TOL
end do
end do
end Function Hessian
Subroutine Guess(x)
real(8), INTENT(OUT):: x(m) ! dummy variables
Integer:: i
Write(*,*) "Nonlinear system of", m, "equations"
do i = 1, m
write(*, '(A, I2, A)', advance="no") "Input the guessed value no.", i, " = "
Read(*,*) x(i)
end do
End Subroutine Guess
Subroutine Solve(x)
real(8), INTENT(INOUT):: x(m) ! dummy variables
! local variables
Integer, parameter:: k_max=50 ! Maximum iteration number
Integer:: k, n
real(8):: err, R(m), H(m,m),Minv(m,m),rs(m)
write(*,'(A10, $)') "Iteration"
do i=1, m
write(*,'(A9, I1,$)') "x", i
end do
write(*,'(A20)') "Norm of errors"
write(*, '(A)') repeat('-',60)
R = Residual(x)
err=sqrt(dot_product(R,R))
k=0
write(*,'(I10, $)') k
do i=1, m
write(*,'(F10.3, $)') x(i)
end do
write(*,'(E20.3)') ERR
do while ( (err>TOL).and.(k<k_max) )
k= k + 1
H=Hessian(x) !Get Hessian matrix as H using Hessian Function
Minv=Minverse(H) !Get inverse of Hessian matrix using Minverse Function
R = Residual(x) !Set Residual matrix as R.
rs=0.0
do i =1,m
do j =1,m
rs(i)=rs(i)+Minv(i,j)*R(j) !Get H(x_k)[-1]*r(x_k)
end do
x(i)=x(i)-rs(i) !Set New x
end do
R = Residual(x) !Set Residual matrix as R again.
err=sqrt(dot_product(R,R))
write(*,'(I10, $)') k
do i=1, m
write(*,'(F10.3, $)') x(i)
end do
write(*,'(E20.3)') err
end do
if (k==k_max) then
print*, "Divergence!"
end if
write(*, '(A)') repeat('-',60)
End Subroutine Solve
function Minverse(A)
real(8), dimension(:,:):: A
real(8), dimension(:,:), allocatable:: Minverse, B, V
real(8), dimension(:), allocatable:: temp
real(8), parameter:: TOL=1D-6
real(8):: pivot
integer:: i, j, jmax, k, n,l
n=size(A)**0.5
allocate(Minverse(n,n),B(n,n),V(n,n),temp(n))
B=A
V(:,:)=0.0
do i=1, n
V(i,i)=1.0
end do
do i=1, n
jmax=i
do j=i+1, n
if (abs(B(j,i))>abs(B(jmax,i))) then
jmax=j
temp=B(i,:);B(i,:)=B(jmax,:);B(jmax,:)=temp
temp=V(i,:);V(i,:)=V(jmax,:);V(jmax,:)=temp
end if
end do
pivot=B(i,i)
if (abs(pivot)<TOL) then
print *, "The matrix is singular!"
end if
B(i,:)=B(i,:)/pivot
V(i,:)=V(i,:)/pivot
do k =1,n
if (i .ne. k) then
pivot=B(k,i)
do l=1,n
B(k,l)=B(k,l)-B(i,l)*pivot
V(k,l)=V(k,l)-V(i,l)*pivot
end do
end if
enddo
end do
Minverse=V
end function Minverse
end Program nonlinear_system
