Alternative to add assumptions to vpasolve
Show older comments
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
More Answers (0)
Categories
Find more on Calculus 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!