Matlab help, shouldn't be hard...

everettp

Ok everyone, I'm pretty rusty with matlab so I was wondering if I could get some help... I have a program that bascially is defining different parts of an equation and requires the input of a temperature and a pressure (P,T) Rather than individually inputing a pressure at a certain temp and finding the roots I want to put in a range of temperatures, say 200 K - 500K in 25 intervals and for each temperature I want to have it solve for the roots when P goes to 1 to 150 bar, and finally plot the roots. I tried doing the whole 1:1:150 for temp but I'm having problems with the matrix not being sqare and the whole bit. Here is what I have so far. Thanks for the help in advance!

w = 0.345;     %accentric factor
M = 220.55;    %critical pressure (bar)
N = 647.1;     %critical temperature (K)
R = 83.14472;  %bar.cm³/(mol.K)      
T = 400;       %Isotherm
P = input('Pressure:');
k = 0.480 + 1.574*w - 0.17*(w^2);
alpha_t = (1+k*(1-sqrt(T/N)))^2;
a_t = 0.42748*((R*N)^2)*alpha_t/M;
A = a_t*(P/((R*T)^2));
b = 0.08664*((R*N)/M);
B = (b*P)/(R*T);
beta = A-B-B^2;
gamma = -A*B;
Z = [1 -1 beta gamma];  %Root finder for cubic RSK EOS
G = roots(Z);
disp(G)

drasnor

I don't see where you're defining a matrix/vector/array in that code snippet but remember when you're doing a element-wise matrix operation you need to use . operators.

e.g. matrix_P .* matrix_T (multiplies each element of matrix_P by the element of the same index in matrix_T as long as both are the same size).

You'll get the matrix not square error if you omit the dot because then it tries to compute the vector product.

-drasnor

shwaip

you'll need a .^ here too, if T is a vector:

A = a_t*(P/((R*T)^2));

everettp

w = 0.345;     %accentric factor
M = 220.55;    %critical pressure (bar)
N = 647.1;     %critical temperature (K)
R = 83.14472;  %bar.cm³/(mol.K)      
T = 200:25:600;       %Isotherm
P = 1:0.5:150;
k = 0.480 + 1.574*w - 0.17*(w^2);
alpha_t = (1+k*(1-sqrt(T/N))).^2;
a_t = 0.42748*((R*N)^2)*alpha_t/M;
A = a_t*(P/((R*T).^2));
b = 0.08664*((R*N)/M);
B = (b.*P)/(R.*T);
beta = A-B-B^2;
gamma = -A*B;
Z = [1 -1 beta gamma];  %Root finder for cubic RSK EOS
G = roots(Z);
disp(G)

Ok here's an updated code. It still gives me the matrix isn't square. Any ideas?

drasnor

alpha_t = (1+k.*(1-sqrt(T./N))).^2;
a_t = 0.42748.*((R.*N).^2).*alpha_t./M;
A = a_t.*(P/((R.*T).^2));
b = 0.08664.*((R.*N)./M);
B = (b.*P)/(R.*T);
beta = A-B-B.^2;
gamma = -A.*B;

. operators can be used in place of normal operators for single-element operations. Since you're not doing any vector operations just brute force it.

-drasnor

shwaip

it should also give an line # for the error. that should help you narrow it down.

everettp

Heres the updated code along with the error.

w = 0.345;     %accentric factor
M = 220.55;    %critical pressure (bar)
N = 647.1;     %critical temperature (K)
R = 83.14472;  %bar.cm³/(mol.K)      
T = 200:25:600;       %Isotherm
P = 1:0.5:150;
k = 0.480 + 1.574*w - 0.17*(w^2);
alpha_t = (1+k.*(1-sqrt(T./N))).^2;
a_t = 0.42748.*((R.*N).^2).*alpha_t./M;
A = a_t.*(P./((R.*T).^2));
b = 0.08664.*((R.*N)./M);
B = (b.*P)./(R.*T);
beta = A-B-B.^2;
gamma = -A.*B;
Z = [1 -1 beta gamma];  %Root finder for cubic RSK EOS
G = roots(Z);
disp(G)
??? Error using ==> ./
Matrix dimensions must agree.

As a reminder the ultimate goal is for the program to... set a temperature, find all the roots at a range of pressures, then set a new temp, find all the roots at the same range of pressure, etc. Ultimately I want to plot it out.

Note: whenever I get the program to "work" I just get a bunch of imaginary values, I should have A LOT more real values, I've put in seperate pressures individually and I get more reals. I don't think it's proccessing the way I want it to

