Several Issues with Van der waals function used to plot P against V

Primary issue- With the for function- Illegal use of reserved keyword "for".

Matlab Script:

function [P,V]=ISOTHERMS_VDWTest(T,Pc,Tc,vmax)
  R=8.314;
  b=(R*Tc)/(8*Pc);
  a=(27/64)*(R^2*Tc^2)/Pc;
  v(i+1)=v(i)+v(i)*0.01
  v(1)=b+b*0.01
  i=1
  while v(i)<vmax
    v(i+1)=v(i)+v(i)*0.01
    i=i+1
  end
  v=[v(1),v(2),.....v(i)]
  for J=i=1
    P(J)=(R*T/V(J)-b)-(a/V(J)^2)
  end
  P=P/1000000
  v=v*1000
  plot(v,P)
end

********************

Any other issues or easier way to do it is also appreciated.

[dpb]

for J=i=1

has two "=" signs instead of a colon for a range expression. Not sure what you had in mind???

Also not sure what your intent/thought is here; you already have a vector v that was made by the while loop.

It'll be a real nuisance that you have left of the trailing ; on so many lines to have all that output every time...

v=[v(1),v(2),.....v(i)

Answer 1

[Pieter Hamming]

Part one: problems with code in the current form

v=[v(1),v(2),.....v(i)]

is incorrect syntax, as it is interpreted as

v=[v(1),v(2),
..v(i)]

v is already a row vector, so this row is redundant. Deleting it brings us to the second problem:

for J=i=1

As dpb pointed out, "|for J is i to 1|" has syntax

for J=i:-1:1

Also, i is used before it is defined. Row 5-7 need to be reversed

i=1 %define i
v(1)=b+b*0.01 %define v(1)
v(i+1)=v(i)+v(i)*0.01 %define v(2) based on v(1) and i

At row 14 you use capital V instead of lowercase v.

Fixing all of these gives me a bug-free code.

Part two: easier ways to do the same

please end your lines with a semicolon to suppress output

v can be rewritten:

v(i+1) = 1.01 * v(i); %original
v(i+1) = 1.01^i * v(1); %rewritten

Now, instead of a while loop over i, you can do it in a single calculation:

imax=5; %example
v(1)=b+b*0.01; %define the first entry
v=1.01.^(0:imax) * v(1); %single calculation to do all remaining values of v

I used imax in the above example; but you only have a vmax. Arithmatic to the rescue: we need to find the highest value imax for which v(imax)<vmax holds true.

v(imax) = vmax
1.01.^imax * v(1) = vmax
vmax / v(1) = 1.01.^imax
imax=log(vmax / v(1)) / log(1.01)

So knowing vmax and v(1), you can just calculate imax; no pesky while loops required. The imax provided by this will not be an integer, so use

imax=ceil(imax) % round-up to next integer

Next we have the calculation of P. This can be simplified using matrix-dot-multiplication (or element-wise multiplication).

P=(R*T./v-b)-(a./v.^2);

This does the same thing as the for loop in your code.

Combining all, we get this example code:

function [P,V]=ISOTHERMS_VDWTest(T,Pc,Tc,vmax)
  %original code
  R=8.314;
  b=(R*Tc)/(8*Pc);
  a=(27/64)*(R^2*Tc^2)/Pc;
  v(1)=b+b*0.01;

   %new code
    imax=log(vmax / v(1)) / log(1.01);
    imax=ceil(imax);
    v=1.01.^(0:imax) * v(1);
    P=(R*T./v-b)-(a./v.^2);

    %original code
    P=P/1000000;
    v=v*1000;
    plot(v,P)
  end

Done. Sorry for the long post, but I hope this helps