Can't solve the equations
1 vue (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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 commentaires
Réponses (3)
Torsten
le 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 commentaires
Torsten
le 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
le 14 Oct 2023
Modifié(e) : Walter Roberson
le 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 commentaires
Walter Roberson
le 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
le 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
le 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 commentaires
Voir également
Catégories
En savoir plus sur Linear Algebra dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!