everettp

Is it possible to run that program for one P, display the results and loop back with a new value of pressure?

shwaip

everettp wrote:

Is it possible to run that program for one P, display the results and loop back with a new value of pressure?

I'll have a good answer for you tonight, after I get off work. Just making sure that you don't need this before then.

everettp

shwaip wrote:

I'll have a good answer for you tonight, after I get off work. Just making sure that you don't need this before then.

Nope I have plenty of time, thank you.

shwaip

w = 0.345;                             %accentric factor
M = 220.55;                            %critical pressure (bar)
N = 647.1;                             %critical temperature (K)
R = 83.14472;                          %bar.cm³/(mol.K)      
T = 200:25:600;                        %Isotherm
P = 1:0.5:150;                         %Pressure Values
P=P';                                  %Transpose to make matrix math nice
k = 0.480 + 1.574*w - 0.17*(w^2);
alpha_t = (1+k*(1-sqrt(T/N))).^2;
a_t = 0.42748*((R*N)^2)*alpha_t/M;
A = P*(a_t.*((R*T).^2).^(-1));
b = 0.08664*((R*N)/M);
B = (b*P)*(R*T).^(-1);
beta = A-B-B.^2;
gamma = -A.*B;
Z=zeros(length(P), 2*length(T)+2);     %Initialize arrays
G=zeros(2*length(T)+1,length(P));
for i=1:length(P)
    Z(i,:) = [1 -1 beta(i,:) gamma(i,:)];  %Root finder for cubic RSK EOS
    G(:,i) = roots(Z(i,:));
end
disp(G)

I'm getting a lot of imaginary values too, but I don't really know what cubic RSK EOS is. The answers are consistant with the values I get when I run single values for pressure through. eg, running P=1 through, then P=1.5 gives me the same values as running P=[1 1.5]. Basically, the matrix math is correct, but I don't know if your formulae are.

everettp

This is definately starting to get somewhere. I am going to mess with my equations some more to make sure they're right. I Appreciate the help!

everettp

I guess my equation is set up wrong. It is the Soave Redlich Kwong equation of state. The graph SHOULD resemble a simple cubic, right now Im getting a mess, with negative values of Z which shouldn't happen except for the imaginaries. Attached is what I am looking for (sorry I used paint lol) and also what I ended up with... which is craziness.

Here's what the new code is. My plot may be wrong. I want to plot Pressure vs the roots.

w = 0.345;                             %accentric factor
M = 220.55;                            %critical pressure (bar)
N = 647.1;                             %critical temperature (K)
R = 83.14472;                          %bar.cm³/(mol.K) 
q1 = input('Temperature Range (1):');
q2 = input('Temperature Range (2):');
p1 = input('Pressure Range (1):');
p2 = input('Pressure Range (2):');
T = q1:10:q2;                                %Isotherm
P = p1:0.5:p2;                         %Pressure Values
P=P';                                  %Transpose to make matrix math nice
k = 0.480 + 1.574*w - (0.17*(w^2));   
alpha_t = (1+k*(1-sqrt(T/N))).^2;
a_t = 0.42748*((R*N)^2)*alpha_t/M;
A = P*(a_t.*((R*T).^2).^(-1));
b = 0.08664*((R*N)/M);
B = (b*P)*(R*T).^(-1);
beta = A-B-B.^2;
gamma = -A.*B;
Z=zeros(length(P), 2*length(T)+2);     %Initialize arrays
G=zeros(2*length(T)+1,length(P));
for i=1:length(P)
    Z(i,:) = [1 -1 beta(i,:) gamma(i,:)];  %Root finder for cubic RSK EOS
    G(:,i) = roots(Z(i,:));
end
disp(G)
plot(G,P)

The code is 100% correct. I compaired values by going to this site:
http://www.cheng.cam.ac.uk/~pjb10/thermo/pure.html

shwaip

I know what the problem is - the code for solving RSK is wrong. It's kinda hard to explain, but suffice it to say that:

Z(i,:) = [1 -1 beta(i,:) gamma(i,:)];

isn't doing what you think. It's not solving (1 -1 beta gamma) for each value, its solving (1 -1 beta(1) beta(2)...beta(n) gamma(1) gamma(2) ...gamma(n)), which has weird roots.

The other problem is that you're trying to graph 3 dimensions of info (P,Z,T) using a 2-d graph (P,Z).

You can try a couple things:
1) only do it for one value of T
3) plot it using a 3-d plot

1 takes less work than 3, 3 is cooler than 1, but I prefer 2.