MATLAB Answers

# Help with trapz(): Getting complex values

7 views (last 30 days)
Shashank on 21 Dec 2012
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
##### 0 CommentsShowHide -1 older comments

Sign in to comment.

### Answers (2)

per isakson on 21 Dec 2012
Edited: per isakson on 21 Dec 2012
• 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.
##### 0 CommentsShowHide -1 older comments

Sign in to comment.

Shashank on 21 Dec 2012
Per, How do I make that value real then?! And why would a simple integration result in a complex value at all?
##### 3 CommentsShowHide 2 older comments
Shashank on 21 Dec 2012
I realized that. Thanks for helping me out Per.

Sign in to comment.

### Community Treasure Hunt

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

Start Hunting!