Maths / ODE / Why do I get complex results?

Hello dear Matlab Community,

I have a question that is a more a mathematics- than a MATLAB question but since I spent whole day on it and don't know who else to ask I thought I would post it here.

When running the code below I get complex values for my 4 Characteristics. variables... any idea why?

Many thanks!

Hélène

PAR.POSITION_IN_PAPERMACHINE=[0;1;2;3;4;5];
for k=2:numel(PAR.POSITION_IN_PAPERMACHINE)
    PAR.POSITION_IN_PAPERMACHINE(k)=PAR.POSITION_IN_PAPERMACHINE(k)+PAR.POSITION_IN_PAPERMACHINE(k-1);
end

PAR.HUM.PAPER_COAT=0.47;
PAR.CONDITIONS.PAPER_INI=[1.1;70];
PAR.TEMP.STEAM_CYL_ALL=[5;2;9;2;8];
PAR.AREA.COND_ALL=[3;2;4;5;6];

Characteristics.paper_all=zeros([],2); Characteristics.paper_begsection_all=[];Support.differentialua=zeros([],1);PAR.CONDITIONS.PAPER_INI_ALL=[];PAR.CONDITIONS.PAPER_COAT_ALL=[];

% Solving loop
for k=1:numel(PAR.POSITION_IN_PAPERMACHINE)-1        
    PAR.TEMP.STEAM_CYL=PAR.TEMP.STEAM_CYL_ALL(k); PAR.AREA.COND=PAR.AREA.COND_ALL(k);
         
    % Solving the ODE system
    if k<3
        [Position.paper,Characteristics.paper]=ode45(@(L,y) myODE(L,y,PAR.AREA.COND,PAR.TEMP.STEAM_CYL),PAR.POSITION_IN_PAPERMACHINE(k:k+1),PAR.CONDITIONS.PAPER_INI); PAR.CONDITIONS.PAPER_INI=Characteristics.paper(end,:);  
        PAR.CONDITIONS.PAPER_COAT(1)=PAR.HUM.PAPER_COAT; 
        PAR.CONDITIONS.PAPER_COAT(2)=PAR.CONDITIONS.PAPER_INI(2);
    else
        [Position.paper,Characteristics.paper]=ode45(@(L,y) myODE(L,y,PAR.AREA.COND,PAR.TEMP.STEAM_CYL),PAR.POSITION_IN_PAPERMACHINE(k:k+1),PAR.CONDITIONS.PAPER_COAT);   PAR.CONDITIONS.PAPER_COAT=Characteristics.paper(end,:); 
    end

% Value of u and Tp at beginning of a cylinder&no-cylinder section 
Characteristics.paper_begsection=Characteristics.paper(end-40,:); 

% Differential of u 
if k<3
    Support.differentialu=PAR.CONDITIONS.PAPER_INI(:,1)-Characteristics.paper_begsection(:,1); 
else; Support.differentialu=PAR.CONDITIONS.PAPER_COAT(:,1)-Characteristics.paper_begsection(:,1); 
end

Characteristics.paper_all=[Characteristics.paper_all;Characteristics.paper];PAR.CONDITIONS.PAPER_INI_ALL=[PAR.CONDITIONS.PAPER_INI_ALL;PAR.CONDITIONS.PAPER_INI];PAR.CONDITIONS.PAPER_COAT_ALL=[PAR.CONDITIONS.PAPER_COAT_ALL;PAR.CONDITIONS.PAPER_COAT];Characteristics.paper_begsection_all=[Characteristics.paper_begsection_all;Characteristics.paper_begsection];Support.differentialua=[Support.differentialua;Support.differentialu];
end

function dudL = myODE1(~,u,Tp,~,~)
Psatp = 0.133322*exp(18.3036-3816.44/(Tp+227.03));  
Phi= 1-exp(-(47.58*u^1.877+0.10085*Tp*u^1.0585));
Ppw= Psatp*Phi;                                     

dudL=-(1/Tp)*log(3/(2-Ppw));  
end

function dTpdL=myODE2(~,u,Tp,Areacond,Ts)
Psatp = 0.133322*exp(18.3036-3816.44/(Tp+227.03));  
Phi= 1-exp(-(47.58*u^1.877+0.10085*Tp*u^1.0585));
Ppw= Psatp*Phi ;

DHs = 0.10085*(u^1.0585)*((Tp+273.15)^2)*(1-Phi)/Phi;
DHvap = 2501-2.3237*Tp;  
DHevap = DHs+DHvap; 
Usp=5*u;

dTpdL=0.001*DHevap*Usp*(Ts-Tp)*Areacond+(100-Tp)-(1/Tp)*log(3/(2-Ppw));                                                                                                                                                                    
end

function dy = myODE(L,y,Areacond,Ts)
u=y(1); Tp=y(2);
dudL = myODE1(L,u,Tp,Areacond,Ts);
dTpdL = myODE2(L,u,Tp,Areacond,Ts);
dy = [dudL;dTpdL];
end

Answer 1

[Steven Lord]

Each of your ODE functions includes this line:

Phi= 1-exp(-(47.58*u^1.877+0.10085*Tp*u^1.0585));

For sake of argument, let's let Tp = 0. Then this simplifies to:

Phi= 1-exp(-(47.58*u^1.877));

What does this return when you evaluate it for a negative value of u?

u = -1;
Phi= 1-exp(-(47.58*u^1.877))

Specifically, what does the term inside the exp call return?

-(47.58*(-1)^1.877)

If you let Tp be non-zero the same argument holds, only you're raising that negative number to two non-integer-valued powers.

[Lenilein]

Thank you Steven! Now I just have to figure out why u is taking negative values :D