Solving non-linear equations (0)

Question

Redwood on 18 Aug 2018

0
Link

Direct link to this question

https://nl.mathworks.com/matlabcentral/answers/415381-solving-non-linear-equations-0

Answered: Alex Sha on 31 Jul 2019

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 Comments
Show 1 older commentHide 1 older comment

Walter Roberson on 18 Aug 2018

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 on 18 Aug 2018

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

Sign in to comment.

Sign in to answer this question.

Answer 1

Redwood on 19 Aug 2018

0
Link

Direct link to this answer

https://nl.mathworks.com/matlabcentral/answers/415381-solving-non-linear-equations-0#answer_333293

Edited: Walter Roberson on 19 Aug 2018

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 Comment
Show -1 older commentsHide -1 older comments

Walter Roberson on 19 Aug 2018

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)) )

Sign in to comment.

Answer 2

Redwood on 19 Aug 2018

0
Link

Direct link to this answer

https://nl.mathworks.com/matlabcentral/answers/415381-solving-non-linear-equations-0#answer_333306

Edited: Walter Roberson on 19 Aug 2018

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 Comments
Show NoneHide None

Walter Roberson on 19 Aug 2018

You accidentally wrote

sol = solve( eq1,eq2, eq3. eq4, eq5, eq6 );

instead of

sol = solve( eq1,eq2, eq3, eq4, eq5, eq6 );

Redwood on 20 Aug 2018

Dear Walter,

Thank you so much!!!

Sincerely yours,

J1

Sign in to comment.

Answer 3

Redwood on 19 Aug 2018

0
Link

Direct link to this answer

https://nl.mathworks.com/matlabcentral/answers/415381-solving-non-linear-equations-0#answer_333307

Edited: Walter Roberson on 19 Aug 2018

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

0 Comments
Show -2 older commentsHide -2 older comments

Sign in to comment.

Answer 4

Redwood on 20 Aug 2018

0
Link

Direct link to this answer

https://nl.mathworks.com/matlabcentral/answers/415381-solving-non-linear-equations-0#answer_333420

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')

0 Comments
Show -2 older commentsHide -2 older comments

Sign in to comment.

Answer 5

Alex Sha on 31 Jul 2019

0
Link

Direct link to this answer

https://nl.mathworks.com/matlabcentral/answers/415381-solving-non-linear-equations-0#answer_385581

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

0 Comments
Show -2 older commentsHide -2 older comments

Sign in to comment.

Solving non-linear equations (0)

3 Comments
Show 1 older commentHide 1 older comment

Answers (5)

1 Comment
Show -1 older commentsHide -1 older comments

2 Comments
Show NoneHide None

0 Comments
Show -2 older commentsHide -2 older comments

0 Comments
Show -2 older commentsHide -2 older comments

0 Comments
Show -2 older commentsHide -2 older comments

See Also

Categories

Tags

Community Treasure Hunt

Solving non-linear equations (0)

3 Comments Show 1 older commentHide 1 older comment

Answers (5)

1 Comment Show -1 older commentsHide -1 older comments

2 Comments Show NoneHide None

0 Comments Show -2 older commentsHide -2 older comments

0 Comments Show -2 older commentsHide -2 older comments

0 Comments Show -2 older commentsHide -2 older comments

See Also

Categories

Tags

Community Treasure Hunt

3 Comments
Show 1 older commentHide 1 older comment

1 Comment
Show -1 older commentsHide -1 older comments

2 Comments
Show NoneHide None

0 Comments
Show -2 older commentsHide -2 older comments

0 Comments
Show -2 older commentsHide -2 older comments

0 Comments
Show -2 older commentsHide -2 older comments