Negative output Fsolve Function
Afficher commentaires plus anciens
Hi all,
can anyone help me understand why I am getting a negative output from this code and how can I modify the code so I do not get a negative output?
function [x,exitflag]= Function_final(c)
c=[10^-7,10^-7,0.5*10^-3,0.5*10^-3,0.25*10^-3,0.25*10^-3,0.25*10^-3,0.25*10^-3,0.5*10^-3,0.5*10^-3,0.5*10^-3,0.5*10^-4,0.25*10^-1,0.5*10^-4,0.5*10^-6];
MaT = c(3)+c(4)+c(9)+c(10)+c(11)+c(12);
MbT = c(5)+c(6)+c(7)+c(8)+c(9)+c(10)+c(11);
McT = c(2)+c(6)+2*c(7)+3*c(8)+c(9)+2*c(10)-c(1)-c(4)+c(13)+2*c(15)+c(12);
MdT = c(12)+c(13)+c(14)+c(15);
z=[1,1,2,1,3,2,1,0,0,1,1,1,1,2,0];
T=36;
I=0.5*sum(c.*z.^2);
A=0.486+(6.07*10^-4)*T+(6.43*10^-6)*T^2;
gamma_initial=zeros(1,15);
for i=1:15
gamma_initial(i)=10^(-A*z(i)^2*(I^0.5/(1+I^0.5)-0.3*I));
end
p = [gamma_initial, c];
options=optimoptions('fsolve','Display','iter');
[x,exitflag]=fsolve(@(x)Approach(x,MaT,MbT,McT,MdT),p,options);
end
function F = Approach(p,MaT,MbT,McT,MdT)
K1 = 1.01*10^-14;
K2 = 2.13*10^-13;
K3 = 1.14*10^12;
K4 = 1.73*10^19;
K5 = 2.71*10^21;
K6 = 5.36*10^14;
K7 = 3.41*10^20;
K8 = 1.18*10^7;
K9 = 4.01;
K10= 1.74*10^10;
K11 = 3.09*10^3;
z=[1,1,2,1,3,2,1,0,0,1,1,1,1,2,0];
T=36;
A=0.486+(6.07*10^-4)*T+(6.43*10^-6)*T^2;
c = p(16:30);
la = log10(p(1:15));
I=0.5*sum(c.*z.^2);
F(1) = la(1)^2 - la(1) * (I*A*z(1)^2) * 2.6 - 1/(1-I)^2 * A^2 * z(1)^4 * (I - 1.69 * I^2 - 0.09 * I^4);
F(2) = la(2)^2 - la(2) * (I*A*z(1)^2) * 2.6 - 1/(1-I)^2 * A^2 * z(1)^4 * (I - 1.69 * I^2 - 0.09 * I^4);
F(3) = la(3)^2 - la(3) * (I*A*z(1)^2) * 2.6 - 1/(1-I)^2 * A^2 * z(1)^4 * (I - 1.69 * I^2 - 0.09 * I^4);
F(4) = la(4)^2 - la(4) * (I*A*z(1)^2) * 2.6 - 1/(1-I)^2 * A^2 * z(1)^4 * (I - 1.69 * I^2 - 0.09 * I^4);
F(5) = la(5)^2 - la(5) * (I*A*z(1)^2) * 2.6 - 1/(1-I)^2 * A^2 * z(1)^4 * (I - 1.69 * I^2 - 0.09 * I^4);
F(6) = la(6)^2 - la(6) * (I*A*z(1)^2) * 2.6 - 1/(1-I)^2 * A^2 * z(1)^4 * (I - 1.69 * I^2 - 0.09 * I^4);
F(7) = la(7)^2 - la(7) * (I*A*z(1)^2) * 2.6 - 1/(1-I)^2 * A^2 * z(1)^4 * (I - 1.69 * I^2 - 0.09 * I^4);
F(8) = la(8)^2 - la(8) * (I*A*z(1)^2) * 2.6 - 1/(1-I)^2 * A^2 * z(1)^4 * (I - 1.69 * I^2 - 0.09 * I^4);
F(9) = la(9)^2 - la(9) * (I*A*z(1)^2) * 2.6 - 1/(1-I)^2 * A^2 * z(1)^4 * (I - 1.69 * I^2 - 0.09 * I^4);
F(10) = la(10)^2 - la(10) * (I*A*z(1)^2) * 2.6 - 1/(1-I)^2 * A^2 * z(1)^4 * (I - 1.69 * I^2 - 0.09 * I^4);
F(11) = la(11)^2 - la(11) * (I*A*z(1)^2) * 2.6 - 1/(1-I)^2 * A^2 * z(1)^4 * (I - 1.69 * I^2 - 0.09 * I^4);
F(12) = la(12)^2 - la(12) * (I*A*z(1)^2) * 2.6 - 1/(1-I)^2 * A^2 * z(1)^4 * (I - 1.69 * I^2 - 0.09 * I^4);
F(13) = la(13)^2 - la(13) * (I*A*z(1)^2) * 2.6 - 1/(1-I)^2 * A^2 * z(1)^4 * (I - 1.69 * I^2 - 0.09 * I^4);
F(14) = la(14)^2 - la(14) * (I*A*z(1)^2) * 2.6 - 1/(1-I)^2 * A^2 * z(1)^4 * (I - 1.69 * I^2 - 0.09 * I^4);
F(15) = la(15)^2 - la(15) * (I*A*z(1)^2) * 2.6 - 1/(1-I)^2 * A^2 * z(1)^4 * (I - 1.69 * I^2 - 0.09 * I^4);
F(16) = (K1*(10^la(2)*p(17))^(-1))-(10^la(1)*p(16));
F(17) = (K2*(10^la(3)*p(18))*(10^la(1)*p(16)))-(10^la(4)*p(19));
F(18) = (K3*(10^la(5)*p(20))*(10^la(2) *p(17)))-(10^la(6)*p(21));
F(19) = (K4*(10^la(5)*p(20))*(10^la(2)*p(17))^2)-(10^la(7)*p(22));
F(20) = (K5*(10^la(5)*p(20))*(10^la(2)*p(17))^3)-(10^la(8)*p(23));
F(21) = (K6*(10^la(3)*p(18))*(10^la(5)*p(20))*(10^la(2)*p(17)))-(10^la(9)*p(24));
F(22) = (K7*(10^la(3)*p(18))*(10^la(5)*p(20))*(10^la(2)*p(17))^2)-(10^la(10)*p(25));
F(23) = (K8*(10^la(3)*p(18))*(10^la(5)*p(20)))-(10^la(11)*p(26));
F(24) = (K9*(10^la(3)*p(18))*(10^la(14)*p(29))*(10^la(2)*p(17)))-(10^la(12)*p(27));
F(25) = (K10*(10^la(2)*p(17))*(10^la(14)*p(29))) - (10^la(13)*p(28));
F(26) = (K11*((10^la(2)*p(17))^2)*(10^la(14)*p(29))) - (10^la(15)*p(30));
F(27) = p(18)+p(19)+p(24)+p(25)+p(26)-MaT;
F(28) = p(20)+p(21)+p(22)+p(23)+p(24)+p(25)+p(26)-MbT;
F(29) = p(17)+p(21)+2*p(22)+3*p(23)+p(24)+2*p(25)-p(16)-p(19)-McT;
F(30) = p(27)+p(28)+p(29)+p(30)-MdT;
F = [F(1);F(2);F(3);F(4);F(5);F(6);F(7);F(8);F(9);F(10);F(11);F(12);F(13);F(14);F(15);F(16);F(17);F(18);F(19);F(20);F(21);F(22);F(23);F(24);F(25);F(26);F(27);F(28);F(29);F(30);
end
Réponse acceptée
Plus de réponses (0)
Catégories
En savoir plus sur Data Exploration dans Centre d'aide et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!