Newton-Raphson Solver

We are going to reuse all of the code from the beginning of Working with numbers. Here we discuss how to code up the Newton-Raphson method for finding the roots of the cubic van der Waals equation of state problem.

.....
  ! NOTE: everything above here is exactly the same in 16-22a.f90 and 16-22b.f90

  write(*,*) "V=?  (type the number and press <enter>)"
  read (*,*)  V            ! read initial guess after you look at the plot
  write (*,*) '# initial guess: V=',V

  write (*,*) '# n','x_n','f(x_n)','fp(x_n)'
  do i=1,20
     f=c3*(V**3) + c2*(V**2) + c1*V + c0
     fp = 3.0*c3*(V**2) + 2.0*c2*V + c1     ! df/dV worked out analytically
     write (*,*) i,V,f,fp
     V = V-f/fp
  end do

  write (*,*) '#V = ',V, '(dm^3/mol)'
  write (*,*) 'molar density = ',(1/V),'mol/dm^3'
  write (*,*) 'mass density = ',(mm/V),'g/dm^3'

---

As usual, let's take it line-by-line:

  read (*,*)  V            ! read initial guess after you look at the plot

reads from the standard input a value and puts it into the variable V

 write (*,*) '# initial guess: V=',V

  write (*,*) '# n','x_n','f(x_n)','fp(x_n)'

prints out that value again, just to confirm that we got it right, then prints out a heading for the chart of numbers will will generate in our loop

  do i=1,20
     f=c3*(V**3) + c2*(V**2) + c1*V + c0
     fp = 3.0*c3*(V**2) + 2.0*c2*V + c1     ! df/dV worked out analytically
     write (*,*) i,V,f,fp
     V = V-f/fp
  end do

Loops over the variable i from 1 to 20. When the loop variable is an integer, the default is to increment it by 1. Next, we calculate f and fp, as described in Math Chapter G. (fp is the derivative of f). We print out the values, and then in the last line calculate our new value for V. The new value for "V" is then computed as shown in eqn (G.1) in McQuarrie and Simon.

  write (*,*) '#V = ',V, '(dm^3/mol)'
  write (*,*) 'molar density = ',(1/V),'mol/dm^3'
  write (*,*) 'mass density = ',(mm/V),'g/dm^3'

Finally, having reached a converged value for "V", we print it out, as well as (for our convenience), print out the molar density and mass density.

---

The smallest root corresponds to the molar volume of the liquid, and the largest root is the molar volume of the vapor. I used the plot that we computed in Working with numbers to make my initial guess for the values of these "zeros" of the vdW equation. Compiling and running the program above (and typing in "0.0" or "0.3" when appropriate, so that the read statement we discussed above can "get" the value), here is my output:

For the liquid

jschrier@ws7:~/pchem1/ps1$ ./a.out
 #conditions =     35.00000        142.6900    
 #vdW eqn: (    1.000000     ) v^3 + (  -0.3663687     ) v^2 + ( 
   3.8020000E-02 ) v + (  -1.2101767E-03 )
0.0
 # initial guess: V=    0.000000    
 # nx_nf(x_n)fp(x_n)
            1    0.000000      -1.2101767E-03   3.8020000E-02
            2   3.1830005E-02  -3.3893756E-04   1.7736409E-02
            3   5.0939709E-02  -9.1941329E-05   8.4791277E-03
            4   6.1782964E-02  -2.3833243E-05   4.2007118E-03
            5   6.7456581E-02  -5.6446297E-06   2.2432059E-03
            6   6.9972903E-02  -1.0224758E-06   1.4368519E-03
            7   7.0684507E-02  -7.8929588E-08   1.2157112E-03
            8   7.0749432E-02  -5.8207661E-10   1.1956841E-03
            9   7.0749916E-02    0.000000       1.1955388E-03
           10   7.0749916E-02    0.000000       1.1955388E-03
           11   7.0749916E-02    0.000000       1.1955388E-03
           12   7.0749916E-02    0.000000       1.1955388E-03
           13   7.0749916E-02    0.000000       1.1955388E-03
           14   7.0749916E-02    0.000000       1.1955388E-03
           15   7.0749916E-02    0.000000       1.1955388E-03
           16   7.0749916E-02    0.000000       1.1955388E-03
           17   7.0749916E-02    0.000000       1.1955388E-03
           18   7.0749916E-02    0.000000       1.1955388E-03
           19   7.0749916E-02    0.000000       1.1955388E-03
           20   7.0749916E-02    0.000000       1.1955388E-03
 #V =    7.0749916E-02 (dm^3/mol)
 molar density =     14.13429     mol/dm^3
 mass density =     564.6367     g/dm^3
FORTRAN STOP

For the vapor

jschrier@ws7:~/pchem1/ps1$ ./a.out
 #conditions =     35.00000        142.6900    
 #vdW eqn: (    1.000000     ) v^3 + (  -0.3663687     ) v^2 + ( 
   3.8020000E-02 ) v + (  -1.2101767E-03 )
0.3
 # initial guess: V=   0.3000000    
 # nx_nf(x_n)fp(x_n)
            1   0.3000000       4.2226380E-03   8.8198751E-02
            2   0.2521236       1.1134231E-03   4.3978527E-02
            3   0.2268062       2.3375126E-04   2.6153758E-02
            4   0.2178686       2.4371315E-05   2.0779714E-02
            5   0.2166958       3.9441511E-07   2.0110056E-02
            6   0.2166762      -4.6566129E-10   2.0098940E-02
            7   0.2166762       4.6566129E-10   2.0098954E-02
            8   0.2166762      -4.6566129E-10   2.0098940E-02
            9   0.2166762       4.6566129E-10   2.0098954E-02
           10   0.2166762      -4.6566129E-10   2.0098940E-02
           11   0.2166762       4.6566129E-10   2.0098954E-02
           12   0.2166762      -4.6566129E-10   2.0098940E-02
           13   0.2166762       4.6566129E-10   2.0098954E-02
           14   0.2166762      -4.6566129E-10   2.0098940E-02
           15   0.2166762       4.6566129E-10   2.0098954E-02
           16   0.2166762      -4.6566129E-10   2.0098940E-02
           17   0.2166762       4.6566129E-10   2.0098954E-02
           18   0.2166762      -4.6566129E-10   2.0098940E-02
           19   0.2166762       4.6566129E-10   2.0098954E-02
           20   0.2166762      -4.6566129E-10   2.0098940E-02
 #V =    0.2166762     (dm^3/mol)
 molar density =     4.615182     mol/dm^3
 mass density =     184.3673     g/dm^3
FORTRAN STOP