Метод Ньютона для нелинейного набора уравнений

Question

Я новичок в программировании в Fortran90. Я использовал метод NR для системы нелинейных уравнений, найденных в Numerical Recipes, и собрал код, который не генерирует никаких ошибок при компиляции с 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

Вы не предоставляют значения для ntrial, tolx или tolf в вызове mnewt. Какие значения вы хотите использовать алгоритм?

Ваше первоначальное предположение с x(1)=0.0 приводит к делению на ноль при вычислении fjac(1,1)=((x(1))**(-0.5)).

Ваши массивы 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) 

Running свою версию пишет это как выход:

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 Fortran с включенной полной проверкой времени выполнения.