Function handle with integrals of multiple equations?

Thank you for answering 

[David Hill]

Is gamma scalar or is it a function of T? If it is a function of T, then what is the function?

[Deema Khunda]

it is a funciton of T, well, gamma is calculated through the following code :

gamma = Cp/Cp-R
R=8 
Cp= A+B*t+C*(t^2)+D*(t^3)+E/(t^2);
t = T/1000 % where t is a function of T ( temperature) 

gamma changes with T

[Deema Khunda]

apologies, didnt give numerical values for A,B,C,D,E :

A= -0.703029;	
B= 108.4773;
C= -42.52157;
D= 5.862788;
E= 0.678565;	

Thank you

[Ameer Hamza]

what are the values of V1, P1, T1 and V2?

[Torsten]

You don't need Cp since gamma=Cp/Cp-R =1-R :-)

I think you meant gamma = Cp/(Cp-R).

[Torsten]

gamma=Cp/(Cp-1) ?

I thought gamma=Cp/(Cp-R) ...

Answer 1

[Star Strider]

I get different result with a strictly numeric version:

V1 = 230;
P1 = 2.7;
T1 = 300;
V2 = 30; 

A = -0.703029;	
B = 108.4773;
C = -42.52157;
D = 5.862788;
E = 0.678565;	

R=8;
Cp = @(t) A+B*t+C*(t.^2)+D*(t.^3)+E./(t.^2);
gamma = @(t) Cp(t)/(Cp(t)-R);

T = @(t) 1000*t;
Eq1 = @(P2,T2) log(P2)-log(P1) - integral(@(t)(gamma(t)/(gamma(t)-1))./T(t), T1, T2);
Eq2 = @(T2) log(V2)-log(V1) - integral(@(t)(gamma(t)/(gamma(t)-1))./T(t), T1, T2);

B0 =  randi(100,2,1);
B = fsolve(@(b) [Eq1(b(1),b(2)); Eq2(b(2))], B0);                          % Choose Appropriate Valuew For ‘B0’ To Get The Correct Result
fprintf(1, '\nP2 = %8.4f\nT2 = %8.4f\n', B)

producing:

P2 =   0.3522
T2 = 299.9684

EDIT — (4 Apr 2020 at 18:23)

V1 = 230;
P1 = 2.7;
T1 = 300;
V2 = 30;

A = -0.703029;
B = 108.4773;
C = -42.52157;
D = 5.862788;
E = 0.678565;

R=8;
Cp = @(t) A+B*t+C*(t.^2)+D*(t.^3)+E./(t.^2);
gamma = @(t) Cp(t)./(Cp(t)-R);

T = @(t) 1000*t;
Vv = V1:-0.5:V2;
B0 =  randi(500,1,2);
for k = 1:numel(Vv)
    Eq1 = @(P2,T2) log(P2)-log(P1) - integral(@(t)(gamma(t)./(gamma(t)-1))./T(t), T1, T2);
    Eq2 = @(T2) log(Vv(k))-log(V1) - integral(@(t)(gamma(t)./(gamma(t)-1))./T(t), T1, T2);
    
    Prms(k,:) = fsolve(@(b) [Eq1(b(1),b(2)); Eq2(b(2))], B0);                          % Choose Appropriate Valuew For ‘B0’ To Get The Correct Result
    Prms = abs(Prms);
