Function handle with integrals of multiple equations?

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

14 个评论
显示 12更早的评论隐藏 12更早的评论

Star Strider 2020-4-4

在 MATLAB Online 中打开

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 2020-4-6

在 MATLAB Online 中打开

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 2020-4-6

在 MATLAB Online 中打开

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.

Star Strider 2020-4-6

在 MATLAB Online 中打开

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 2020-4-7

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 2020-4-4

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/515337-function-handle-with-integrals-of-multiple-equations#answer_423981

编辑：Ameer Hamza 2020-4-5

在 MATLAB Online 中打开

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.

10 个评论
显示 8更早的评论隐藏 8更早的评论

Deema Khunda 2020-4-5

在 MATLAB Online 中打开

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 2020-4-5

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.

请先登录，再进行评论。

Function handle with integrals of multiple equations?

6 个评论
显示 4更早的评论隐藏 4更早的评论

采纳的回答

14 个评论
显示 12更早的评论隐藏 12更早的评论

更多回答（1 个）

10 个评论
显示 8更早的评论隐藏 8更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Function handle with integrals of multiple equations?

6 个评论 显示 4更早的评论隐藏 4更早的评论

采纳的回答

14 个评论 显示 12更早的评论隐藏 12更早的评论

更多回答（1 个）

10 个评论 显示 8更早的评论隐藏 8更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

WeChat

6 个评论
显示 4更早的评论隐藏 4更早的评论

14 个评论
显示 12更早的评论隐藏 12更早的评论

10 个评论
显示 8更早的评论隐藏 8更早的评论