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```
```