10.

34 Numerical Methods Applied to Chemical Engineering Fall 2015

Homework #0: MATLAB®and Linear Algebra

Please reference the homework guidelines for how to format and submit this assignment.

Problem 1. Cylinder with heat source

Crystalline silicon is annealed to remove defects, before it is sliced to make the substrate for
computer chips and photovoltaic devices. One way to do this is to run an electric current down
the long cylindrical crystal rods. In this problem you will compute the steady-state temperature
distribution inside the rod during this process.
Consider a cylinder of radius R with uniformly distributed heat sources with heat generation
q̇ and temperature-dependent thermal conductivity k(T ). If the cylinder is sufficiently long so
that the temperature can be considered a function of radius only, the differential equation for
temperature can be written as:

d2 T 1 dT q̇
2
+ + =0
dr r dr k
This second-order differential equation can be converted into a system of two coupled first-order
differential equations by defining new variables u1 and u2 as follows:

u1 = T
dT
u2 =
dr
The two differential equations thus obtained are:
du1
= u2
dr
du2 u2 q̇
=− −
dr r k
The surface temperature of the rod T (R) = u1 (r = R) is measured to be 393 Kelvin. If the
system is at steady state, and the rod is so long that axial heat flow can be neglected, the heat flow
to the surface must equal the total heat generated inside the rod, i.e. at the surface,

dT
−2πRLk(393 Kelvin) = πR2 Lq̇
dr
q̇
i.e., u2 (r = R) = −R
2k(393 Kelvin)

Standard ODE solvers expect to integrate starting from zero, not from R. To accomplish this,

1
we can replace r with a new variable y defined this way:

y =R−r
u1 (y = 0) = 393 Kelvin
q̇
u2 (y = 0) = −R
2k(393 Kelvin)
du1
= −u2
dy
du2 u2 q̇
= +
dy R − y k(u1 )

A function you can modify to make your solution username_HW0_P1.m is provided to you, WITH
AN IMPORTANT CATCH. The input and outputs for the function ode23 have been removed. Your
task is to fill in the blanks.
Assume R = 0.20 m and q̇ = 20 kW/m3 . The thermal conductivity of silicon is given by
 −1.3
148 W T
k(T ) = .
m · K 300 K

1. In username_HW0_P1.m, the ode23 function calls an imbedded (pre-written) function dudy

with the entry udash(2) blank. Fill in udash(2) with the expression provided above for dudy .
2

Type help ode23 in the MATLAB command window to better understand the usage and
syntax of the ode23 function. This will be useful in the remaining parts of the problem.

2. Write an additional function in username_HW0_P1.m that plots T and dT /dr versus r using
the subplot function. Include comments at the start of your plotting function by describing
its inputs and outputs briefly. Attach the plot generated from the function.

3. How would you alter the function username_HW0_P1.m if the surface temperature was T =
290 K, and q̇ was 34 kW/m3 ? Copy-paste the original and modified lines of code into your
report. Attach plots of T and dT /dr versus r for the new parameters.

4. How would you have the ode23 function report the temperature in the cylinder at radius =
0.156 meters? Copy-paste the original and modified lines of code into your report.

5. How would you change the absolute or relative error tolerance used by ode23? Copy-paste
the original and modified lines of code into your report.

6. If the function dudy was not imbedded in the function username_HW0_P1.m (i.e., dudy was
stored in a separate m-file, dudy.m), would you have to change anything in the code
username_HW0_P1.m? If so, what?

2
Problem 2. Bessel Approximation with Finite Sums
Bessel functions Jν (r) arise in many problems with cylindrical or spherical symmetry (see for
example Deen’s text, section 4.7).
A series expansion is known for these Bessel functions:
∞
X (−1)k (r/2)ν+2k
Jν (r) =
k!Γ (ν + k + 1)
k=0

where Γ is the Gamma function. For integer ν < 0, the first term of this series is not defined, and
those Bessel functions are computed using a different formula. However, this function works for
both positive and negative non-integer values of ν.

1. Write a MATLAB function my_bessel which takes in nu, r, and a tolerance tol as inputs
and returns an estimate of Jν (r) based on a finite truncated series expansion (i.e., run the
sum for k up to k = kmax rather than k = ∞). The estimate should converge to the true
value of the Bessel function as tol gets very small. If the input nu is a negative integer,
your function should return an error message. Your code should be sufficiently flexible and
vectorized to evaluate for one r value or a vector of r values, returning either a single value or
the corresponding vector of evaluations. For testing your function: it is known that J0 (0) = 1,
J1/2 (r) = sqrt(2/πr) sin r, J−1/2 (r) = sqrt(2/πr) cos r; and J2 (11.6198) ≈ 0. Paste your error
checking and Bessel calculation loop in your solution.

2. Write a function plot_bessel which reads in nu, r, and a vector tolvec and generates a
nicely labeled plot with several traces overlapped:

(a) First trace is a plot of Jν (r) from r = 0 to r = r_max, for tol=tolvec(1)

(b) Second trace is a plot of Jν (r) from r = 0 to r = r_max, for tol=tolvec(2)
(c) Etc.

