Solving non-linear equations (0)

1 次查看（过去 30 天）
Redwood
Redwood 2018-8-18
回答： Alex Sha 2019-7-31
Dear all,
I solved non-linear equations using Matlab, but all values are close to 0.
I think this solution is not incorrect, and I would like to know my coding is correct or not.
In addition, I would like to get the value of Me1. Please let me know how to get this value.
Thank you very much in advance.
Sincerely yours,
J1
function F = costate2(x)
t = 0;
tlast = 30;
phi = 3.14;
ep = 0.1;
L1 = 157.5;
L2 = 799.5;
GHGs = 2954;
delta1 = 0.005;
gam1 = 0.10957;
eta11 = 0.40022;
eta21 = 0.99931;
ag1 = 20.49921;
alpha1 = 15.93413;
rhoq1 = 0.00710;
beta11 = 0.35;
beta21 = 0.63;
beta31 = 0.02;
g1 = 2.6005;
chi1 = 0.01;
xi1 = 0.02;
rhom1 = 0.0041773;
delta2 = 0.00581;
gam2 = 0.14990;
eta12 = 0.08;
eta22 = 0.59998;
ag2 = 21.49967;
alpha2 = 3.25702;
rhoq2 = 0.005;
beta12 = 0.60;
beta22 = 0.38;
beta32 = 0.02;
g2 = 2.8643;
chi2 = 0.01;
xi2 = 0.02;
rhom2 = -0.00800;
d = 0.05;
psi = 0.1;
K1 = 3486;
K2 = 1951;
Kg1 = 3.49;
Kg2 = 1.95;
Pw = 0.1;
Pg = 0.05;
GHG = 3162;
Q1 = 1496;
Q2 = 604;
E1 = 2216;
E2 = 2417;
a1 = 0.1;
a2 = 0.5;
b1 = 0.5;
b2 = 0.5;
Gbar = 12000;
P1 = 309;
P2 = 1338;
gam3 = 0.0492288;
gam4 = 0.0243931;
d1 = 1.2521405;
delta3 = 0.0250461;
lamda1 = x(1);
qg1 = x(2);
qp1 = x(3);
qe1 = x(4);
qc1 = x(5);
qk1 = x(6);
Gbar1 = (P1/(P1+P2))*Gbar;
Gbar2 = (P2/(P1+P2))*Gbar;
M1 = (g1*exp(-rhom1*t)-chi1*(Kg1/(K1+Kg1))^xi1)*E1;
W1 = (Pw-E1-Pg*(P1./(P1+P2))*Gbar)^2+4*beta31*Pw*Q1*M1*Pg;
J1 = (Pg*(P1./(P1+P2))*Gbar-Pw*E1)^2+4*beta31*Pw*Q1*M1*Pg;
Pe1 = ((Pw*E1-Pg*Gbar1) + sqrt(W1))/(2*E1);
Me1 = ((Pg*(P1./(P1+P2))*Gbar-Pw*E1)+sqrt(J1))/(2*Pg);
Me2 = (g2*exp(-rhom2*t)-chi2*(Kg2./(K2+Kg2))^xi2)*beta32*Q2;
Y1 = E1^2*Pe1^2+beta31*Pw*Q1*M1*Pg;
QQe1 = alpha1*exp(rhoq1*t)*(K1+Kg1)^beta11*L1^beta21*E1^beta31;
Op = Me1 - (P1./(P1+P2))*Gbar;
PP = P1*P2;
F(1) = x(1)*( phi/(2*tlast) )
-x(3)*( psi*(P1+P2)*Pg/(sqrt(PP)*Gbar) )*( (chi1*xi1*beta31*Pw*Q1*E1^2*Pe1*Kg1^xi1)/(Y1*(K1+Kg1)^(xi1+1)) )
-x(4)*chi1*xi1*E1*(Kg1^xi1)/(K1+Kg1)^(xi1+1)
-x(5)*a1*beta11*QQe1/(K1+Kg1)
-x(6)*b1*(beta31*Pw*Q1*Pg/(E1*Pe1^2))*( (chi1*xi1*beta31*Pw*Q1*E1^2*Pe1*Kg1^xi1)/(Y1*(K1+Kg1)^(xi1+1)) );
F(2) = x(2)*( phi/(2*tlast) + 1/ag1 )
-x(2).^2/(x(1)*2*ag1)
-x(1)*( 1/(2*ag1) - Pg*( (chi1*xi1*beta31*Pw*Q1*E1^2*Pe1*Kg1^xi1)*(1/(K1+Kg1)-1/Kg1)/(Y1*(K1+Kg1)^xi1) ) )
-x(3)*psi*(P1+P2)*Pg*(1/(sqrt(PP)*Gbar))* ( (chi1*xi1*beta31*Pw*Q1*E1^2*Pe1*Kg1^xi1)*(1/(K1+Kg1)-1/Kg1)/(Y1*(K1+Kg1)^xi1) )
-x(4)*chi1*xi1*E1*((Kg1/(K1+Kg1))^xi1)*(1/(K1+Kg1)-1/Kg1)
-x(5)*a1*beta11*QQe1/(K1+Kg1)
+x(6)*b1*beta31*Pw*Q1*Pg*(1/(E1*Pe1^2))*( (chi1*xi1*beta31*Pw*Q1*E1^2*Pe1*Kg1^xi1)*(1/(K1+Kg1)-1/Kg1)/(Y1*(K1+Kg1)^xi1) )
;
F(3) = x(3)*( phi/(2*tlast) )
+x(1)*( (Me1-Gbar1) - Pg*beta31*Pw*Q1*M1*(Me1-Gbar1)/Y1 )
+x(6)*b1*beta31*Pw*Q1*M1*(Me1-Gbar1)/Y1;
F(4) = x(4)*( phi/(2*tlast) )
+x(1).^((eta11*gam1+gam1)/(eta11*gam1+gam1-1)) * eta21*(ep/eta11)^(eta11*gam1/(eta11+gam1-1)) * GHGs^(-eta21*gam1/(eta11*gam1+gam1-1)) * GHG^(((eta21-eta11)*gam1-gam1+1)/(eta11*gam1+gam1-1));
F(5) = x(5)*( phi/(2*tlast) )
-x(1)*(1-Pg*beta31*Pw+M1*E1*Pe1/Y1)
-x(3)*psi*(P1+P2)*Pg*(1/(sqrt(PP)*Gbar))*beta31*Pw+M1*E1*Pe1/Y1
-x(6)*b1*beta31*Pw+M1*E1*Pe1/Y1;
F(6) = x(6)*( phi/(2*tlast) )
+x(1)*Pg*M1*(1/E1)*beta31*Pw*Q1*Pg*(Me1-Gbar1)/Y1
-x(3)*psi*(P1+P2)*Pg*(1/(sqrt(PP)*Gbar))*M1*(1/E1)*beta31*Pw*Q1*Pg*(Me1-Gbar1)/Y1
-x(4)*M1/E1
-x(5)*a1*beta31*QQe1/E1;
fun = @costate2;
x0 = [1,1,1,1,1,1];
x = fsolve(fun,x0)
  3 个评论
