I run fsolve for interval and only for one x it gives me wrong value

Hi,

so I have this unexplainalbe thing for me. I have system of equations.

%Posylyt
res=zeros(1,2*(Ncell-1+Ncell-1+Ncell));
res(1)=QTOT-qIM_N(1)-qC_N(1);
for i=1:Ncell-2
    res(1+i)=qIM_N(i)-qC_N(i+1)-qIM_N(i+1);
    res(Ncell-1+i)=qC_N(1+i)-qOM_N(i)+qOM_N(i+1); 
end
res(2*Ncell-2)=qIM_N(Ncell-1)-qC_N(Ncell);
res(2*Ncell-1)=qC_N(Ncell)-qOM_N(Ncell-1); 
for i=1:Ncell-1
    res(2*Ncell-1+i)=DpOM_N(i)+DpCtotal_N(1+i)+DpIM_N(i)-DpCtotal_N(i);
end

%Negalyt
res(3*Ncell-1)=QTOT-qIM_P(1)-qC_P(1);
for i=1:Ncell-2
    res(3*Ncell-1+i)=qIM_P(i)-qC_P(i+1)-qIM_P(i+1);
    res(4*Ncell-3+i)=qC_P(1+i)-qOM_P(i)+qOM_P(i+1);
end
res(4*Ncell-3+Ncell-1)=qIM_P(Ncell-1)-qC_P(Ncell);
res(4*Ncell-3+Ncell)=qC_P(Ncell)-qOM_P(Ncell-1);
for i=1:Ncell-1
    res(4*Ncell-3+Ncell+i)=DpOM_P(i)+DpCtotal_P(1+i)+DpIM_P(i)-DpCtotal_P(i);
end

where I have closed system and balance inflow and outflow flowrate + I need extra equations where I say that pressure drop in closed loops is zero

pressures are also functions of flowrate + some of them depends on length of pipe

Posylyt and Negalyt are separate independent identical systems for now and shares only constants

I run a for cycle with fsolve where I calculate fsolve for different length of pipes - so only some of pressure drops will be influenced

I got a nice chart, but here comes weird thing! For only one value and only for system Negalyt, I got completely wrong value. And then it goes as should be.

I have no idea why, but it happens always only for one value and for specific initial input flow rate

The system posylyt goes exactly the same behaviour, but without this error. And I have no idea, why, because to me all equations and other set ups seems completely identical. And I assume if there would be an error in some pressure equation it woul get then different behaviour so maybe it is just fsolve setting? Any experience?

options = optimoptions('fsolve','Display','off','TolFun',1e-15,'TolX',1e-15,'MaxFunEvals',100000,'MaxIter',10000);

[Walter Roberson]

It would help if we had your full code and data to test with.

[Martin Pecha]

So I ran a debug and found out, that my fsolve is extremely sensitive to x0 initial values

x0=zeros(Ncell,1);
%Posylyt
for i=1:Ncell
   x0(i)=QTOT/Ncell;
end
for i=1:Ncell-1
   x0(Ncell+i)=QTOT/(Ncell-i);
end
for i=1:Ncell-1
    x0(2*Ncell-1+i)=QTOT/i;
end
%NEGALYT
for i=1:Ncell
   x0(3*Ncell-2+i)=QTOT/Ncell; %I had here QTOT/i 
end
for i=1:Ncell-1
   x0(4*Ncell-2+i)=QTOT/(Ncell-i);
end
for i=1:Ncell-1
    x0(5*Ncell-3+i)=QTOT/i;
end

so maybe some idea for someone in the future - CHECK YOUR initial X0