Help with trapz(): Getting complex values

Here is my code: I am getting complex values for "T_integral" and that is messing up my temperature profile TSil. Any help will be appreciated.

 %%Initializing Variables
Lz = 0.00038; % thickness of the silicon wafer
Ft = 1; % simulation time
Nz = 10; % number of grid points
Nt = 100; % number of time steps
dz = Lz/Nz; % grid cell size
dt = Ft/Nt; % time step size
z = linspace(0,Lz,Nz+1);

 %%Thermal Properties
alpha_Sil =  0.00008; % alpha for silicon (m^2/sec)
k_Sil = 149;  % thermal conductivity for silicon (W/m K)
d = alpha_Sil*dt/(dz^2);
Ta = 298;   % temperature of ambient air (K)
ha = 15;    % heat transfer coefficient of air at Ta (W/m^2 K)

 %%Optical Properties
a_Sil = 52.6; % absorption coefficient (m^-1)
r_SilAir = 0.34; % refelctivity of the Si-air interface (experimentally measured)
t_Sil = 0.9802;  % transmissivity of the silicon layer due to absorption alone
e_Sil = 1 - r_SilAir;   % emissivity of silicon surface

 %%Constants
sigma = 5.67037E-08;    % Stefan-Boltzmann constant
F_12 = 0.017517;    % energy emitted by a black body contained within the two wavelength intervals of the IR camera
Ec = 10.3439;   % energy measured by the IR camera

 %%Finite Difference Scheme (Crank Nicholson in space & F.E in time)
TSil = 329*ones(1,Nz+1);    % initial condition (arbitrary temperature profile)
TSilnew = TSil;
t = 0;
for kt = 1:Nt+1
    t = t + dt;
    for kz = 1:Nz+1
        Y(kz) = (TSil(1,kz)^4).*exp(-a_Sil*z(kz));
    end
    T_integral = trapz(z,Y(:));
    TSilnew(1) = (((Ec - (r_SilAir*F_12*sigma*(Ta^4)) - (e_Sil*a_Sil*sigma*F_12*T_integral))/(e_Sil*t_Sil*F_12*sigma)).^0.25);
    for kz = 2:Nz
        TSilnew(kz) = (d/(2*(1+d)))*(TSilnew(kz-1) + TSilnew(kz+1) + TSil(kz-1) + TSil(kz+1)) + ((1-d)/(1+d))*TSil(kz);
    end
    TSilnew(kz+1) = (d/(1+d-(dz*d*ha/k_Sil)))*(TSilnew(Nz) + TSil(Nz)) + ((1-d-(dz*d*ha/k_Sil))/(1+d-(dz*d*ha/k_Sil)))*TSil(Nz+1) - ((2*dz*ha*d/k_Sil)/(1+d-(dz*d*ha/k_Sil)))*Ta;
    TSil = TSilnew;
end

Answer 1

[per isakson]

• Put the code in a function
• Set a conditional, not(isreal(Y)), breakpoint on the line

     T_integral = trapz(z,Y(:));

• Run the function and you will see that eventually Y is complex
• Set a conditional, not(isreal(TSilnew)), breakpoint on the line

    TSilnew(1) = (((Ec - (r_SilAir*F_12*sigma*(Ta^4)).....

• Run the function and you'll see that TSilnew(1) is set to a complex value for kt=56. You have a negative value raised to 1/4.

Answer 2

[Shashank]

Per, How do I make that value real then?! And why would a simple integration result in a complex value at all?

[Shashank]

Thanks for your help btw.

[per isakson]

It isn't the integral that produce the the complex value in the first place. The integrand is complex because

    TSilnew(1) = (((Ec - ....

returns a complex value. You have a negative value raised to 1/4. Thus, there is an error in your code.

[Shashank]

I realized that. Thanks for helping me out Per.