Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Newton Raphson method (convergence problem)

  1. Nov 26, 2015 #1
    Hello, I am new in the field of FORTRAN. I started to write a code using Newton Raphson method. Below is the code of main program portion only. During my calculation, I have found that my results are not converging, Any help or advice will be highly appreciated in this case.
    Thank you for the effort of helping in advance.


    Program Free_energy
    Implicit None
    Integer:: i,j,l,k
    Double Precision Gmic,U,Phi
    Double Precision s,Uderv,cc,DG0,cf,ctot1,
    & Keff,Cagg,ctot,sum,Cagg1,cnew,
    & Cagg2,cctot1,kapp,error,Csalt,Csalt1,
    & cext
    Double precision :: Nelec
    Double precision :: Nwat
    Double precision :: Nlig

    ctot = 1D0
    cf=ctot/2
    cext = 0.01

    sum=0.0
    DG0 = -100.0

    Do While(DG0.le.20)

    do i=0,2
    do j=0,10
    do l=0,10

    Nelec=i
    Nwat=j
    Nlig=l

    Keff=dexp(-Gmic(Nelec,Nwat,Nlig)-Nelec*DG0)
    Cagg=Keff*((cf)**Nlig)
    Cagg1=sum+Nlig*Cagg
    Cagg2=sum+(Nlig**2)*Cagg

    enddo
    enddo
    enddo

    ctot1=ctot-Cagg1
    cctot1=ctot-Cagg2

    cnew = cf-(ctot1/cctot1)

    error = abs((cnew-cf)/cnew)

    cf = cnew

    write(72,*) cf, error
    if(error.lt.1D-7) exit

    Csalt1 = 0.0

    do i=0,2
    do j=0,10
    do l=0,10

    Nelec=i
    Nwat=j
    Nlig=l

    Keff = dexp(-Gmic(Nelec,Nwat,Nlig))
    Csalt = Csalt1+Keff*((cf)**Nlig)*Nelec

    enddo
    enddo
    enddo

    Kapp = Csalt/cext

    write(70,*) DG0, -log(kapp)/2
    write(72,*) cf, -log(kapp)/2

    DG0 = DG0 + 0.1

    enddo
     
  2. jcsd
  3. Nov 26, 2015 #2
    I find that graphing the results of the iterations can be helpful.
     
  4. Nov 26, 2015 #3

    Dr Transport

    User Avatar
    Science Advisor
    Gold Member

    Up the number of iterations and check for their convergence......it seems to me at first glance that you are not iterating enough.
     
  5. Nov 27, 2015 #4
    First I would like to thank you both for the reply,

    Since you have mentioned that increasing the iteration number may solve the problem, but my question is how can I do that?

    Also in my code, the function has summation over an exponential series of the Gimc, I am afraid of that also because I think that summation is also not happening there properly.
    Further, I made two changes in the expression of the derivative of my function (Cagg2) and the expression of Keff in my code. The rest of the code is unchangeable and almost like the previous one except that two lines. I ran this one and I find this plot (-log(Kapp)/2 vs DG0) attached with this reply. My aim was to calculate the "cf" which has also the negative values in the calculated results.

    I would like to thank you again for the help. Any advice or direction in this case will become helpful to me in near future.

    Note: For your interest, the new corrected version of the code (only the main program portion only) has been given below,

    Program Free_Energy
    Implicit None
    Integer:: i,j,l,k
    Double Precision Gmic,U,Phi
    Double Precision s,Uderv,cc,DG0,cf,ctot1,
    & Keff,Cagg,ctot,sum,Cagg1,cnew,
    & Cagg2,cctot1,kapp,error,Csalt,Csalt1,
    & cext
    Double precision :: Nelec,kB,T
    Double precision :: Nwat
    Double precision :: Nlig

    ctot = 1D0
    cf=ctot/2
    cext = 0.01

    DG0 = -100.0
    Do While(DG0.le.20)

    do i=0,2
    do j=0,10
    do l=0,10

    Nelec=i
    Nwat=j
    Nlig=l

    Keff=dexp(-Gmic(Nelec,Nwat,Nlig)-Nelec*DG0)
    Cagg=Keff*((cf)**Nlig)
    Cagg1=sum+Nlig*Cagg
    Cagg2=sum+(Nlig**2)*Keff*((cf)**(Nlig-1))


    enddo
    enddo
    enddo

    ctot1=ctot - Cagg1
    cctot1=ctot - Cagg2
    cnew = cf-(ctot1/cctot1)

    error = abs((cnew-cf)/cnew)

    cf = cnew

    write(74,*) cf,error

    if(error<1D-7)exit

    Csalt1 = 0.0

    do i=0,2
    do j=0,10
    do l=0,10

    Nelec=i
    Nwat=j
    Nlig=l

    Keff = dexp(-Gmic(Nelec,Nwat,Nlig))
    Csalt = Csalt1+Keff*((cf)**Nlig)*Nelec

    write(73,*) cf
    enddo
    enddo
    enddo

    Kapp = Csalt/cext

    write(70,*) DG0, -log(Kapp)/2

    write(72,*) cf, -log(Kapp)/2

    DG0 = DG0 + 0.1

    enddo
    End
     

    Attached Files:

  6. Nov 28, 2015 #5

    Dr Transport

    User Avatar
    Science Advisor
    Gold Member

    part of your issue is that for example you define Nwat as a double, then assign an integer to it, yes it theoretically works, but better programming style is

    Nwat = double(j)

    it makes a double out of the j so that it is j.00000000
     
  7. Nov 29, 2015 #6
    Thank you for the reply, I will follow your suggestion and let you know about the result.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Newton Raphson method (convergence problem)
  1. Newton Method (Replies: 1)

Loading...