Can't solve the equations
1 view (last 30 days)
Show older comments
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 Comments
Answers (3)
Torsten
on 14 Oct 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 Comments
Torsten
on 14 Oct 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
on 14 Oct 2023
Edited: Walter Roberson
on 15 Oct 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 Comments
Walter Roberson
on 15 Oct 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
on 15 Oct 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
on 15 Oct 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 Comments
See Also
Categories
Find more on Linear Algebra in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!