显示 1更早的评论
Walter Roberson
Walter Roberson 2018-8-18
You did not present the equations in symbolic form such as an image of the equations, so as far as we as onlookers are concerned, the equations are defined by your current code, and so must be correct. We would need some other representation of the equations in order to be able to check the code against the equations.
Star Strider
Star Strider 2018-8-18
Since ‘Me1’ does not appear to be a function of ‘x’, to display its value, put this line after the ‘Me1’ assignment:
fprintf(1, '\n\t\tMe1 = %23.15E\n', Me1)
in context:
Me1 = ((Pg*(P1./(P1+P2))*Gbar-Pw*E1)+sqrt(J1))/(2*Pg);
fprintf(1, '\n\t\tMe1 = %23.15E\n', Me1)
It will display the value of ‘Me1’ (several times) in the Command Window.
The result will be:
Me1 = 1.476154199966165E+02

请先登录，再进行评论。

请先登录，再回答此问题。

回答（5 个）

Redwood
Redwood 2018-8-19
编辑：Walter Roberson 2018-8-19
Dear all,
Thank you very much for your help.
I revised my code to solve non-linear equations.
I got some errors, and I would like to know how to fix these errors.
Thank you very much for your help.
Sincerely yours,
J1
t = 0;
tlast = 30;
phi = 3.14;
ep = 0.1;
L1 = 157.5;
L2 = 799.5;
GHGs = 2954;
delta1 = 0.005;
gam1 = 0.10957;
eta11 = 0.40022;
eta21 = 0.99931;
ag1 = 20.49921;
alpha1 = 15.93413;
rhoq1 = 0.00710;
beta11 = 0.35;
beta21 = 0.63;
beta31 = 0.02;
g1 = 2.6005;
chi1 = 0.01;
xi1 = 0.02;
rhom1 = 0.0041773;
delta2 = 0.00581;
gam2 = 0.14990;
eta12 = 0.08;
eta22 = 0.59998;
ag2 = 21.49967;
alpha2 = 3.25702;
rhoq2 = 0.005;
beta12 = 0.60;
beta22 = 0.38;
beta32 = 0.02;
g2 = 2.8643;
chi2 = 0.01;
xi2 = 0.02;
rhom2 = -0.00800;
d = 0.05;
psi = 0.1;
K1 = 3486;
K2 = 1951;
Kg1 = 3.49;
Kg2 = 1.95;
Pw = 0.1;
Pg = 0.05;
GHG = 3162;
Q1 = 1496;
Q2 = 604;
E1 = 2216;
E2 = 2417;
a1 = 0.1;
a2 = 0.5;
b1 = 0.5;
b2 = 0.5;
Gbar = 12000;
P1 = 309;
P2 = 1338;
gam3 = 0.0492288;
gam4 = 0.0243931;
d1 = 1.2521405;
delta3 = 0.0250461;
Gbar1 = (P1/(P1+P2))*Gbar;
Gbar2 = (P2/(P1+P2))*Gbar;
M1 = (g1*exp(-rhom1*t)-chi1*(Kg1/(K1+Kg1))^xi1)*E1;
W1 = (Pw-E1-Pg*(P1./(P1+P2))*Gbar)^2+4*beta31*Pw*Q1*M1*Pg;
J1 = (Pg*(P1./(P1+P2))*Gbar-Pw*E1)^2+4*beta31*Pw*Q1*M1*Pg;
Pe1 = ((Pw*E1-Pg*Gbar1) + sqrt(W1))/(2*E1);
Me1 = ((Pg*(P1./(P1+P2))*Gbar-Pw*E1)+sqrt(J1))/(2*Pg );
Me2 = (g2*exp(-rhom2*t)-chi2*(Kg2./(K2+Kg2))^xi2)*beta32*Q2;
Y1 = E1^2*Pe1^2+beta31*Pw*Q1*M1*Pg;
QQe1 = alpha1*exp(rhoq1*t)*(K1+Kg1)^beta11*L1^beta21*E1^beta31;
Op = Me1 - (P1./(P1+P2))*Gbar;
PP = P1*P2;
syms lamda1 qg1 qp1 qn1 qc1 qk1
eq1 = lamda1*( phi/(2*tlast) )
-qp1*( psi*(P1+P2)*Pg/(sqrt(PP)*Gbar) )*( (chi1*xi1*beta31*Pw*Q1*E1^2*Pe1*Kg1^xi1)/(Y1*(K1+Kg1)^(xi1+1)) )
-qn1*chi1*xi1*E1*(Kg1^xi1)/(K1+Kg1)^(xi1+1)
-qc1*a1*beta11*QQe1/(K1+Kg1)
-qk1*b1*(beta31*Pw*Q1*Pg/(E1*Pe1^2))*( (chi1*xi1*beta31*Pw*Q1*E1^2*Pe1*Kg1^xi1)/(Y1*(K1+Kg1)^(xi1+1)) );
eq2 = qg1*( phi/(2*tlast) + 1/ag1 )
-qg1^2/(lamda1*2*ag1)
-lamda1*( 1/(2*ag1) - Pg*( (chi1*xi1*beta31*Pw*Q1*E1^2*Pe1*Kg1^xi1)*(1/(K1+Kg1)-1/Kg1)/(Y1*(K1+Kg1)^xi1) ) )
-qp1*psi*(P1+P2)*Pg*(1/(sqrt(PP)*Gbar))* ( (chi1*xi1*beta31*Pw*Q1*E1^2*Pe1*Kg1^xi1)*(1/(K1+Kg1)-1/Kg1)/(Y1*(K1+Kg1)^xi1) )
-qn1*chi1*xi1*E1*((Kg1/(K1+Kg1))^xi1)*(1/(K1+Kg1)-1/Kg1)
-qc1*a1*beta11*QQe1/(K1+Kg1)
+qk1*b1*beta31*Pw*Q1*Pg*(1/(E1*Pe1^2))*( (chi1*xi1*beta31*Pw*Q1*E1^2*Pe1*Kg1^xi1)*(1/(K1+Kg1)-1/Kg1)/(Y1*(K1+Kg1)^xi1) );
eq3 = qp1*( phi/(2*tlast) )
+lamda1*( (Me1-Gbar1) - Pg*beta31*Pw*Q1*M1*(Me1-Gbar1)/Y1 )
+qk1*b1*beta31*Pw*Q1*M1*(Me1-Gbar1)/Y1;
eq4 = qn1*( phi/(2*tlast) )
+lamda1^((eta11*gam1+gam1)/(eta11*gam1+gam1-1)) * eta21*(ep/eta11)^(eta11*gam1/(eta11+gam1-1)) * GHGs^(-eta21*gam1/(eta11*gam1+gam1-1)) * GHG^(((eta21-eta11)*gam1-gam1+1)/(eta11*gam1+gam1-1));
eq5 = qc1*( phi/(2*tlast) )
-lamda1*(1-Pg*beta31*Pw+M1*E1*Pe1/Y1)
-qp1*psi*(P1+P2)*Pg*(1/(sqrt(PP)*Gbar))*beta31*Pw+M1*E1*Pe1/Y1
-qk1*b1*beta31*Pw+M1*E1*Pe1/Y1;
eq6 = qk1*( phi/(2*tlast) )
+lamda1*Pg*M1*(1/E1)*beta31*Pw*Q1*Pg*(Me1-Gbar1)/Y1
-qp1*psi*(P1+P2)*Pg*(1/(sqrt(PP)*Gbar))*M1*(1/E1)*beta31*Pw*Q1*Pg*(Me1-Gbar1)/Y1
-qn1*M1/E1
-qc1*a1*beta31*QQe1/E1;
sol = solve( eq1,eq2, eq3. eq4, eq5, eq6 );
sol. lamda1;
sol. qg1;
sol. qp1;
sol. qn1;
sol. qc1;
sol. qk1;
  1 个评论
