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
Best Answer