How to solve these non-linear equations?

syms a b c
i1 = 310;
i2 = 349.64;
i3 = 353;
c1 = 11.1984;
n1 = 0.5067;
c2 = 15.9867;
n2 = 0.4271;
c3 = 8.6028;
n3 = 0.2449;
mpl = 308.5;
eqn1 = 48.62236629 - (a + b + c + mpl)/mpl == 0;
eqn2 = i1*(1 - n1 * c1* a^(n1-1))/(1-c1*a^(n1-1)) *(1-c1*a^(n1-1)*((a+b+c+mpl)/(c1*a^n1 + b+c+mpl))) - i2 * (1-n2*c2*b^(n2-1)*((b+c+mpl)/(c2*b^n2+c+mpl))) == 0;
eqn3 = i2*(1 - n2 * c2* b^(n2-1))/(1-c2*b^(n2-1)) *(1-c2*b^(n2-1)*((b+c+mpl)/(c2*b^n2+c+mpl))) - i3 * (1-n3*c3*c^(n3-1)*((c+mpl)/(c3*c^n3+mpl))) == 0;
system = [eqn1,eqn2,eqn3];

Answer 1

[Yash]

Hi Samir,

To solve nonlinear equations in MATLAB, you can utilize the 'fsolve' function from the Optimization Toolbox. This function is specifically designed to find the roots of a system of nonlinear equations. By providing an initial guess, 'fsolve' attempts to converge to a solution that satisfies the equations.

The 'fsolve' function will attempt to find a solution for the system of equations starting from the initial guess provided. It will return the solution vector 'x' that satisfies the equations, or an error if it fails to converge.

To know more about the 'fsolve' function, refer to this documentation: https://in.mathworks.com/help/optim/ug/fsolve.html

Hope this helps!