Walter Roberson
Walter Roberson 2018-8-19
You are missing continuation marks for your eq* variables.
eq1 = lamda1*( phi/(2*tlast) )
-qp1*( psi*(P1+P2)*Pg/(sqrt(PP)*Gbar) )*( (chi1*xi1*beta31*Pw*Q1*E1^2*Pe1*Kg1^xi1)/(Y1*(K1+Kg1)^(xi1+1)) )
is not considered all one long expression: it is considered an assignment, followed by the calculation and display (but no assignment) of the expression on the second line.
You would need
eq1 = lamda1*( phi/(2*tlast) ) ...
-qp1*( psi*(P1+P2)*Pg/(sqrt(PP)*Gbar) )*( (chi1*xi1*beta31*Pw*Q1*E1^2*Pe1*Kg1^xi1)/(Y1*(K1+Kg1)^(xi1+1)) )

请先登录，再进行评论。

Redwood
Redwood 2018-8-19
编辑：Walter Roberson 2018-8-19
Dear Walter,
Thank you very much for your comment.
I revised my coding, and I still got some errors.
Please let me know how to fix this problem.
Thank you very much in advance.
Sincerely yours,
J1
t = 0;
tlast = 30;
phi = 3.14;
ep = 0.1;
L1 = 157.5;
L2 = 799.5;
GHGs = 2954;
delta1 = 0.005;
gam1 = 0.10957;
eta11 = 0.40022;
eta21 = 0.99931;
ag1 = 20.49921;
alpha1 = 15.93413;
rhoq1 = 0.00710;
beta11 = 0.35;
beta21 = 0.63;
beta31 = 0.02;
g1 = 2.6005;
chi1 = 0.01;
xi1 = 0.02;
rhom1 = 0.0041773;
delta2 = 0.00581;
gam2 = 0.14990;
eta12 = 0.08;
eta22 = 0.59998;
ag2 = 21.49967;
alpha2 = 3.25702;
rhoq2 = 0.005;
beta12 = 0.60;
beta22 = 0.38;
beta32 = 0.02;
g2 = 2.8643;
chi2 = 0.01;
xi2 = 0.02;
rhom2 = -0.00800;
d = 0.05;
psi = 0.1;
K1 = 3486;
K2 = 1951;
Kg1 = 3.49;
Kg2 = 1.95;
Pw = 0.1;
Pg = 0.05;
GHG = 3162;
Q1 = 1496;
Q2 = 604;
E1 = 2216;
E2 = 2417;
a1 = 0.1;
a2 = 0.5;
b1 = 0.5;
b2 = 0.5;
Gbar = 12000;
P1 = 309;
P2 = 1338;
gam3 = 0.0492288;
gam4 = 0.0243931;
d1 = 1.2521405;
delta3 = 0.0250461;
Gbar1 = (P1/(P1+P2))*Gbar;
Gbar2 = (P2/(P1+P2))*Gbar;
M1 = (g1*exp(-rhom1*t)-chi1*(Kg1/(K1+Kg1))^xi1)*E1;
W1 = (Pw-E1-Pg*(P1./(P1+P2))*Gbar)^2+4*beta31*Pw*Q1*M1*Pg;
J1 = (Pg*(P1./(P1+P2))*Gbar-Pw*E1)^2+4*beta31*Pw*Q1*M1*Pg;
Pe1 = ((Pw*E1-Pg*Gbar1) + sqrt(W1))/(2*E1);
Me1 = ((Pg*(P1./(P1+P2))*Gbar-Pw*E1)+sqrt(J1))/(2*Pg );
Me2 = (g2*exp(-rhom2*t)-chi2*(Kg2./(K2+Kg2))^xi2)*beta32*Q2;
Y1 = E1^2*Pe1^2+beta31*Pw*Q1*M1*Pg;
QQe1 = alpha1*exp(rhoq1*t)*(K1+Kg1)^beta11*L1^beta21*E1^beta31;
Op = Me1 - (P1./(P1+P2))*Gbar;
PP = P1*P2;
syms lamda1 qg1 qp1 qn1 qc1 qk1
eq1 = lamda1*( phi/(2*tlast) )...
-qp1*( psi*(P1+P2)*Pg/(sqrt(PP)*Gbar) )*( (chi1*xi1*beta31*Pw*Q1*E1^2*Pe1*Kg1^xi1)/(Y1*(K1+Kg1)^(xi1+1)) )
-qn1*chi1*xi1*E1*(Kg1^xi1)/(K1+Kg1)^(xi1+1)
-qc1*a1*beta11*QQe1/(K1+Kg1)
-qk1*b1*(beta31*Pw*Q1*Pg/(E1*Pe1^2))*( (chi1*xi1*beta31*Pw*Q1*E1^2*Pe1*Kg1^xi1)/(Y1*(K1+Kg1)^(xi1+1)) );
eq2 = qg1*( phi/(2*tlast) + 1/ag1 )...
-qg1^2/(lamda1*2*ag1) -lamda1*( 1/(2*ag1) - Pg*( (chi1*xi1*beta31*Pw*Q1*E1^2*Pe1*Kg1^xi1)*(1/(K1+Kg1)-1/Kg1)/(Y1*(K1+Kg1)^xi1) ) )
-qp1*psi*(P1+P2)*Pg*(1/(sqrt(PP)*Gbar))* ( (chi1*xi1*beta31*Pw*Q1*E1^2*Pe1*Kg1^xi1)*(1/(K1+Kg1)-1/Kg1)/(Y1*(K1+Kg1)^xi1) )
-qn1*chi1*xi1*E1*((Kg1/(K1+Kg1))^xi1)*(1/(K1+Kg1)-1/Kg1)
-qc1*a1*beta11*QQe1/(K1+Kg1)
+qk1*b1*beta31*Pw*Q1*Pg*(1/(E1*Pe1^2))*( (chi1*xi1*beta31*Pw*Q1*E1^2*Pe1*Kg1^xi1)*(1/(K1+Kg1)-1/Kg1)/(Y1*(K1+Kg1)^xi1) );
eq3 = qp1*( phi/(2*tlast) )...
+lamda1*( (Me1-Gbar1) - Pg*beta31*Pw*Q1*M1*(Me1-Gbar1)/Y1 )
+qk1*b1*beta31*Pw*Q1*M1*(Me1-Gbar1)/Y1;
eq4 = qn1*( phi/(2*tlast) )...
+lamda1^((eta11*gam1+gam1)/(eta11*gam1+gam1-1)) * eta21*(ep/eta11)^(eta11*gam1/(eta11+gam1-1)) * GHGs^(-eta21*gam1/(eta11*gam1+gam1-1)) * GHG^(((eta21-eta11)*gam1-gam1+1)/(eta11*gam1+gam1-1));
eq5 = qc1*( phi/(2*tlast) )...
-lamda1*(1-Pg*beta31*Pw+M1*E1*Pe1/Y1)
-qp1*psi*(P1+P2)*Pg*(1/(sqrt(PP)*Gbar))*beta31*Pw+M1*E1*Pe1/Y1
-qk1*b1*beta31*Pw+M1*E1*Pe1/Y1;
eq6 = qk1*( phi/(2*tlast) )...
+lamda1*Pg*M1*(1/E1)*beta31*Pw*Q1*Pg*(Me1-Gbar1)/Y1
-qp1*psi*(P1+P2)*Pg*(1/(sqrt(PP)*Gbar))*M1*(1/E1)*beta31*Pw*Q1*Pg*(Me1-Gbar1)/Y1
-qn1*M1/E1 -qc1*a1*beta31*QQe1/E1;
sol = solve( eq1,eq2, eq3. eq4, eq5, eq6 );
sol. lamda1;
sol. qg1;
sol. qp1;
sol. qn1;
sol. qc1;
sol. qk1;
  2 个评论
