Solving non-linear equations (0)
1 view (last 30 days)
Show older comments
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
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
Answers (5)
Redwood
on 19 Aug 2018
Edited: Walter Roberson
on 19 Aug 2018
1 Comment
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)) )
Redwood
on 19 Aug 2018
Edited: Walter Roberson
on 19 Aug 2018
2 Comments
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 );
Alex Sha
on 31 Jul 2019
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
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!