Alternative to add assumptions to vpasolve
4 views (last 30 days)
Show older comments
vaggelis papasot
on 19 Feb 2022
Commented: vaggelis papasot
on 22 Feb 2022
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);
0 Comments
Accepted Answer
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)
More Answers (0)
See Also
Categories
Find more on Gamma Functions in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!