Can't solve the equations

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})

Answer 1

[Torsten]

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 = struct with fields:
    pb: [22×1 sym]
    pg: [22×1 sym]

sol.pb

ans = 

sol.pg

ans = 

[Della]

Thanks. There is a range for sb. I want to create a chart for different values of sb.

[Torsten]

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.

[Della]

for sb=1,pb and pg can be pb=0.204, pg=0.569.

Answer 2

[Walter Roberson]

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

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)

[Walter Roberson]

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

ans = 

results(end).pb

ans = 

results(1).pg

ans = 

results(end).pg

ans = 

[Walter Roberson]

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]

ans = 2×20

1.0e+00 * 0.0291447919677483 0.0253732170108281 0.0225845056045588 0.0204201246960775 0.0186813852321928 0.017247971886532 0.0160421733371156 0.0150112624587187 0.0141180439986893 0.0133354371911016 0.01264320221554 0.0120258743986351 0.011471412414968 0.0109702860280629 0.010514843909762 0.0100988653213731 0.00971723567039779 0.00936570745894238 0.00858862305891743 0.00780732288795605 1335.14028993321 1974.97009576186 2145.78729108529 1511.14961528718 1182.27703715526 986.799195296022 852.743665892573 753.481228112936 676.341488793207 614.342423177441 563.251375567278 520.323911377975 483.689016894767 452.020425979667 424.348168910396 399.944402661918 378.251092930994 358.832487943541 341.342901911571 325.50429056599

[pg_min; pg_max]

ans = 2×20

1.0e+00 * 0.169444542378108 0.0510938957530228 0.0210878605540615 0.010322326167264 0.00564394458555274 0.00333999026206204 0.002098852177844 0.00138297117091768 0.000947069253688153 0.00066965071291072 0.000486465283261733 0.000361663165518009 0.000274321261451029 0.000211750975666033 0.000165998240266018 0.000131930430457683 0.000106149867321332 8.63556982317282e-05 7.095781435412e-05 5.88369829950541e-05 702.692953408748 1291.87212611237 2300.74505672041 2954.42302230835 3163.9182232134 3275.11623634661 3347.01826529229 3398.28760833292 3437.0798587401 3467.63772085638 3492.42591698404 3512.99039458511 3530.35729751647 3545.23813004112 3558.14368294439 3569.45121442342 3579.44603036459 3588.34826647822 3596.33072574168 3603.53110195361

Answer 3

[Sam Chak]

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