Walter Roberson
Walter Roberson 2018-8-19
You accidentally wrote
sol = solve( eq1,eq2, eq3. eq4, eq5, eq6 );
instead of
sol = solve( eq1,eq2, eq3, eq4, eq5, eq6 );
Redwood
Redwood 2018-8-20
Dear Walter,
Thank you so much!!!
Sincerely yours,
J1

请先登录，再进行评论。

Redwood
Redwood 2018-8-19
编辑：Walter Roberson 2018-8-19
Dear Walter,
I tried to put "..."
For example,
eq1 = lamda1*( phi/(2*tlast) )...
-qp1*( psi*(P1+P2)*Pg/(sqrt(PP)*Gbar) )*( (chi1*xi1*beta31*Pw*Q1*E1^2*Pe1*Kg1^xi1)/(Y1*(K1+Kg1)^(xi1+1)) )...
-qn1*chi1*xi1*E1*(Kg1^xi1)/(K1+Kg1)^(xi1+1)...
-qc1*a1*beta11*QQe1/(K1+Kg1)...
-qk1*b1*(beta31*Pw*Q1*Pg/(E1*Pe1^2))*( (chi1*xi1*beta31*Pw*Q1*E1^2*Pe1*Kg1^xi1)/(Y1*(K1+Kg1)^(xi1+1)) );
But I also got some errors. Please let me know how to fix this problem
Thank you very much in advance.
Sincerely yours,
J1
Redwood
Redwood 2018-8-20
Dear Walter,
I got some answers and I think my code has no errors.
But I cannot see the values of the unknown variables(lamda1, qg1, qp1, qn1, qc1, qk1), and I am not sure if I solver non-linear equations or not.
Please let me know how to get the values of the unknown variables(lamda1, qg1, qp1, qn1, qc1, qk1).
Thank you very much in advance.
Sincerely yours,
J1
t = 0; tlast = 30; phi = 3.14;
ep = 0.1; L1 = 157.5; L2 = 799.5; GHGs = 2954;
delta1 = 0.005; gam1 = 0.10957; eta11 = 0.40022; eta21 = 0.99931;
ag1 = 20.49921;
alpha1 = 15.93413; rhoq1 = 0.00710; beta11 = 0.35; beta21 = 0.63; beta31 = 0.02;
g1 = 2.6005; chi1 = 0.01; xi1 = 0.02; rhom1 = 0.0041773;
delta2 = 0.00581; gam2 = 0.14990; eta12 = 0.08; eta22 = 0.59998;
ag2 = 21.49967;
alpha2 = 3.25702; rhoq2 = 0.005; beta12 = 0.60; beta22 = 0.38; beta32 = 0.02;
g2 = 2.8643; chi2 = 0.01; xi2 = 0.02; rhom2 = -0.00800; d = 0.05; psi = 0.1;
K1 = 3486; K2 = 1951; Kg1 = 3.49; Kg2 = 1.95; Pw = 0.1; Pg = 0.05; GHG = 3162; Q1 = 1496; Q2 = 604; E1 = 2216; E2 = 2417;
a1 = 0.1; a2 = 0.5; b1 = 0.5; b2 = 0.5;
Gbar = 12; P1 = 309; P2 = 1338;
gam3 = 0.0492288; gam4 = 0.0243931; d1 = 1.2521405; delta3 = 0.0250461;
Gbar1 = (P1/(P1+P2))*Gbar; Gbar2 = (P2/(P1+P2))*Gbar; M1 = (g1*exp(-rhom1*t)-chi1*(Kg1/(K1+Kg1))^xi1)*E1; W1 = (Pw-E1-Pg*(P1./(P1+P2))*Gbar)^2+4*beta31*Pw*Q1*M1*Pg; J1 = (Pg*(P1./(P1+P2))*Gbar-Pw*E1)^2+4*beta31*Pw*Q1*M1*Pg; Pe1 = ((Pw*E1-Pg*Gbar1) + sqrt(W1))/(2*E1);
Me1 = ((Pg*(P1./(P1+P2))*Gbar-Pw*E1)+sqrt(J1))/(2*Pg ); Me2 = (g2*exp(-rhom2*t)-chi2*(Kg2./(K2+Kg2))^xi2)*beta32*Q2; Y1 = E1^2*Pe1^2+beta31*Pw*Q1*M1*Pg; QQe1 = alpha1*exp(rhoq1*t)*(K1+Kg1)^beta11*L1^beta21*E1^beta31; Op = Me1 - (P1./(P1+P2))*Gbar; PP = P1*P2;
syms lamda1 qg1 qp1 qn1 qc1 qk1
eq1 = lamda1*( -phi/(2*tlast) )... +qp1*( psi*(P1+P2)*Pg/(sqrt(PP)*Gbar) )*( (chi1*xi1*beta31*Pw*Q1*E1^2*Pe1*Kg1^xi1)/(Y1*(K1+Kg1)^(xi1+1)) ) +qn1*chi1*xi1*E1*(Kg1^xi1)/(K1+Kg1)^(xi1+1) +qc1*a1*beta11*QQe1/(K1+Kg1) -qk1*b1*(beta31*Pw*Q1*Pg/(E1*Pe1^2))*( (chi1*xi1*beta31*Pw*Q1*E1^2*Pe1*Kg1^xi1)/(Y1*(K1+Kg1)^(xi1+1)) );
eq2 = qg1*( -phi/(2*tlast) )... +lamda1*( 1/(2*ag1) - Pg*( (chi1*xi1*beta31*Pw*Q1*E1^2*Pe1*Kg1^xi1)*(1/(K1+Kg1)-1/Kg1)/(Y1*(K1+Kg1)^xi1) ) ) +qp1*psi*(P1+P2)*Pg*(1/(sqrt(PP)*Gbar))* ( (chi1*xi1*beta31*Pw*Q1*E1^2*Pe1*Kg1^xi1)*(1/(K1+Kg1)-1/Kg1)/(Y1*(K1+Kg1)^xi1) ) +qn1*chi1*xi1*E1*((Kg1/(K1+Kg1))^xi1)*(1/(K1+Kg1)-1/Kg1) +qc1*a1*beta11*QQe1/(K1+Kg1) -qk1*b1*beta31*Pw*Q1*Pg*(1/(E1*Pe1^2))*( (chi1*xi1*beta31*Pw*Q1*E1^2*Pe1*Kg1^xi1)*(1/(K1+Kg1)-1/Kg1)/(Y1*(K1+Kg1)^xi1) );
eq3 = qp1*( -phi/(2*tlast) )... -lamda1*( (Me1-Gbar1) - Pg*beta31*Pw*Q1*M1*(Me1-Gbar1)/Y1 ) -qk1*b1*beta31*Pw*Q1*M1*(Me1-Gbar1)/Y1;
eq4 = qn1*( phi/(2*tlast) )... +lamda1^((eta11*gam1+gam1)/(eta11*gam1+gam1-1)) * eta21*(ep/eta11)^(eta11*gam1/(gam1+gam1*eta11-1)) * GHGs^(-eta21*gam1/(eta11*gam1+gam1-1)) * GHG^(((eta21-eta11)*gam1-gam1+1)/(eta11*gam1+gam1-1));
eq5 = qc1*( -phi/(2*tlast) )... +lamda1*(1-Pg*beta31*Pw+M1*E1*Pe1/Y1) +qp1*psi*(P1+P2)*Pg*(1/(sqrt(PP)*Gbar))*beta31*Pw+M1*E1*Pe1/Y1 +qk1*b1*beta31*Pw+E1^2*Pe1/Y1;
eq6 = qk1*( phi/(2*tlast) )... +lamda1*Pg*M1*(1/E1)*beta31*Pw*Q1*Pg*(Me1-Gbar1)/Y1 -qp1*psi*(P1+P2)*Pg*(1/(sqrt(PP)*Gbar))*M1*(1/E1)*beta31*Pw*Q1*Pg*(Me1-Gbar1)/Y1 -qn1*M1/E1 -qc1*a1*beta31*QQe1/E1;
sol = solve( eq1,eq2, eq3, eq4, eq5, eq6 ); sol. lamda1; sol. qg1; sol. qp1; sol. qn1; sol. qc1; sol. qk1;
disp('lamda1, qg1, qp1, qn1, qc1, qk1')
Alex Sha
Alex Sha 2019-7-31
try numerical solution:
lamda1: 1.48312225679967E-24
qp1: 54185.6537859035
qn1: -133.888674663691
qc1: 0.188502621518643
qk1: -6630.74381075639
qg1: 282.666798246712
Me1=76.2725002858213

类别

在 Help Center 和 File Exchange 中查找有关 Mathematics 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by