Newton 방정식 비선형 방정식

Question

저는 Fortran90에서 프로그래밍하는 신참입니다. Numerical Recipes에있는 비선형 방정식 시스템에 대해 NR 메서드를 사용하고 GFortran으로 컴파일 할 때 오류를 생성하지 않는 코드를 조합했습니다. 문제는 출력에 대한 값도 생성되지 않는다는 것입니다. 내 초기 추측 루트 값이 실제 루트에서 멀리 떨어져 있거나이 코드에서 오류가 발생했기 때문일 수 있습니까?

이 문제에 대한 도움이나 조언을 주시면 감사하겠습니다.

program main 
    implicit real*8(a-h,o-z) 
    parameter(n=4) 
    logical check 
    real*8 x(n),fvec(n),fjac(n) 
    open(20,file="output1.txt",status="unknown") 
    do i=1,n 
     x(i)= 4. 
    enddo 
    call mnewt(ntrial,x,n,tolx,tolf)  
     call usrfun(x,n,fvec,fjac) 
    do i=1,n 
     write(20,*) 'x(',i,')=',x(i),fvec(i) 
    enddo 
    end 

subroutine mnewt(ntrial,x,n,tolx,tolf) 

integer n,ntrial,np 
real*8 tolf,tolx,x(n) 
parameter (np=15) 
!uses lubksb, ludcmp, usrfun 
! Given an initial guess x for a root in n dimensions, take ntrial Newton-Raphson steps to 
! improve the root. Stop if the root converges in either summed absolute variable increments 
! tolx or summed absolute function values tolf. 
integer i,k,indx(np) 
real*8 d,errf,errx,fjac(np,np),fvec(np),p(np) 

do 14 k=1,ntrial 
    call usrfun(x,n,fvec,fjac) 
    errf=0. 
    do 11 i=1,n 
     errf=errf+abs(fvec(i)) 
    11 continue 
    if(errf.le.tolf)return 
    do 12 i=1,n 
     p(i)=-fvec(i) 
    12 continue 

    call ludcmp(fjac,n,np,indx,d) 
    call lubksb(fjac,n,np,indx,p) 

    errx=0. 
    do 13 i=1,n 
    errx=errx+abs(p(i)) 
    x(i)=x(i)+p(i) 
    13 continue 
    if(errx.le.tolx)return 

14 continue 
return 
end 

subroutine usrfun(x,n,fvec,fjac) 

    implicit none 

    integer n 
    real*8 x(n),fvec(n),fjac(n,n), hl, ul, br, bl 

hl=1.00 
ul=1.00 
br=0.20 
bl=0.00 

! Initial guesses 

     x(1)=0.0 
     x(2)=1.5 
     x(3)=0.5 
     x(4)=0.5 

    fvec(1)=(x(2))+(2*sqrt((x(1))))-ul-(2*(sqrt(hl))) 
    fvec(2)=((x(3))*(x(4)))-((x(1))*(x(2))) 
    fvec(3)=((x(3))*(x(4))*(x(4)))+(0.5*(x(3))*(x(3)))-((x(1))*(x(2))*(x(2)))-(0.5*(x(1))*(x(1)))+(0.5*(br-bl)*x(1)+x(3)) 
    fvec(4)=(x(4))-sqrt((x(3))) 

    fjac(1,1)=((x(1))**(-0.5)) 
    fjac(1,2)=1 
    fjac(1,3)=0 
    fjac(1,4)=0 
    fjac(2,1)=(-x(2)) 
    fjac(2,2)=(-x(1)) 
    fjac(2,3)=x(4) 
    fjac(2,4)=x(3) 
    fjac(3,1)=((x(2))**2)-(x(1))+(0.5)*(br-bl) 
    fjac(3,2)=-2*((x(1))*(x(2))) 
    fjac(3,3)=((x(4))*(x(4)))+(x(3))+(0.5)*(br-bl)*(x(3)) 
    fjac(3,4)=2*((x(3))*(x(4))) 
    fjac(4,1)=0 
    fjac(4,2)=0 
    fjac(4,3)=-0.5*((x(3))**(-0.5)) 
    fjac(4,4)=1 