The last trace is the corresponding plot of Jν (r) generated by the built-in MATLAB function
besselj.m. Please limit the number of traces to 10. Use this function to generate a pretty
plot which shows how your Bessel function routine converges as you tighten the tolerance
using nu = 0 and r_max = 30.

3
Problem 3. The paper “A colloidal quantum dot spectrometer” by Bao and Bawendi (Nature
523:67, 2015) describes a 2-D array of thin films built from different colloidal quantum dots and
used as a filter for a CCD camera such as the one in your cell phone. The device is a micro-
spectrometer that utilizes the the tunable photonic band gap of quantum dots to filter the light
transmitted to pixels of the CCD camera. If each element of the 2-D array of filter films is designated
by a number, i = 1, 2, ..., nf , where nf is the number of films in the array, then the intensity of
light transmitted to the CCD under filter i can be approximated as
nλ
X
Ji = = Ti (λj )Φ(λj ), for each i = 1, 2, . . . , nf . (1)
j=1

where Φ(λ) is the spectrum of light incident on the filter, Ti (λ) is the transmission spectrum of filter
i, and λj with j = 1, 2, ..., nλ is one of a discrete set of wavelengths of light over the range of interest,
say 400-600 nm. nλ is the number of wavelengths at which the incident light will be sampled. This
device aims to use the intensities measured by the CCD under different filter elements to infer the
spectrum of the incident light Φ(λ) in a process called spectral reconstruction.
You will simulate the transmission spectrum of each filter using a simple sigmoidal function
that mimics the properties of a colloidal quantum dot film. For a device having nf filters, filter i
has the simulated transmission spectrum:
1
Ti (λ) =
1+ exp−(λ−si )/w

with λ measured in nanometers, w = 50 nm and si = 400 + 200(i − 1)/(nf − 1) nm. You may use
any built-in MATLAB functions needed during this problem.

1. Plot the transmission spectra of the filters assuming nf = 30 and nλ = 100 with the sampled
wavelengths evenly spaced between 400 and 600 nm.

2. For the special case of a square system with nλ = nf = n, there will be values of Φ :=
[Φ(λ1 ), Φ(λ2 ), · · · , Φ(λn )]T that fit the data exactly. Write a system of equations in the form
AΦ = J where Φ contains the unknown incident spectrum at different discrete wavelengths:
{Φj }nj=1 and J := [J1 , J2 , · · · , Jn ]T contains the integrated intensity measured by the CCD
under different filters. Explicitly write out the form of A.

3. Let the λ samples be evenly spaced between 400 and 600 nm. Using nf = nλ = n =
2, 3, · · · , 20, plot cond(A), kAk, and kA−1 k on the same axes as a function of n. How does
cond(A) compare to kAkkA−1 k for each value of n. Explain any differences.

4. Assume that each of the measured intensities has a relative error no bigger than 0.1%. Derive
an expression giving an upper bound for the relative error in Φ in terms of the relative error
in the measured intensities and any of cond(A), kAk, and kA−1 k. Let the λ samples be
evenly spaced between 400 and 600 nm. Using nf = nλ = n = 2, 3, · · · , 20, compute this
error bound for several values of n. Plot these relative error bounds for Φ as a function of n.

5. Assume that each of the transmission spectra for the filters was measured with a relative
error no bigger than 0.1% at each sampled wavelength. Derive an expression giving an upper
bound for the relative error in Φ in terms of the relative error in the transmission spectra
and any of cond(A), kAk, and kA−1 k. Be sure to note the conditions under which this upper
bound is valid. Let the λ samples be evenly spaced between 400 and 600 nm. Plot the values

4
 of the bound in nf = nλ = n = 2, 3, · · · , 30. For which values of n is the relative error in
Φ bounded from above? Explain what it means for the relative error in Φ to have no upper
bound.

6. Consider an alternative case in which nf > nλ for which there does not exist an exact solution
for Φ. A linear least squares solution can be constructed from the system of equations:
AT AΦ = AT J, however. We wish to estimate the sensitivity of the least squares solution to
perturbations in the measured intensities, J. Assume that the relative error in J is no bigger
than 0.1%. Derive an upper bound on the relative error in the least squares solution. Is this
bound independent of the measured values J? You may use the fact that the matrix A has
full column rank when nf > nλ and the Ti (λ) are linearly independent so that AT A has an
inverse. Explain what factors will enhance the sensitivity of the least squares solution.

7. Explain the intuitive meaning of the ill-conditioning for these systems of equations. Examine
extended data figure 1 of the cited Nature article and compare the simulated transmission
spectra to those of the experimental colloidal quantum dot thin films. Suggest properties for
the basis functions (transmission spectra) needed for accurate calculation of the incident light
spectrum in the range of wavelengths: 400-600 nm.

5
MIT OpenCourseWare
https://ocw.mit.edu

10.34 Numerical Methods Applied to Chemical Engineering

Fall 2015

For information about citing these materials or our Terms of Use, visit: https://ocw.mit.edu/terms.