Alternative to add assumptions to vpasolve

4 views (last 30 days)
I need to numerically solve the following system of equations, in order to calculate the parameters w1,w2,w3,w4,a1,a2,a3,a4,mu1,mu2,mu3, and mu4. Also, I need the constraint that w1+w2+w3+w4=1. I know that vpasolve does not use assumptions and that one can only define ranges for the parameters that need to be calculated. I have the following question. Is there a way to define the constaint w1+w2+w3+w4=1 to get a numerical solution with vpasolve?
C1=0.9945;
C2=1;
C3=1.0162;
C4=1.0431;
C5=1.0809;
C6=1.13;
C7=1.1911;
C8=1.2650;
C9=1.3530;
C10=1.4564;
C11=1.5769;
C12=1.7165;
syms w1 a1 mu1 w2 a2 mu2 w3 a3 mu3 w4 a4 mu4
eq1=w1*( gamma( mu1+(1/a1) )/( ( mu1^(1/a1) )*gamma(mu1) ) )...
+w2*( gamma(mu2+(1/a2))/( ( mu2^(1/a2) )*gamma(mu2) ) )...
+w3*( gamma( mu3+(1/a3) )/( ( mu3^(1/a3) )*gamma(mu3) ) )...
+w4*( gamma( mu1+(1/a4) )/( ( mu4^(1/a4) )*gamma(mu4) ) )==C1;
eq2=w1*( gamma( mu1+(2/a1) )/( ( mu1^(2/a1) )*gamma(mu1) ) )...
+w2*( gamma(mu2+(2/a2))/( ( mu2^(2/a2) )*gamma(mu2) ) )...
+w3*( gamma( mu3+(2/a3) )/( ( mu3^(2/a3) )*gamma(mu3) ) )...
+w4*( gamma( mu1+(2/a4) )/( ( mu4^(2/a4) )*gamma(mu4) ) )==C2;
eq3=w1*( gamma( mu1+(3/a1) )/( ( mu1^(3/a1) )*gamma(mu1) ) )...
+w2*( gamma(mu2+(3/a2))/( ( mu2^(3/a2) )*gamma(mu2) ) )...
+w3*( gamma( mu3+(3/a3) )/( ( mu3^(3/a3) )*gamma(mu3) ) )...
+w4*( gamma( mu1+(3/a4) )/( ( mu4^(3/a4) )*gamma(mu4) ) )==C3;
eq4=w1*( gamma( mu1+(4/a1) )/( ( mu1^(4/a1) )*gamma(mu1) ) )...
+w2*( gamma(mu2+(4/a2))/( ( mu2^(4/a2) )*gamma(mu2) ) )...
+w3*( gamma( mu3+(4/a3) )/( ( mu3^(4/a3) )*gamma(mu3) ) )...
+w4*( gamma( mu1+(4/a4) )/( ( mu4^(4/a4) )*gamma(mu4) ) )==C4;
eq5=w1*( gamma( mu1+(5/a1) )/( ( mu1^(5/a1) )*gamma(mu1) ) )...
+w2*( gamma(mu2+(5/a2))/( ( mu2^(5/a2) )*gamma(mu2) ) )...
+w3*( gamma( mu3+(5/a3) )/( ( mu3^(5/a3) )*gamma(mu3) ) )...
+w4*( gamma( mu1+(5/a4) )/( ( mu4^(5/a4) )*gamma(mu4) ) )==C5;
eq6=w1*( gamma( mu1+(6/a1) )/( ( mu1^(6/a1) )*gamma(mu1) ) )...
+w2*( gamma(mu2+(6/a2))/( ( mu2^(6/a2) )*gamma(mu2) ) )...
+w3*( gamma( mu3+(6/a3) )/( ( mu3^(6/a3) )*gamma(mu3) ) )...
+w4*( gamma( mu1+(6/a4) )/( ( mu4^(6/a4) )*gamma(mu4) ) )==C6;
eq7=w1*( gamma( mu1+(7/a1) )/( ( mu1^(7/a1) )*gamma(mu1) ) )...
+w2*( gamma(mu2+(7/a2))/( ( mu2^(7/a2) )*gamma(mu2) ) )...
+w3*( gamma( mu3+(7/a3) )/( ( mu3^(7/a3) )*gamma(mu3) ) )...
+w4*( gamma( mu1+(7/a4) )/( ( mu4^(7/a4) )*gamma(mu4) ) )==C7;
eq8=w1*( gamma( mu1+(8/a1) )/( ( mu1^(8/a1) )*gamma(mu1) ) )...
+w2*( gamma(mu2+(8/a2))/( ( mu2^(8/a2) )*gamma(mu2) ) )...
+w3*( gamma( mu3+(8/a3) )/( ( mu3^(8/a3) )*gamma(mu3) ) )...
+w4*( gamma( mu1+(8/a4) )/( ( mu4^(8/a4) )*gamma(mu4) ) )==C8;
eq9=w1*( gamma( mu1+(9/a1) )/( ( mu1^(9/a1) )*gamma(mu1) ) )...
+w2*( gamma(mu2+(9/a2))/( ( mu2^(9/a2) )*gamma(mu2) ) )...
+w3*( gamma( mu3+(9/a3) )/( ( mu3^(9/a3) )*gamma(mu3) ) )...
+w4*( gamma( mu1+(9/a4) )/( ( mu4^(9/a4) )*gamma(mu4) ) )==C9;
eq10=w1*( gamma( mu1+(10/a1) )/( ( mu1^(10/a1) )*gamma(mu1) ) )...
+w2*( gamma(mu2+(10/a2))/( ( mu2^(10/a2) )*gamma(mu2) ) )...
+w3*( gamma( mu3+(10/a3) )/( ( mu3^(10/a3) )*gamma(mu3) ) )...
+w4*( gamma( mu1+(10/a4) )/( ( mu4^(10/a4) )*gamma(mu4) ) )==C10;
eq11=w1*( gamma( mu1+(11/a1) )/( ( mu1^(11/a1) )*gamma(mu1) ) )...
+w2*( gamma(mu2+(11/a2))/( ( mu2^(11/a2) )*gamma(mu2) ) )...
+w3*( gamma( mu3+(11/a3) )/( ( mu3^(11/a3) )*gamma(mu3) ) )...
+w4*( gamma( mu1+(11/a4) )/( ( mu4^(11/a4) )*gamma(mu4) ) )==C11;
eq12=w1*( gamma( mu1+(12/a1) )/( ( mu1^(12/a1) )*gamma(mu1) ) )...
+w2*( gamma(mu2+(12/a2))/( ( mu2^(12/a2) )*gamma(mu2) ) )...
+w3*( gamma( mu3+(12/a3) )/( ( mu3^(12/a3) )*gamma(mu3) ) )...
+w4*( gamma( mu1+(12/a4) )/( ( mu4^(12/a4) )*gamma(mu4) ) )==C12;
eqn=[eq1 eq2 eq3 eq4 eq5 eq6 eq7 eq8 eq9 eq10 eq11 eq12];
vars=[w1,a1,mu1,w2,a2,mu2,w3,a3,mu3,w4,a4,mu4];
range=[0 1; 0 Inf; 0 Inf; 0 1; 0 Inf; 0 Inf; 0 1; 0 Inf;0 Inf; 0 1; 0 Inf; 0 Inf];
S=vpasolve(eqn,vars,range,'Random',true);

Accepted Answer

Torsten
Torsten on 19 Feb 2022
I don't know if your system of equations has a solution, but you can enforce sum(w(i))=1 by replacing
wi by wi/(w1+w2+w3+w4) (i=1,...,4)
  1 Comment
vaggelis papasot
vaggelis papasot on 22 Feb 2022
Thank you @Torsten. This is a good workaround for this problem. However, the system now yields no solution.

Sign in to comment.

More Answers (0)

Products


Release

R2020a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!