Can't solve the equations
1 visualización (últimos 30 días)
Mostrar comentarios más antiguos
I want to solve for pb and pg and then make a chart of pb for different sb values. However, I am unable to solve these equations. Could you please help me?
syms pb pg
ds=0.25;
dm=0.2;
k=1;
sg=1;
m=0.01;
r=0.06;
f=0.075;
ms=0.2;
lg=1;
lb=1;
a=0.25;
b=0.25;
ab=((-ds+dm)*pb⁄2*k);
ag=((-ds+dm)*pg⁄2*k);
nb=((1⁄((r+f+sb+ab*ds-(m+ab*dm-2*(ms+(((k*ab^2))⁄((2*pb))))+lg*((pg⁄pb)-1))-ms)))^((1⁄((1-a-b)))));
ng=((1⁄((r+f+sg+ag*ds-(m+ag*dm-2*(ms+(((k*ag^2))⁄((2*pg))))+lb*((pb⁄pg)-1))-ms)))^((1⁄((1-a-b)))));
pbb=pb-nb;
pgg=pg-ng;
solve({pbb==0,pgg==0},{pb,pg})
0 comentarios
Respuestas (3)
Torsten
el 14 de Oct. de 2023
You must tell us which solution you want.
syms pb pg
sb = 1;
ds=0.25;
dm=0.2;
k=1;
sg=1;
m=0.01;
r=0.06;
f=0.075;
ms=0.2;
lg=1;
lb=1;
a=0.25;
b=0.25;
ab=((-ds+dm)*pb/2*k);
ag=((-ds+dm)*pg/2*k);
nb=((1/((r+f+sb+ab*ds-(m+ab*dm-2*(ms+(((k*ab^2))/((2*pb))))+lg*((pg/pb)-1))-ms)))^((1/((1-a-b)))));
ng=((1/((r+f+sg+ag*ds-(m+ag*dm-2*(ms+(((k*ag^2))/((2*pg))))+lb*((pb/pg)-1))-ms)))^((1/((1-a-b)))));
pbb=pb-nb;
pgg=pg-ng;
sol = vpasolve([pbb==0,pgg==0],[pb,pg])
sol.pb
sol.pg
3 comentarios
Torsten
el 14 de Oct. de 2023
You didn't answer the question. For sb = 1, MATLAB returns 21 solutions for pb and pg. You must have a criterion to sort out the solution of the 21 that you need.
Walter Roberson
el 14 de Oct. de 2023
Editada: Walter Roberson
el 15 de Oct. de 2023
You had a number of places where you were using the character ⁄ which is the "fraction slash" https://www.compart.com/en/unicode/U+2044
Also, it is not a good idea to use floating point numbers with solve(). Floating point numbers are in-exact, but solve() is designed to look for indefinitely-precise solutions. For example is your f intended to be exactly ![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1511834/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1511834/image.png)
syms sb %undefined in original
Q = @(v) sym(v);
syms pb pg
ds = Q(0.25);
dm = Q(0.2);
k = Q(1);
sg = Q(1);
m = Q(0.01);
r = Q(0.06);
f = Q(0.075);
ms = Q(0.2);
lg = Q(1);
lb = Q(1);
a = Q(0.25);
b = Q(0.25);
ab=((-ds+dm)*pb/2*k);
ag=((-ds+dm)*pg/2*k);
nb=((1/((r+f+sb+ab*ds-(m+ab*dm-2*(ms+(((k*ab^2))/((2*pb))))+lg*((pg/pb)-1))-ms)))^((1/((1-a-b)))));
ng=((1/((r+f+sg+ag*ds-(m+ag*dm-2*(ms+(((k*ag^2))/((2*pg))))+lb*((pb/pg)-1))-ms)))^((1/((1-a-b)))));
pbb=pb-nb;
pgg=pg-ng;
sol = solve([pbb==0,pgg==0],[pb,pg])
With this symbolic variable for sb (because you did not define sb in the original) then MATLAB is not able to find a solution.
If sb is given specific numeric values, MATLAB is able to solve the equations numerically (but not symbolically)
2 comentarios
Walter Roberson
el 15 de Oct. de 2023
Even if you restrict to positives, you are dealing with polynomials of degree 12 or 15 that have variable numbers of real roots depending on the value of sb. Why should you choose one value over another?
syms sb %undefined in original
Q = @(v) sym(v);
syms pb pg positive
ds = Q(0.25);
dm = Q(0.2);
k = Q(1);
sg = Q(1);
m = Q(0.01);
r = Q(0.06);
f = Q(0.075);
ms = Q(0.2);
lg = Q(1);
lb = Q(1);
a = Q(0.25);
b = Q(0.25);
ab=((-ds+dm)*pb/2*k);
ag=((-ds+dm)*pg/2*k);
nb=((1/((r+f+sb+ab*ds-(m+ab*dm-2*(ms+(((k*ab^2))/((2*pb))))+lg*((pg/pb)-1))-ms)))^((1/((1-a-b)))));
ng=((1/((r+f+sg+ag*ds-(m+ag*dm-2*(ms+(((k*ag^2))/((2*pg))))+lb*((pb/pg)-1))-ms)))^((1/((1-a-b)))));
pbb=pb-nb;
pgg=pg-ng;
sb_vals = 1:0.5:10;
num_sb = length(sb_vals);
for sb_idx = num_sb : -1 : 1
results(sb_idx) = solve(subs([pbb==0,pgg==0], sb, sb_vals(sb_idx)),[pb,pg],'maxdegree',4,'real',true);
end
results(1).pb
results(end).pb
results(1).pg
results(end).pg
Walter Roberson
el 15 de Oct. de 2023
format long g
syms sb %undefined in original
Q = @(v) sym(v);
syms pb pg positive
ds = Q(0.25);
dm = Q(0.2);
k = Q(1);
sg = Q(1);
m = Q(0.01);
r = Q(0.06);
f = Q(0.075);
ms = Q(0.2);
lg = Q(1);
lb = Q(1);
a = Q(0.25);
b = Q(0.25);
ab=((-ds+dm)*pb/2*k);
ag=((-ds+dm)*pg/2*k);
nb=((1/((r+f+sb+ab*ds-(m+ab*dm-2*(ms+(((k*ab^2))/((2*pb))))+lg*((pg/pb)-1))-ms)))^((1/((1-a-b)))));
ng=((1/((r+f+sg+ag*ds-(m+ag*dm-2*(ms+(((k*ag^2))/((2*pg))))+lb*((pb/pg)-1))-ms)))^((1/((1-a-b)))));
pbb=pb-nb;
pgg=pg-ng;
sb_vals = linspace(0,10,20);
num_sb = length(sb_vals);
ax1 = subplot(2,1,1);
ax2 = subplot(2,1,2);
for sb_idx = num_sb : -1 : 1
this_sb = sb_vals(sb_idx);
results(sb_idx) = solve(subs([pbb==0,pgg==0], sb, this_sb),[pb,pg],'maxdegree',4,'real',true);
scatter(ax1, this_sb, double(results(sb_idx).pb), []);
hold(ax1, 'on')
scatter(ax2, this_sb, double(results(sb_idx).pg), []);
hold(ax2, 'on');
end
title(ax1, 'pbb');
title(ax2, 'pgg');
hold(ax1, 'off');
hold(ax2, 'off');
pb_min = arrayfun(@(S) min(double(S.pb)), results);
pb_max = arrayfun(@(S) max(double(S.pb)), results);
pg_min = arrayfun(@(S) min(double(S.pg)), results);
pg_max = arrayfun(@(S) max(double(S.pg)), results);
[pb_min; pb_max]
[pg_min; pg_max]
Sam Chak
el 15 de Oct. de 2023
Hi @Della
I assume you want to plot something similar to this. However, please note that the solution set depends on the initial guess values. It's not something you can randomly select from the solution pool, such as pb = 0.204 and pg = 0.569. These values should match, for example, like {pb = 0.5699, pg = 0.5699} or {pb = 0.02283, pg = 0.20418}.
sb = 1:0.5:10;
options = optimoptions('fsolve', 'Display', 'none');
for j = 1:length(sb)
x0 = [0.2, 0.6];
x = fsolve(@(x) DellaFcn(x, sb(j)), x0, options);
pb = x(1);
plot(sb(j), pb, 'bo'), hold on, grid on
end
xlabel('sb'), ylabel('pb')
title('Solution depends on the initial guess x0')
function F = DellaFcn(x, param)
% definition
pb = x(1);
pg = x(2);
% parameter
sb = param;
% constants
ds=0.25;
dm=0.2;
k=1;
sg=1;
m=0.01;
r=0.06;
f=0.075;
ms=0.2;
lg=1;
lb=1;
a=0.25;
b=0.25;
% equations
ab=((-ds+dm)*pb/2*k);
ag=((-ds+dm)*pg/2*k);
nb=((1/((r+f+sb+ab*ds-(m+ab*dm-2*(ms+(((k*ab^2))/((2*pb))))+lg*((pg/pb)-1))-ms)))^((1/((1-a-b)))));
ng=((1/((r+f+sg+ag*ds-(m+ag*dm-2*(ms+(((k*ag^2))/((2*pg))))+lb*((pb/pg)-1))-ms)))^((1/((1-a-b)))));
pbb=pb-nb;
pgg=pg-ng;
% the system
F(1) = pbb;
F(2) = pgg;
end
0 comentarios
Ver también
Categorías
Más información sobre Linear Algebra en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!