Newton-Raphson Method and the Van der Waal Equation Coding question

Question

So I need to find Specific Volumes using the Newton Method from the Van der Waal equation over a change of constant Temperature but variant Pressure.

I have been able to make a list of however many iterations of the altered Van der Waal equation for the root finding method from Pressure Min to Pressure Max (3.9atm to 59.2atm). Further, I am able to select one element/expression from that list and find the root/Specific Volume for that iteration.

However, I cannot seem to figure out how to make a list of all the Specific Volumes, ie solve for the root in more than one expression at a time.

This is my code solving for the 10th element in the list (about 12.2atm):

---

a = 1.36;

b = .003183;

R = .0820578;

T = 333;

inc = (59.2 - 3.9)/60;

c = Table[ vi - ((((atm) + (a/vi^2)) (vi - b) - (R*T))/((atm) - (a/ vi^2) + (2 a*b/vi^3))), {atm, 3.9, 59.2, inc}]

FindRoot[c[[10]] == 0, {vi, 4}, WorkingPrecision -> 20]

---

The output:

---

{-((-27.3252 + (3.9 + 1.36/vi^2) (-0.003183 + vi))/(
   3.9 + 0.00865776/vi^3 - 1.36/vi^2)) + 
  vi,

 -((-27.3252 + (4.82167 + 1.36/vi^2) (-0.003183 + vi))/(
   4.82167 + 0.00865776/vi^3 - 1.36/vi^2)) + 
  vi, 

-((-27.3252 + (5.74333 + 1.36/vi^2) (-0.003183 + vi))/(
   5.74333 + 0.00865776/vi^3 - 1.36/vi^2)) + vi,

AND SO ON.............

Also:

---

{vi -> 2.1930639624928306932}

---

Which is a correct answer.

I've tried this fix to solve for more than one element in the list(elements one to ten):

---

FindRoot[c[[1;;10]] == 0, {vi, 4}, WorkingPrecision -> 20]

---

but this yields an error:

---

FindRoot::nveq: The number of equations does not match the number of variables in FindRoot[c[[1;;10]]==0,{vi,4},WorkingPrecision->20]. >>

---

Any ideas? Thanks!

EDIT: I'm sorry, it's working again for the Van der Waals but neither of the other 2. No matter how I change the starting point the output now is the same, odd number repeated and I don't know how to manipulate either FindMinimum nor Foldlist to change that.

My code is

R = 8.314;
A = 507.28;
B = 0.10476;
a = 0.07132;
b = 0.07235;
c = 660000;
T = 333;
iter = 60;

FX = -P - (A (1 - a/v))/v^2 + (R T (1 - c/(T^3 v)) (B (1 - b/v) + v))/
   v^2;
FXX = -((a A)/v^4) + (2 A (1 - a/v))/v^3 + (
   R T (1 + (b B)/v^2) (1 - c/(T^3 v)))/v^2 + (
   c R (B (1 - b/v) + v))/(T^2 v^4) - (
   2 R T (1 - c/(T^3 v)) (B (1 - b/v) + v))/v^3;

FXXX = v - FX/FXX;

nn = Table[FXXX, {P, 400, 6000, iter}];

Ps = Table[P, {P, 400, 6000, iter}]

vs = Rest@FoldList[v /. Last@FindMinimum[#2, {v, #1}] &, 12, nn]

ListPlot[Thread[{Ps, vs}], Joined -> True, PlotRange -> {0, 7}]

---

The output is

{1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77,
  1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 
 1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77,
  1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 
 1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77,
  1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 
 1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77,
  1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 
 1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77,
  1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 
 1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77,
  1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 
 1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77,
  1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 
 1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77,
  1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 
 1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77,
  1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 
 1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77,
  1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 
 1.3665*10^77, 1.3665*10^77, 1.3665*10^77, 1.3665*10^77}

---

I'm sorry I basically know nothing haha, but I need this done and I just don't know how to do it. Thanks so much for your help.

Answer 1

To solve a list of equations separately, you could just use Table again:

Table[FindRoot[c[[i]] == 0, {vi, 4}], {i, Length@c}]

Or without an index:

Table[FindRoot[eq == 0, {vi, 4}], {eq, c}]

Or you could Map a function over your list:

FindRoot[# == 0, {vi, 4}]& /@ c

But I'm not so sure it's the root you are after. If you plot one of your functions:

Plot[c[[10]], {vi, 0, 10}]

What you are finding is actually a local minimum:

because the root finding is getting stuck there.

If that's what you are really after, you should probably be using FindMinimum instead of FindRoot. In any case, you are probably going to have problems with the initial starting point.

Edit: After a bit of messing around, perhaps this is what you are after:

atms = Table[atm, {atm, 3.9, 59.2, inc}];
vis = Table[vi /. Last@FindMinimum[c[[i]], {vi, 7 i^(-1/2) - 0.5}], {i, Length@c}];
ListPlot[Thread[{atms, vis}], Joined -> True, PlotRange -> {0, 7}]

Edit 2: A better way to apply starting conditions is just to use the result from the previous minimum. Starting with 7 for the first minimum, FoldList allows you to pass this on to each successive minimisation.

vis = Rest@FoldList[vi /. Last@FindMinimum[#2, {vi, #1}] &, 7, c];