complex numbers as function output

3 vues (au cours des 30 derniers jours)
Sunetra CV
Sunetra CV le 22 Oct 2019
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 commentaires
dpb
dpb le 22 Oct 2019
I'd guess z goes negative...use debugger to step through and watch what happens where to find your logic error...
Sunetra CV
Sunetra CV le 23 Oct 2019
Thank You. I'll look into it.

Connectez-vous pour commenter.

Réponses (1)

Hari Krishna Ravuri
Hari Krishna Ravuri le 7 Nov 2019
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.

Community Treasure Hunt

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

Start Hunting!

Translated by