Solve a System of five non linear equation

I want to solve numerically a system of five non linear equations. I use the solve equation but instead of numerical expressions for the parameters, matlab returns a symbol for each parameter. For the parameters a1, mu1, a2, mu2, and w1 I need them to be positive and real, also 0<=w1<=1. My question is, can the following system return numerical expressions for the parameters of interest? If so, is there an option in the solve() function that I need to include?

C1=2.083333333333334;
C2=5.514861889710719;
C3=15.477300064209032;
C4=44.665706966377420;
C5=1.313098986326810e+02;

syms a1 mu1 a2 mu2 w1

E_i=[( w1*gamma(mu1+1/a1)/( ( mu1^(1/a1) )*gamma(mu1) ) )...
    +( (1-w1)*gamma(mu2+1/a2)/( ( mu2^(1/a2) )*gamma(mu2) ) )==C1,...
    ( w1*gamma(mu1+2/a1)/( ( mu1^(2/a1) )*gamma(mu1) ) )...
    +( (1-w1)*gamma(mu2+2/a2)/( ( mu2^(2/a2) )*gamma(mu2) ) )==C2,...
    ( w1*gamma(mu1+3/a1)/( ( mu1^(3/a1) )*gamma(mu1) ) )...
    +( (1-w1)*gamma(mu2+3/a2)/( ( mu2^(3/a2) )*gamma(mu2) ) )==C3,...
    ( w1*gamma(mu1+4/a1)/( ( mu1^(4/a1) )*gamma(mu1) ) )...
    +( (1-w1)*gamma(mu2+4/a2)/( ( mu2^(4/a2) )*gamma(mu2) ) )==C4,...
    ( w1*gamma(mu1+5/a1)/( ( mu1^(5/a1) )*gamma(mu1) ) )...
    +( (1-w1)*gamma(mu2+5/a2)/( ( mu2^(5/a2) )*gamma(mu2) ) )==C5
    ];

assume(a1>0 & mu1>0 & a2>0 & mu2>0 & 0<=w1 & w1<=1) % assumptions to narrow results
S=solve(E_i,'Real',true,'ReturnConditions', true);

[Torsten]

What is the output if you type

S.a1

S.a2

S.mu1

S.mu2

S.w1

after the solve command ?

[vaggelis papasot]

@Torsten thank you for replying. The solve command gives the following struct

S = 

struct with fields:

a1: [32×1 sym]
a2: [32×1 sym]
mu1: [32×1 sym]
mu2: [32×1 sym]
w1: [32×1 sym]
parameters: [1×5 sym]
conditions: [32×1 sym]

Where for example

S.a1(1) 
ans =z
and 
S.a2(4) 
ans =z1

[Torsten]

And what is z and z1 ?

What is in the struct "conditions" ?

Could you output one coherent solution S.a1(1),S.a2(1),...

[vaggelis papasot]

S.a1(1)
ans =z
S.a2(1)
ans =z1
S.mu1(1)
ans =z2
S.mu2(1)
ans =z3
S.w1(1)
ans =z4

The condition struct contains 32 expressions on of them is

S.conditions(1)
ans =
gamma(z3 + 1/z1)/(z3^(1/z1)*gamma(z3)) - 25/12 == 0 & gamma(z3 + 5/z1)/(z3^(5/z1)*gamma(z3)) - 4620056332439061/35184372088832 == 0 & gamma(z3 + 4/z1)/(z3^(4/z1)*gamma(z3)) - 6286139414063035/140737488355328 == 0 & gamma(z2 + 2/z) == 0 & gamma(z2 + 3/z) == 0 & gamma(z2 + 4/z) == 0 & gamma(z2 + 5/z) == 0 & z4 == 0 & 0 < z & 0 < z1 & 0 < z2 & 0 < z3 & 0 <= z4 & z4 <= 1 & gamma(z3 + 2/z1)/(z3^(2/z1)*gamma(z3)) - 1552295621968809/281474976710656 == 0 & gamma(z3 + 3/z1)/(z3^(3/z1)*gamma(z3)) - 272279542194817/17592186044416 == 0

I suppose if you run the script you will get the same values in the fields of the S struct.

[Torsten]

Remove all the conditions you imposed and simply try

S = vpasolve(E_i)

Answer 1

[John D'Errico]

E_i=[( w1*gamma(mu1+1/a1)/( ( mu1^(1/a1) )*gamma(mu1) ) )+( (1-w1)*gamma(mu2+1/a2)/( ( mu2^(1/a2) )*gamma(mu2) ) )==C1,...
    ( w1*gamma(mu1+2/a1)/( ( mu1^(2/a1) )*gamma(mu1) ) )+( (1-w1)*gamma(mu2+2/a2)/( ( mu2^(2/a2) )*gamma(mu2) ) )==C2,...
    ( w1*gamma(mu1+3/a1)/( ( mu1^(3/a1) )*gamma(mu1) ) )+( (1-w1)*gamma(mu2+3/a2)/( ( mu2^(3/a2) )*gamma(mu2) ) )==C3,...
    ( w1*gamma(mu1+4/a1)/( ( mu1^(4/a1) )*gamma(mu1) ) )+( (1-w1)*gamma(mu2+4/a2)/( ( mu2^(4/a2) )*gamma(mu2) ) )==C4,...
    ( w1*gamma(mu1+5/a1)/( ( mu1^(5/a1) )*gamma(mu1) ) )+( (1-w1)*gamma(mu2+5/a2)/( ( mu2^(5/a2) )*gamma(mu2) ) )==C5
    ];

numel(E_i)
ans =
5

So 5 equations, in the 5 unknowns a1,a2,mu1,mu2,w1.

S=vpasolve(E_i)
S = 
struct with fields:

a1: [0×1 sym]
a2: [0×1 sym]
mu1: [0×1 sym]
mu2: [0×1 sym]
w1: [0×1 sym]

But nothing insures that a solution always exists for any set of equations you may choose to pose. vpasolve did not find one. Nor did solve. Solve just returned a result that effectively said the solution is the solution to the system of equations listed in the conditions. vpasolve tried, but gave up.

In fact, I would postulate that certainly there will be no explicit (algebraic) solution that solve could find. you cannot plead with a solver. You cannot demand that it do better.

That does not mean no solution exists, only that a numerical solver was unable to find it based on the default starting values. It might have been successful, it you give it better guesses (than the default) for those starting values.

Answer 2

[William Rose]

@vaggelis papasot,

You can get numeric values for the unknown quatities by using fsolve() and by not using "syms". Rewrite your equations, and instead of uing a1, mu1, a2, mu2, w1, use x(1:5). Use fsolve to estimate the vector of unknowns, x. See the help for fsolve().