Can't solve the equations
1 visualizzazione (ultimi 30 giorni)
Mostra commenti meno recenti
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 Commenti
Risposte (3)
Torsten
il 14 Ott 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 Commenti
Torsten
il 14 Ott 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
il 14 Ott 2023
Modificato: Walter Roberson
il 15 Ott 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 Commenti
Walter Roberson
il 15 Ott 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
il 15 Ott 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
il 15 Ott 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 Commenti
Vedere anche
Categorie
Scopri di più su Linear Algebra in Help Center e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!