Programming Example 12 日本語解説版（未完）
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
!以降、Residual(x)を使ったときはその時点でのx(1),x(2),x(3)が
!それぞれ代入されて、新しいResiudal(1),(2),(3)が作られます。
!目標はこの方程式を解くこと。つまり、Residual＝０になるx(1),(2),(3)を
!見つけることです。以降の操作で少しづつResidualは０に近づくはずです。
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 ! 1
v_ = Residual(x_) ! 2
v = Residual(x) ! 3
Hessian(i,j) = ( v_(i) - v(i) ) / TOL ! 4
!ここに書かれていることはいたってシンプルで簡単。
!f' = (f(x+TOL)-f(x))/TOL を計算するだけ
!1. x_(1)～x_(3)はx(1)~x(3)に1D-6を足す。
!2. それぞれResidualに入れる。
!3. vとv_を作る
!4. (v_(i)-v(i))/TOL = 微分結果
! 上の微分結果をHessian行列に格納。
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
! 繰り返し作業において、初期値は自分で用意しなきゃいけないので、ここで
! x(1)~(3)の初期値を読み込みます。
End Subroutine Guess
Subroutine Solve(x) !ここはコアの計算部分
real(8), INTENT(INOUT):: x(m) ! dummy variables
Integer, parameter:: k_max=50 ! Maximum iteration number
Integer:: k, n
real(8):: err, R(m), H(m,m),Minv(m,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) !一項目を最初の方程式に入れ、Rに代入。
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) !1. Get Hessian matrix as H using Hessian Function
Minv=Minverse(H) !2. Get inverse of Hessian matrix using Minverse Function
x=x-matmul(Minv,R) !3. Update x
R = Residual(x) !4. Set Residual matrix as R again.
err=sqrt(dot_product(R,R)) !5.Set error.
write(*,'(I10, $)') k !6. あとは書くだけ
do i=1, m
write(*,'(F10.3, $)') x(i)
end do
write(*,'(E20.3)') err
end do
if (k==k_max) then
print*, "Divergence!";stop
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
