complex numbers as function output

The functions used to solve my ODE's using Euler's method give real values in the first iteration ; where as the function returns complex variables from the third iteration. I'm not able to figure out the reason behind this.
function file
function f =f_UP(K)
global M Cp R P z t A tda
P = (R*K(1)*K(2));
% z = (1.4888*(t^0.75)*(Cp^0.25)*(A^0.25)*(K(3)^0.25)*(K(2)^0.25)*((K(2)-tda)^0.25))
% z = 1.4888*(t^0.75)*(Cp^0.25)*(A^0.25)*(K(3)^0.25)*(K(2)^0.25)*((K(2)-tda)^0.25);
f(1) = ((-K(1))*(1.113-2.14*z +2.81*z^2-1.76*z^3+0.396*z^4+ ...
0.08694*z^5-0.0579*z^6+0.0077*z^7))-((K(3))*(0.4235+1.5*z-1.133*z^2+0.3741*z^3-...
0.04537*z^4));
f(2)= (((-1)*R*(K(2)))./(P*M*Cp))*(760.91-1608.08*(z)+2449.22*z^2-...
(2001.61*z^3)+ 904.51*z^4 -212.087*z^5+20.095*z^6);
f(3)= ((-1)/((K(1))*(K(3))))*((1.113-2.14*z+2.81*z^2-1.76*z^3+ ...
0.396*z^4+0.08694*z^5-0.0579*z^6 +0.0077*z^7));
f=f'
Main file
clear all
close all
global M Cp R A tda z t
deltt =0.1;
t_span = 0:1;
x_H2O = 0.5;
x_CO2 = 0.5;
M_CO2 = 44*(10^-3); %kg/mol
M_H2O = 18*(10^-3); %kg/mol
M = (x_CO2*M_CO2) + (x_H2O*M_H2O);
A = 3.14*(0.091)^2;
tda = 300;
R = 8.314; %m3.Pa/mol/K
cp_CO2 = 1168; %J/kg/K
cp_H2O = 2147; %J/kg/K
Cp = (x_CO2*cp_CO2)+(x_H2O*cp_H2O)
K(1,:) = [0.43545,800,1.2]
t = 0;
tend =10;
td = 0:tend;
or i= 1:length(td)
z = 1.4888*(td(i)^0.75)*(Cp^0.25)*(A^0.25)*(K(3)^0.25)*(K(2)^0.25)*((K(2)-tda)^0.25)
K(i+1,:)= K(i,:)' + deltt*f_UP(K(i,:))
g(i)= z
% t = t+1
end

2 Commenti

dpb
dpb il 22 Ott 2019
I'd guess z goes negative...use debugger to step through and watch what happens where to find your logic error...
Thank You. I'll look into it.

Accedi per commentare.

Risposte (1)

Hari Krishna Ravuri
Hari Krishna Ravuri il 7 Nov 2019

0 voti

You may consider debugging your script using MATLAB Debugger. Please refer https://in.mathworks.com/help/matlab/debugging-code.html for more information on debugging a MATLAB script.

Categorie

Scopri di più su Statistics and Machine Learning Toolbox in Centro assistenza e File Exchange

Richiesto:

il 22 Ott 2019

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by