Newton-Raphson Method and the Van der Waal Equation Coding 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[] == 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.

• May I suggest you to rewrite the equation of state in terms of the compressibility factor (Z=PV/RT). This reduces the equation of state to a cubic ploynomial the roots of which are easy to obtain. – Claude Leibovici Oct 17 '13 at 10:44

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[], {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];
• This works! Thanks so much. How did you find the starting condition? The 7i^-.5 - .5 ? I'm trying to apply this method to two other equations, Ideal Gas Law and the Beattie Bridgeman but I can't get it working. – scootingblenny12 Oct 17 '13 at 3:12
• Starting conditions were found by trial and error: minimise with a starting value function, plot the results, make up a new function that approximates the good bits of the plot, repeat. See the edit for a better way to do it. – wxffles Oct 17 '13 at 3:27
• 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. – scootingblenny12 Oct 17 '13 at 15:29
• I've made an edit with the new code and output for the Beattie-Bridgeman – scootingblenny12 Oct 17 '13 at 15:33
• Looks like you actually want a root this time instead of the minimum. Change FindMinimum back to FindRoot. Also don't forget to adjust your PlotRange. – wxffles Oct 17 '13 at 19:24