end subroutine usrfun 


    subroutine ludcmp(a,n,np,indx,d) !fjac=a 
     integer n,np,indx(n),nmax 
     real*8 d,a(np,np),tiny 
     parameter (nmax=2500,tiny=1.0e-20) 
     integer i,imax,j,k 
     real*8 aamax,dum,sum,vv(nmax) 

     d=1. 
     do 12 i=1,n 
     aamax=0. 
     do 11 j=1,n 
      if (abs(a(i,j)).gt.aamax) aamax=abs(a(i,j)) 
! print*,a(21,j) 
11  continue 
!  print*,i,aamax 
! pause 
     if (aamax.eq.0.) pause 'singular matrix in ludcmp' 
     vv(i)=1./aamax 
12 continue 
     do 19 j=1,n 
     do 14 i=1,j-1 
      sum=a(i,j) 
      do 13 k=1,i-1 
      sum=sum-a(i,k)*a(k,j) 
13  continue 
      a(i,j)=sum 
14  continue 
     aamax=0. 
     do 16 i=j,n 

      sum=a(i,j) 
      do 15 k=1,j-1 
      sum=sum-a(i,k)*a(k,j) 
15  continue 
      a(i,j)=sum 
      dum=vv(i)*abs(sum) 
      if (dum.ge.aamax) then 
      imax=i 
      aamax=dum 
      endif 
16  continue 
     if (j.ne.imax)then 
      do 17 k=1,n 
      dum=a(imax,k) 
      a(imax,k)=a(j,k) 
      a(j,k)=dum 
17  continue 
      d=-d 
      vv(imax)=vv(j) 
     endif 
     indx(j)=imax 
     if(a(j,j).eq.0.)a(j,j)=tiny 
     if(j.ne.n)then 
      dum=1./a(j,j) 

      do 18 i=j+1,n 
      a(i,j)=a(i,j)*dum 
18  continue 
     endif 
19 continue 
     return 
     end 

!lubksb 

    subroutine lubksb(a,n,np,indx,b) 
     integer n,np,indx(n) 
     real*8 a(np,np),b(n) 
     integer i,ii,j,ll 
     real*8 sum 
     ii=0 
     do 12 i=1,n 
     ll=indx(i) 
     sum=b(ll) 
     b(ll)=b(i) 
     if (ii.ne.0)then 
      do 11 j=ii,i-1 
      sum=sum-a(i,j)*b(j) 
11  continue 
     else if (sum.ne.0.) then 
      ii=i 
     endif 
     b(i)=sum 
12 continue 
     do 14 i=n,1,-1 
     sum=b(i) 
     do 13 j=i+1,n 
      sum=sum-a(i,j)*b(j) 
13  continue 
     b(i)=sum/a(i,i) 
14 continue 
     return 
     end 

Answer 1

당신은 mnewt로 호출 ntrial, tolx 또는 tolf에 대한 어떤 값을 제공하지 않습니다. 알고리즘에 어떤 값을 사용 하시겠습니까?

을 계산할 때 x(1)=0.0의 초기 추측은 0으로 나눕니다.

배열 fjac 및 a은 전체 크기가 잘못 표시됩니다. 확실히 [n,n] 모양이 있어야합니까?

> diff mine.f90 yours.f90 
5c5 
<  real*8 x(n),fvec(n),fjac(n,n) 
--- 
>  real*8 x(n),fvec(n),fjac(n) 
10c10 
<  call mnewt(10,x,n,1.d-6,1.d-6)  
--- 
>  call mnewt(ntrial,x,n,tolx,tolf)  
68c68 
<   x(1)=0.1 
--- 
>   x(1)=0.0 
100c100 
<  real*8 d,a(n,n),tiny 
--- 
>  real*8 d,a(np,np),tiny 
165c165 
<  real*8 a(n,n),b(n) 
--- 
>  real*8 a(np,np),b(n) 

내 버전을 실행이 출력으로 기록 :

x(1)= 0.1000000014901161 -0.8675444632541631 
x(2)= 1.5000000000000000 9.9999997764825821E-02 
x(3)= 0.5000000000000000 0.5299999967962503 
x(4)= 0.5000000000000000 -0.2071067811865476 

당신이 무엇을 기대인가요

나는 코드를 다음과 같이 변경했다?

전체 런타임 검사가 활성화 된 NAG 포트란 컴파일러를 사용하여 컴파일하여 이러한 문제를 진단하는 데 약 5 분이 소요되었습니다.