end
Out = table(Vv.', Prms(:,1), Prms(:,2), 'VariableNames',{'V','P2','T2'})

Sample output:

Out =

  401×3 table

      V        P2         T2  
    _____    _______    ______

      230        2.7       300
    229.5     2.6941       300
      229     2.6883       300
    228.5     2.6824       300
      228     2.6765       300
                . . .
       31    0.36391    299.97
     30.5    0.35804    299.97
       30    0.35217    299.97

[Deema Khunda]

Hello

thanks so much for the code, one issue though, its that with decreasing volume we expect both the temperature and pressure to increase not decrease, and increase by large percentage. I have been trying to figure out why your code is giving these values as it is not the correct answer to the equation.

[Star Strider]

I missed the negative sign in the second equation. The coirrect code for ‘Eq2’ should be:

Eq2 = @(T2) log(Vv(k)./V1) + integral(@(t)(gamma(t)./(gamma(t)-1))./T(t), T1, T2);

(dividing the log argument rather than subtracting logs makes the code slightly more efficient), with the loop corrected to:

for k = 1:numel(Vv)
    Eq1 = @(P2,T2) log(P2./P1) - integral(@(t)(gamma(t)./(gamma(t)-1))./T(t), T1, T2);
    Eq2 = @(T2) log(Vv(k)./V1) + integral(@(t)(gamma(t)./(gamma(t)-1))./T(t), T1, T2);
    Prms(k,:) = fsolve(@(b) [Eq1(b(1),b(2)); Eq2(b(2))], B0);                          % Choose Appropriate Valuew For ‘B0’ To Get The Correct Result
    Prms = abs(Prms);
end

and the sample table now being:

Out =

  401×3 table

      V        P2        T2  
    _____    ______    ______

      230       2.7       300
    229.5    2.7059       300
      229    2.7118       300
    228.5    2.7177       300
               . . .
       32    19.406    300.03
     31.5    19.714    300.03
       31    20.032    300.03
     30.5    20.361    300.03
       30      20.7    300.03

My apologies for the initial error.

[Deema Khunda]

Thanks so much, can we try and add how gamma is changing with each loop to the table of V, P2 and T2 as well?

[Star Strider]

My pleasure!

Since ‘gamma’ is part of the integrals, and has only ‘t’ as an argument, and since only ‘T2’ changes, it can be coded as:

gammav = gamma(Prms(:,2)/1000);

since ‘Prms(:,2)’ is ‘T2’ in my code, and included in the ‘Out’ table as:

Out = table(Vv.', Prms(:,1), Prms(:,2), gammav, 'VariableNames',{'V','P2','T2','gamma'})

It doesn¹t change significantly:

Out =
      V        P2        T2      gamma 
    _____    ______    ______    ______

      230       2.7       300    1.2887
    229.5    2.7059       300    1.2887
      229    2.7118       300    1.2887
                    . . .
       31    20.032    300.03    1.2887
     30.5    20.361    300.03    1.2887
       30      20.7    300.03    1.2887

It may change more if viewed in extended precision, however not with the default precision I use.

[Ameer Hamza]

There is still an issue with the equation. In the expression of Eq2, there is no gamma in the numerator.

Eq2 = @(T2) log(Vv(k)./V1) + integral(@(t)(1./(gamma(t)-1))./T(t), T1, T2);
                                          %^ it should be 1

I, too, feel for this trap in my first attempt.

[Star Strider]

Thank you! I misread that.

Correcting it does not change the output in my code.

[Deema Khunda]

Hello Stat Strider,

I have been rediting your code to apply the following equaiton instead: pressure calculation as funciton of volume while temperature is calculated as function of pressure :

I have beeb getting unrealstic results as the following code is applied :

V1 = 230;
P1 = 2.7;
T1 = 300;
V2 = 30;
A = -0.703029;
B = 108.4773;
C = -42.52157;
D = 5.862788;
E = 0.678565;
R=8;
Cp = @(t) A+B*t+C*(t.^2)+D*(t.^3)+E./(t.^2);
gamma = @(t) Cp(t)./(Cp(t)-R);
T = @(t) 1000*t;
Vv = V1:-0.5:V2;
B0 =  randi(500,2,1);
for k = 1:numel(Vv)

    Eq1 = @(P2) log(Vv(k)./V1) + integral(@(t)(1./(gamma(t)))./T(t), P1,P2);
    Eq2 = @(P2,T2) log(P2)-log(P1) - integral(@(t)(gamma(t)./(gamma(t)-1))./T(t), T1, T2)
    Prms(k,:) = fsolve(@(b) [Eq1(b(2)); Eq2(b(1),b(2))], B0);                          % Choose Appropriate Valuew For bB0( To Get The Correct Result
    Prms = abs(Prms);
    gammav = gamma(Prms(:,2)/1000);
end

Out = table(Vv.', Prms(:,1), Prms(:,2), gammav, 'VariableNames',{'V','P2','T2','gamma'})

Would you please be able to tell me what I may have got wrong?

Thank you

[Star Strider]

My pleasure.

This is different than your initial system. Note that ‘Eq1’ is in ‘P’, not ‘T’, so the integrand should be dividing by ’P’ instead of ‘T’. Also, ‘gamma’ is a function of ‘t’ (or ‘T’), not ‘P’, so here it is simply an unintegrated constant, and the first equation evaluates to:

I am not certain how to code that as a function of ‘T’ since it is not on a function of ‘T’ (or ‘t’), and that is required to evaluate γ.

I decided to let the Symbolic Math Toolbox have a go at solving for ‘P2’ and ‘T2’ since this system is significantly different than the first system (and the symbolic use here is a one-off derivation and does not involve iterative calculations).

This code:

syms gamma P P1 P2 Cp(t) R V V1 V2 T T1 T2
assume(P1 > 0)
assume(P2 > 0)
assume(T1 > 0)
assume(T2 > 0)
t = T/1000;
gamma = Cp(t)/(Cp(t)-R);
Eq1 = int(1/V, V, V1, V2) == -int((1/gamma)*(1/P), P, P1, P2);
Eq2 = int(1/P, P, P1, P2) == int((gamma/(gamma-1))*(1/T), T, T1, T2);
[P2S,T2S,Prms,Cnds] = solve([Eq1,Eq2],[P2,T2], 'ReturnConditions',1)

produced:

P2S =
z
 
T2S =
z1
 
Prms =
[ z, z1]
 
Cnds =
log(z) - log(P1) - int(Cp(T/1000)/(T*(R - Cp(T/1000))*(Cp(T/1000)/(R - Cp(T/1000)) + 1)), T, T1, z1) == 0 & 0 < z & 0 < z1 & piecewise(V1 <= 0 & 0 <= V2, int(1/V, V, V1, V2) - log(P1/z) + (R*log(P1/z))/Cp(T/1000), 0 < V1 | V2 < 0, log(V2) - log(V1) - log(P1/z) + (R*log(P1/z))/Cp(T/1000)) == 0

So there does not appear to be a symbolic solution for it (so I feel vindicated), and I am not certain how to code a numerical solution to it for the reasons I already mentioned.

So I am not certain where to go with this. I would appreciate any clarifications on the functional relationship of ‘P’ and ‘T’ in the first equation, since that appears to be important in expressing it as a function of ‘t’ that could be integrated numerically.

[Deema Khunda]

do you mean the expeted trend between 'T' and 'P'? they are proportional , as Pressure increases due to increase in volume, the temperature also increase, as the T temperature increases, gamma increases as its a function of temperature,

can we start the solutions as assuming an initial value of gamma calculated at T1=300, after obtaining this gamma value , solve the equation :

then after obtaiing a value for P2, apply this value of P2 into the following equation of P and T

After obtaining T2, use this input of T2 to calculate gamma at T2, then calculate another value of P2m then another T22 based on the new gamma, in this case we have an initial valye of gamma to start the solution with. so this will be implemented in the following loop :

[Star Strider]

‘do you mean the expeted trend between 'T' and 'P'? they are proportional , as Pressure increases due to increase in volume, the temperature also increase, as the T temperature increases, gamma increases as its a function of temperature’

I understand that. However they need to be related as mathematical functions of ‘T’. That mathematical relationship appears to be missing here. It is necessary in order to calculate γ, at least as described in these equations.

[Deema Khunda]

I do not have any other mathmatical relationship between T and P other than eq 2 the integral , I can give you the following relationship used for solving for P2 when gamma is constant

P2 = P1*(T2/T2)^(gama/(gamma-1))

but you are right in the fact that it doesnt seemt to have a symbolic solution, therefor other means of integration should be used.

would it be easier if Cp was just coded as function of T, bypassing the use of t , so the relation of Cp is :

Cp = @(T) A+B*T+C*((T/1000)^2)+D*((T/1000)^3)+E/((T/1000)^2);

[Torsten]

Beforehand, you could make a table of P2-T2 values from the second equation.

Then you can solve the two equations as you did before by interpolating gamma in the first equation as a function of P.

But I'd be surprised if the results of the two sets of equations were different because for constant gamma, they are equivalent.

[Star Strider]

Should this:

P2 = P1*(T2/T2)^(gamma/(gamma-1))

be ‘T1/T2’ or ‘T2/T1’? I doubt that would work anyway, because ‘gamma’ needs to be a function of ‘T’ for the rest of it to work.

I was an undergraduate Chemistry major (long ago), so encountered the gas laws and the approximations to deal with non-ideal gases, however I admit to being lost with this latest set of equations. The problem is that ‘gamma’ is a function only of ‘T’, and there is no existing relation that would work in the first equation, making it integrable as a function of ‘P’. In the absence of such a relation, I doubt that it is currently possible to procede further with this.

[Torsten]

Isn't it true that gamma(P) in the first equation has to be evaluated at the value of T which satisfies the second equation with P2 being replaced by P and T2 being replaced by T ?

Answer 2

[Ameer Hamza]

This solves the equations using symbolic equations

syms V P T P2 T2 
A= -0.703029;
B= 108.4773;
C= -42.52157;
D= 5.862788;
E= 0.678565;
P1= 2.7;
T1= 300;
V2= 30;
R=8;

t = T/1000;
Cp= A+B*t+C*(t^2)+D*(t^3)+E/(t^2);
gamma = Cp/(Cp-R);
V_val = 230:-0.5:V2;
P2_sol = zeros(size(V_val));
T2_sol = zeros(size(V_val));
gamma_vec = zeros(size(V_val));
for i=1:numel(V_val)
    V1 = V_val(i);
    
    eq1_lhs = int(1/P, P1, P2);
    eq1_rhs = int(gamma/(gamma-1)/T, T1, T2);
    eq2_lhs = int(1/V, V1, V2);
    eq2_rhs = -int(1/(gamma-1)/T, T1, T2);
    sol = vpasolve([eq1_lhs==eq1_rhs, eq2_lhs==eq2_rhs], [P2 T2]);
    
    P2_sol(i) = sol.P2; % solution for P2
    T2_sol(i) = sol.T2; % solution for T2
    gamma_vec(i) = subs(gamma, sol.T2);
end

%%
T = table(V_val', real(P2_sol)', T2_sol', gamma_vec', 'VariableNames', {'V1', 'P2', 'T2', 'gamma'});

Since this is using the symbolic toolbox, so the speed of execution can be slow. It will take a few minutes to finish.

[Note] As you mentioned in comment to Star Strider's comment, the volume is decreasing, but remember we start with V1 = 230, and end at V2 = 30, so in that case, we will maximum change in volume, and hence the maximum change in temperature and pressure. Now suppose we start at V1 = 150 and end at V1 = 30, the difference in volume is small, and therefore the change in temperature and pressure will also be small. I hope this clearify the confusion about decreasing values of P2 and T2 when we decrease V1.

[Torsten]

eq1_rhs = int(gamma/(gamma-1)/T,T1,T2)

eq2_rhs = int(-1/(gamma-1)/T,T1,T2)

[Ameer Hamza]

Torsten, thanks for pointing out. I overlooked it.

[Deema Khunda]

Thanks so much Ammer Hamza, can I ask how can I modify the answer to give me an output vector of solutions for V increasing incrementally from V1 to V2 in steps of 0.5?

Thank you

[Ameer Hamza]

Do you mean increase from V2 to V1 in steps of 0.5? V1 is larger, and V2 is smaller.

[Deema Khunda]

no I mean decrease from V1 to V2, can we implement it for decrease in steps from V1 to V2 ?

[Ameer Hamza]

So, if i understand correctly, you want to keep V2 constant, and vary the value of V1 from 230 to 30 in steps of 0.5?

[Deema Khunda]

yes, final V2 is fixed, I want to see a table of values as V1 decreases ( varies) incrementally from 230 to 30 in steps of 0.5

[Ameer Hamza]

Check the updated answer. Also refer to the note at the end for confusion about the values of P2 and T2.

[Deema Khunda]

Hey ,

I got an error in line 15 as

Unrecognized function or variable 'V_val'.

Error in combustion_trial (line 15)

P2_sol = zeros(size(V_val));

I think line line 15 V_val needs to be updated to V1_val as its defined earlier in the code as V1_val , I have updated it to the following :

syms V P T P2 T2 
A= -0.703029;
B= 108.4773;
C= -42.52157;
D= 5.862788;
E= 0.678565;
P1= 2.7;
T1= 300;
V2= 30;
R=8;
t = T/1000;
Cp= A+B*t+C*(t^2)+D*(t^3)+E/(t^2);
gamma = Cp/(Cp-R);
V_val = 230:-0.5:V2;
P2_sol = zeros(size(V_val));
T2_sol = zeros(size(V_val));
for i=1:numel(V_val)
    V1 = V_val(i);
    
    eq1_lhs = int(1/P, P1, P2);
    eq1_rhs = int(gamma/(gamma-1)/T, T1, T2);
    eq2_lhs = int(1/V, V1, V2);
    eq2_rhs = -int(1/(gamma-1)/T, T1, T2);
    sol = vpasolve([eq1_lhs==eq1_rhs, eq2_lhs==eq2_rhs], [P2 T2]);
    
    P2_sol(i) = sol.P2; % solution for P2
    T2_sol(i) = sol.T2; % solution for T2
end
%%
T = table(V_val', P2_sol', T2_sol', 'VariableNames', {'V1', 'P2', 'T2'});

Can I ask if I could add gamma to the table to see how its changing with temperature?

Thank you

[Ameer Hamza]

Deema, thats correct. That was a mistake in my code. Please check the updated answer. The value of gamma at T=T2 is also added to the table.