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);endPAR.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); endCharacteristics.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];endfunction 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)); endfunction 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)); endfunction 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
Best Answer