Thread Subject: non-linear system : can't get a proper solution

T = input(' T = '); % T should be set at around 1000
T_aux=1/T*1e04;
ln_Kp1 = 2.5904*T_aux - 29.077; % meglio interpolare in un range di T + ristretto tipo 600-1200K
Kp1 = exp(ln_Kp1);
ln_Kp2 = 0.4694*T_aux - 4.3578;
Kp2 = exp(ln_Kp2);
ln_Kp3 = -0.9795*T_aux + 12.078;
Kp3 = exp(ln_Kp3);
ln_Kp4 = -1.4449*T_aux + 20.379;
Kp4 = exp(ln_Kp4);
SBR = input(' SBR = '); % SBR should be set anywhere in (0,1]
q = 2.1587*SBR;
qtot = q+1.505/11;
ni0 = [1;0;qtot;0;12.5/11;5/11;0];
a=[-1,-1,0,1; 0,1,0,0; 1,-1,0,-1; 1,0,-1,0; -3,1,2,1.4; 0,0,1,0; 0,0,0,-0.1];
options = optimset('TolFun',1e-14);

[x]=fsolve(@(x)sist4eq(x,ni0,Kp1,Kp2,Kp3,Kp4),[0.05;0.02;0.02;0.027],options);

disp(x);

ni = a*x + ni0;
disp(ni); % every component of vector ni must be real positive

% if T=1073.15 and SBR near to 0 , the result for ni should be similar
% to ni = ( 17.7 , 4.37 , 4.47 , 4.44 , 18.13 , 7.7 , 0.36 )' .

% Using fsolve I only get complex solutions ( for x as for ni )...

PLEASE HELP ME...SHOULD I CHANGE THE STARTING VECTOR IN fsolve? I TRIED NEWTON METHOD AS WELL AND THE RESULTS IS NaN ( not a number )...Jacobian badly scaled or near to singular...PLEASE TELL ME WHAT TO CHANGE TO MAKE IT WORK

I write here below the file named sist4eq that figures in fsolve call :


function f=sist4eq(x,ni0,Kp1,Kp2,Kp3,Kp4)

f=[ (((x(1)-x(3))*(ni0(3)+x(1)-x(2)-x(4))*((sum(ni0)-2*x(1)+2*x(3)+1.3*x(4))^2))/...
        ((ni0(1)-x(1)-x(2)+x(4))*((ni0(5)-3*(x(1)+x(2)+2*x(3)+1.4*x(4)))^3))-Kp1 );
 ( (x(2)*(ni0(5)-3*(x(1)+x(2)+2*x(3)+1.4*x(4)) ))/...
          ((ni0(1)-x(1)-x(2)+x(4))*(ni0(3)+(x(1)-x(2)-x(4))))-Kp2 );
( (((ni0(5)-3*(x(1)+x(2)+2*x(3)+1.4*x(4)))^2)*(1/(sum(ni0)-2*x(1)+2*x(3)+1.3*x(4))))/...
         ((x(1)-x(3)))-Kp3);
     ( ((ni0(1)-x(1)-x(2)+x(4))*((ni0(5)-3*x(1)+x(2)+2*x(3)+1.4*x(4))^1.4)*(1/((sum(ni0)-2*x(1)+2*x(3)+1.3*x(4))^1.3)))/...
         (((-0.1*x(4))^0.1)*(ni0(3)+x(1)-x(2)-x(4)))-Kp4);]

SORRY IF THE LAST 2 eq. ARE ON 3 ROWS BUT I CAN'T DO BETTER WITH THIS EDITOR .