[John D'Errico]

Not completely correct. fsolve is not built in. It is part of the optimization toolbox, and only available if you have that toolbox. In my humble opinion, it is one of the most useful toolboxes I have, but not everyone will have it.

[Yash]

Even I had that toolbox installed, updated the answer, thanks!

Answer 2

[Sam Chak]

Hi @Samir Thapa

If you have the Optimization Toolbox installed, then you can use the 'fsolve' function to solve the system of nonlinear equations. However, some nonlinear systems can have multiple solutions, depending on the initial guess values that are chosen.

% Solution set #1
x0a       = 1*[1, 1, 1];
[x, fval] = fsolve(@nonlinfcn, x0a)

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

x = 1×3

1.0e+04 * 1.4656 0.0032 0.0003

fval = 1×3

1.0e-12 * -0.0071 0.3837 -0.2025

% Solution set #2
x0b       = 2*[1, 1, 1];
[x, fval] = fsolve(@nonlinfcn, x0b)

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

x = 1×3

1.0e+04 * 0.0689 1.3999 0.0004

fval = 1×3

1.0e-10 * 0.0053 -0.6560 0.0006

function F = nonlinfcn(x)
    i1   = 310;
    i2   = 349.64;
    i3   = 353;
    c1   = 11.1984;
    n1   = 0.5067;
    c2   = 15.9867;
    n2   = 0.4271;
    c3   = 8.6028;
    n3   = 0.2449;
    mpl  = 308.5;
    F(1) = 48.62236629 - (x(1) + x(2) + x(3) + mpl)/mpl;
    F(2) = i1*(1 - n1*c1*x(1)^(n1 - 1))/(1 - c1*x(1)^(n1 - 1))*(1 - c1*x(1)^(n1 - 1)*((x(1) + x(2) + x(3) + mpl)/(c1*x(1)^n1 + x(2) + x(3) + mpl))) - i2*(1 - n2*c2*x(2)^(n2 - 1)*((x(2) + x(3) + mpl)/(c2*x(2)^n2 + x(3) + mpl)));
    F(3) = i2*(1 - n2*c2*x(2)^(n2 - 1))/(1 - c2*x(2)^(n2 - 1))*(1 - c2*x(2)^(n2 - 1)*((x(2)        + x(3) + mpl)/(c2*x(2)^n2        + x(3) + mpl))) - i3*(1 - n3*c3*x(3)^(n3 - 1)*((x(3)        + mpl)/(c3*x(3)^n3        + mpl)));
end

[Dyuman Joshi]

However, some nonlinear systems can have multiple solutions, and the output from fsolve will depend on the initial guess values that are chosen.

[Sam Chak]

@Dyuman Joshi, It's clearer now. Thanks!

Answer 3

[Walter Roberson]

There might be additional solutions.

Q = @(v) sym(v);

syms a b positive
syms c real
i1 = Q(310);
i2 = Q(34964) / Q(10)^2;
i3 = Q(353);
c1 = Q(111984) / Q(10)^4;
n1 = Q(5067) / Q(10)^4;
c2 = Q(159867) / Q(10)^4;
n2 = Q(4271) / Q(10)^4;
c3 = Q(86028) / Q(10)^4;
n3 = Q(2449) / Q(10)^4;
mpl = Q(308.5);
eqn1 = Q(4862236629) / Q(10)^8 - (a + b + c + mpl)/mpl == 0;
eqn2 = i1*(1 - n1 * c1* a^(n1-1))/(1-c1*a^(n1-1)) *(1-c1*a^(n1-1)*((a+b+c+mpl)/(c1*a^n1 + b+c+mpl))) - i2 * (1-n2*c2*b^(n2-1)*((b+c+mpl)/(c2*b^n2+c+mpl))) == 0;
eqn3 = i2*(1 - n2 * c2* b^(n2-1))/(1-c2*b^(n2-1)) *(1-c2*b^(n2-1)*((b+c+mpl)/(c2*b^n2+c+mpl))) - i3 * (1-n3*c3*c^(n3-1)*((c+mpl)/(c3*c^n3+mpl))) == 0;
system = ([eqn1; eqn2; eqn3])

system = 

start1 = [0.0689    1.3999    0.0004].' * 1e-4;
start2 = [9000 5000 4];
start3 = [14000 32 2];
sol1 = vpasolve(system, [a b c], start1)

sol1 = struct with fields:
    a: 689.05838496729307867145108434607
    b: 13998.717161477911435393065197254
    c: 3.7244540197954859354837183999334

sol2 = vpasolve(system, [a b c], start2)

sol2 = struct with fields:
    a: 9064.135929359262874364288421267
    b: 5623.1253454523787386377782441771
    c: 4.238725653358386997933334555938

sol3 = vpasolve(system, [a b c], start3)

sol3 = struct with fields:
    a: 14656.340055193821157737413004269
    b: 32.365749358899772536957181841301
    c: 2.7941959122790697256298138895309

[Sam Chak]

Thank you for providing the additional solutions. When dealing with complicated nonlinear systems, is there a recommended approach for determining the existence of a specific number of solutions? For instance, in this case, should we stop at three solution sets or continue searching for more?

[Walter Roberson]

When there are polynomials, you can get a hint as the product of the maximum orders of each variable.

However in this particular equation there are values raised to non-integer powers. Those require numeric solutions for most practical purposes. If you convert the non-integer power into a rational (as I do above) then the polynomial orders get pretty large ... you would have to be working with polynomials of degree 10000 or more for this system. If you do manage to get a root() output then you could at least cross-check each of them... might take time though.

[John D'Errico]

@Sam Chak

Sorry, there is no hard rule that says you can positively stop looking at some point. Yes, a cubic polynomial has exactly 3 solutions. These are not cubics though. These equations are implicitly equivalent to a VERY high order polynomial. Those fractional powers make it so, if you could actually reduce the problem to a single equation in one unknown. You can't do so. And you can't even really know what order that impicit polynomial would be in such a case.

[Sam Chak]

Thanks @Walter Roberson and @John D'Errico for the explanations. If you look at my code, you can already guess that my initial values were purely lucky guesses. I tried searching randomly, but complex-valued solutions were returned.