Homework 9

November 13, 2006

Problem1
1. Let us first derive an analytical espression for Smodel (λ, t, θ) using both
the models mentioned in the problem statement.

• Model 1: Sequential Reaction Model

According this model C is formed from B through an intermediate
species X, we can write the rate of formation of B, X and C as
shown in Equation(1). In writing the rate rules in Equation(1) we
have made use of the fact that the volume is constant.
d[B]
= −k1 [B]
dt
d[X]
= k1 [B] − k2 [X]
dt
d[C]
= k2 [X] (1)
dt
We can solve the equations simultaneously to get the expressions
for the concentrations of [B], [X] and [C] shown in Equation(2).

[B] = [B]0 e−k1 t

� �
e−k1 t e−k2 t
[X] = [B]0 k1 +
k2 − k1 k1 − k2
� �
k2 e−k1 t k1 e−k2 t
[C] = [B]0 1 + + (2)
k1 − k2 k2 − k1

If we look at the absoption spectra of B and C in Figure 1 as

given by SVD we realise that B has two distinct peaks and C has

Cite as: William Green, Jr., course materials for 10.34 Numerical Methods Applied to Chemical Engineering, Fall 2006.

MIT OpenCourseWare ([Link] Massachusetts Institute of Technology. Downloaded on [DD Month YYYY].

 The absoption spectrum of species B and C

0.6
SVD−absorption of B
SVD−absorption of C
0.5

The normalized absorption

0.4

0.3

0.2

0.1

−0.1
430 440 450 460 470 480 490

The wavelength (nm)

Figure 1: Normalized Absoption Spectrum of B and C obtained from SVD

only one peak, which suggests that B can be written as a sum of

two different gaussians and C can be written as a sum of only one
gaussian as shown in Equation(3).

� � � �
(λ − λB,1 )2 (λ − λB,2 )2
AB (λ) = AB,1 exp − 2
+ AB,2 exp − 2

wB,1 wB,2

� �
(λ − λC,1 )2
AC (λ) = AC,1 exp − 2
(3)
wC,1

Putting everything together, we can get an analytical expression

for Smodel (λ, t, θ) as shown in Equation(4).

� � � � ��
−k1 t (λ − λB,1 )2 (λ − λB,2 )2
Smodel (λ, t, θ) = [B]0 e AB,1 exp − 2
+ AB,2 exp − 2
wB,1 wB,2
� �� � ��
k2 e−k1 t k1 e−k2 t (λ − λC,1 )2
+[B]0 1+ + AC,1 exp − 2
(4)
k1 − k2 k2 − k1 wC,1

As is mentioned in the problem statement we cannnot determine

[B]0 or the magnitude of absoptions independently and we have

Cite as: William Green, Jr., course materials for 10.34 Numerical Methods Applied to Chemical Engineering, Fall 2006.

MIT OpenCourseWare ([Link] Massachusetts Institute of Technology. Downloaded on [DD Month YYYY].

 to club some of the parameters together. If we do that we get

Equation(5)

� � � � ��
−k1 t (λ − λB,1 )2 (λ − λB,2 )2
Smodel (λ, t, θ) = e KB,1 exp − 2
+ KB,2 exp − 2
wB,1 wB,2
� �� � ��
k2 e−k1 t k1 e−k2 t (λ − λC,1 )2
+ 1+ + KC,1 exp − 2
(5)
k1 − k2 k2 − k1 wC,1
The parameters that we need to determined in the above model

are k1 , k2 , KB,1 , KB,2 , KC,1 , λB,1 , λB,2 , λC,1 , wB,1 , wB,2 and wC,1 .

• Model 2: Independent Reaction Model

In this reaction model the rate of change of B and C are indepen­

dent of each other. The final expressions for the concentrations of

B and C at any given time t is given in Equation(6).

[B] = [B]0 e−k1 t

� �
[C] = [Y ]0 1 − e−k3 t (6)
The absoption spectrum will be assumed to be similar to the

one provided in Equation(3) and the analytical expression for

Smodel (λ, t, θ) can be written as shown in Equation(7).

� � � � ��
(λ − λB,1 )2 (λ − λB,2 )2
Smodel (λ, t, θ) = [B]0 e−k1 t AB,1 exp − 2
+ AB,2 exp − 2
wB,1 wB,2
� � ��
� � (λ − λC,1 )2
−k3 t
+[Y ]0 1 − e AC,1 exp − 2
(7)
wC,1
Again we can rewrite the model above in terms of the least num­

ber of determinable parameters as shown in Equation(8) and the

parameters that we need to determine are k1 , k3 , KB,1 , KB,2 , KC,1 ,

λB,1 , λB,2 , λC,1 , wB,1 , wB,2 and wC,1 .

� � � � ��
(λ − λB,1 )2 (λ − λB,2 )2
Smodel (λ, t, θ) = e−k1 t KB,1 exp − 2
+ KB,2 exp − 2
wB,1 wB,2
� � ��
� � (λ − λC,1 )2
−k3 t
+ 1 − e KC,1 exp − 2
(8)
wC,1
2. The estimate of noise level given by SVD is 4.857 × 10−4 . This noise
level is calculated by summing up the signal corresponding to all the
small eigenvalues. Thus the noise matrix is given as
�
S= σi xi aTi
i=3

Cite as: William Green, Jr., course materials for 10.34 Numerical Methods Applied to Chemical Engineering, Fall 2006.

MIT OpenCourseWare ([Link] Massachusetts Institute of Technology. Downloaded on [DD Month YYYY].

 where xi and ai are respectively the concentration and absoptions of

the species. We are told that the noise at each wavelength and each
time comes from the same normal distribution. Thus we can calculate
the standard deviation of the noise using the values obtained from the
noise matix.
Since we dont know anything about the reactions taking place and the
abosorption spectrum of the species we will just assume that all the
parameters are uniformly distributed in some arbitrarily large range
of values chosen. For example the values of the rate constant will be
positive and have values less than 1015 which is the range of the rate
constant for a unimolecular reaction. Also since we assume that all the
parameters are independent of each other
�
n
p(Θ) = pi (θi )
i=1
.

3. • Model 1: Sequential Reaction Model

At t=0 we are told that the concentration of C is 0. Thus all
the absoption signal that we observe is due to species B. Using
Equation(5) we realise that the absoption at t=0 can be used to
determine the value of KB,1 , KB,2 , λB,1 , λB,2 , wB,1 and wB,2 . The
matlab code to do the above fitting is provided in file fitAbsForB.m
and the parameters that it calculates are KB,1 = 0.3000, KB,2 =
0.0501, λB,1 = 450.6204, λB,1 = 469.2274,wB,1 = 3.1041 and
wB,2 = 9.7590. The comparison of the Best Fit Absorption and
the Signal Strength at t=0 are provided in Figure 2.
Also if compare the absoption spectrum of B and C using Figure
1, we realise that the absorption signal at wavelength 470 nm
is solely due to B. Thus we can calculate the rate constant k1
using the absoption signal at 470 nm and using the parameters
KB,1 , KB,2 , λB,1 , λB,2 , wB,1 and wB,2 which we have calculated
previously. The matlab code bestfitk1.m is used to get a best
fit value of k1 . The best fit value of k1 = 0.0497 and the graphical
comparison of the model vs data is shown in Figure 3.
From Figure 3 it is apparent that the concentration of B is close to
zero at time t=100 s. Thus the at t=100, all the absoption is due to
species C. We dont know the concentration of C and so we will be
able to calculate the absoption profile of C to within a constant fac­
tor. This fitting is performed using the matlab file fitAbsForC.m.

Cite as: William Green, Jr., course materials for 10.34 Numerical Methods Applied to Chemical Engineering, Fall 2006.
MIT OpenCourseWare ([Link] Massachusetts Institute of Technology. Downloaded on [DD Month YYYY].
 The absoption spectrum of species B
0.35
Best Fit Absorption of B
Signal at t=0
0.3

0.25

0.2
The Absorption

0.15

0.1

0.05

−0.05
430 440 450 460 470 480 490

The wavelength (nm)

Figure 2: BestFit Absoptions of B compared to the Signal Strength at t=0

The normalized concentration of B

1.2
Best Fit [B]
[B] obtained due to signal at 470 nm
1

0.8
Concentration

0.6

0.4

0.2

−0.2
0 20 40 60 80 100

time (sec)

Figure 3: Normalized Absoption Spectrum of B

Cite as: William Green, Jr., course materials for 10.34 Numerical Methods Applied to Chemical Engineering, Fall 2006.

MIT OpenCourseWare ([Link] Massachusetts Institute of Technology. Downloaded on [DD Month YYYY].

 The parameter values are KC,1 = 0.7339, λC,1 = 439.844 and

wC,1 = 6.94. The comaprison of fitted value of absoption of C and
signal strength at t=100 s is given in Figure 4.

The absoption spectrum of species C

0.8
Best Fit Absorption of C
0.7 Signal at t=100 s

0.6

0.5
The Absorption

0.4

0.3

0.2

0.1

−0.1
430 440 450 460 470 480 490
The wavelength (nm)

Figure 4: BestFit Absoptions of C compared to the Signal Strength at t=101

Finally we see that C has a absorption maximum at 439.844 nm

and B does not absorb at that concentration. Thus we can use
the signal strength at 439.844 nm to calculate the best fit value of
k2 . This is done in matalb file bestfitk2 model1.m. The best fit
value of k2 = 0.0889 and the graphical comparison of normalized
concentration of B with the signal at 440 nm is given in Figure 5.
We can now use all the values, estimated so far, as the initial guess
for a routine which tries to fit the entire matrix of signal. The mat-
lab function that accomplishes this task is p3 model1 error bestfit.m.
The deviation between model prediction and the data and the er­
ror is plotted in Figure 6. The best fit values of the parameters is
given in Table 1.

Cite as: William Green, Jr., course materials for 10.34 Numerical Methods Applied to Chemical Engineering, Fall 2006.

MIT OpenCourseWare ([Link] Massachusetts Institute of Technology. Downloaded on [DD Month YYYY].

 The normalized concentration of C

1.2
Best Fit [C]
[C] obtained due to signal at 440 nm
1

0.8
Concentration

0.6

0.4

0.2

−0.2
0 20 40 60 80 100
time (sec)

Figure 5: Normalized Absoption Spectrum of C

Table 1: Best Fit Values for Model 1

k1 0.0286

k2 0.8123

KB,1 0.2259

KB,2 0.0369

KC,1 0.7733

λB,1 450.53

λB,2 469.23

λC,1 439.79

wB,1 3.2342

wB,2 9.7594

wC,1 6.7872

Cite as: William Green, Jr., course materials for 10.34 Numerical Methods Applied to Chemical Engineering, Fall 2006.

MIT OpenCourseWare ([Link] Massachusetts Institute of Technology. Downloaded on [DD Month YYYY].

 Distribution of error due to Model−2

6000

Expt. Error distribution

Exact normal distribution

5000

Frequency of the Error

4000

3000

2000

1000

0
−0.04 −0.02 0 0.02 0.04 0.06 0.08
Magnitude of Error

Figure 6: Error Distribution due to Model 1

• Model 2: Independent Reaction Model

For this model the initial guesses for all the parameters besides
k3 is calculated in exactly the same way as the previous Model.
For calculating the value of k3 , we take the wavelength 440 nm
and fit the concentration of C to the analytical formula given
in Equation(6). This fit is performed using the matlab function
given in command bestfitk2 model2.m. The comparison of the
fit with the Signal strength at 440 nm is given in Figure 7 and the
value of k3 = 0.0303. Again we use the values, calculated so far,
as the initial guesses and perform a full fledged optimization to
calculate the values of all the parameters using the matlab func­
tion p3 model2 error bestfit.m. The values of the parameters
calculated are given in Table 2. The distribution of the error is
shown in Figure 8.

Cite as: William Green, Jr., course materials for 10.34 Numerical Methods Applied to Chemical Engineering, Fall 2006.

MIT OpenCourseWare ([Link] Massachusetts Institute of Technology. Downloaded on [DD Month YYYY].

 The normalized concentration of C

1.2
Best Fit [C]
[C] obtained due to signal at 440 nm
1

0.8
Concentration

0.6

0.4

0.2

−0.2
0 20 40 60 80 100
time (sec)

Figure 7: Normalized Absoption Spectrum of C

Table 2: Best Fit Values for Model 2

k1 0.0500
k3 0.0250
KB,1 0.3000
KB,2 0.0501
KC,1 0.7999
λB,1 450.62
λB,2 469.24
λC,1 439.84
wB,1 3.0993
wB,2 9.7915
wC,1 6.9303

Cite as: William Green, Jr., course materials for 10.34 Numerical Methods Applied to Chemical Engineering, Fall 2006.

MIT OpenCourseWare ([Link] Massachusetts Institute of Technology. Downloaded on [DD Month YYYY].

 Distribution of error due to Model−2

400
Expt. Error distribution
Exact normal distribution
350

300

Frequency of the Error 250

200

150

100

50

0
−2 −1.5 −1 −0.5 0 0.5 1 1.5 2
Magnitude of Error −3
x 10

Figure 8: Error Distribution due to Model 2

4. To apply the χ2 test we assume that the underlying distribution of

errors is a normal distribution. Then we calculate the χ2 value of the
distribution using the Equation(9)
� δi2
χ2 = (9)
i σ

In the above equation δi are the individual errors and σ is the standard
deviation of the error around 0 that we have calculated in part 2 of
the problem. The value χ2 will be distributed according to the χ2 ­
distribution if the number of sample points in large. In our case the
number of sample points are large (5151 points). The number of fitted
parameters are 11 and thus the degree of freedom of the χ2 -distribution
is 5140 (=5151-11).

• Model 1: Sequential Reaction Model

The χ2 value of model 1 = 1.3 × 106. The p-value is 0. Thus there
is 100% probability that the model is incorrect.
• Model 1: Sequential Reaction Model
The χ2 value of model 1 = 5.480 × 103 . The p-value is 0.0005.
Thus there is 99.95% probability that the model is incorrect.

10

Cite as: William Green, Jr., course materials for 10.34 Numerical Methods Applied to Chemical Engineering, Fall 2006.

MIT OpenCourseWare ([Link] Massachusetts Institute of Technology. Downloaded on [DD Month YYYY].

 The values of χ2 of model 1 ad model 2 are calculated using the func­

tions p3 model1 error bestfit.m and p3 model2 error bestfit.m
respectively.

5. For Model 2 we have made an array of k1 and k3 and calculated χ2

values for all the combinations of k1 and k3 . The graph is plotted using
the matlab function p5 contourplots.m and the graph is shown in
Figure 9.

5515
0.025010

0.025008 5510

0.025006 5505

0.025004 5500
k3

0.025002 5495

0.025000 5490

0.024998 5485

0.024996
0.04996 0.04998 0.05000 0.05002 0.05004 0.05006 0.05008
k1

Figure 9: χ2 with different values of kf and kr

To convert the χ2 plot to a plot of probability we just calculate the

value of Δχ2 = χ2 (k1 , k3 ) − χ2min and this value of Δχ2 is distributed
with a χ2 distribution of 2 degrees of freedom. The contour plot of 95%
and 90% confidence interval of k1 and k2 are plotted in Figure 10.
From Figure 10 we can read off the 95% confidence interval for k1 and
k2 as 0.04998 ≤ k1 ≤ 0.050065 and 0.024995 ≤ k3 ≤ 0.025008. Both
the figures are generated using the matlab code p5 countourplots.m.
We can also use the built-in matlab functions nlinfit and nlparci
to calculate the confidence intervals on k1 and k2 and the confidence
intervals given by these functions are 0.049968 ≤ k1 ≤ 0.050074 and
0.024993 ≤ k3 ≤ 0.025014. The matlab function which generates these
confidence inervals is p5 model2 param est.m. The confidence inter­
vals provided by the matlab function are very close to the values given

11

Cite as: William Green, Jr., course materials for 10.34 Numerical Methods Applied to Chemical Engineering, Fall 2006.
MIT OpenCourseWare ([Link] Massachusetts Institute of Technology. Downloaded on [DD Month YYYY].
 0.025010

0.025008 548
698 5484.7 6.1698
6.1 835
0.025006 548 35
8
.7
84
54

5486 835
8
0.025004

.169
k3

.7
5484
54
0.025002

84
.78
54

35
86
0.025000 .16
98 5
5484.783.1698
5486
0.024998

0.024996
0.04996 0.04998 0.05000 0.05002 0.05004 0.05006 0.05008
k1

Figure 10: 90% and 95% confidence contours

by plotting the graph. Due to the large number of data points (5151
data points) and relatively fewer fitted parameters, the confidence inter­
vals of all the parameters are small. The matlab functions approximate
the bottom of the well of the χ2 distribution as a parabola. If the con­
fidence intervals are small then the parabola approximation is usually
very good.

6. In part 2 of the problem we assumed that the distribution of θ was

uniform. Specifically for k1 and k2 we assumed the probability distri­
bution was uniform between 0 and 1014 . If one were to normalize the
probability distribution the constant multiplying factor of pprior (k1 , k2 )
will get eliminated. Thus for a uniform prior probability distribution
the answer to this part will look exactly like Figure 10. The confidence
intervals calculated will also remain unchanged. The values of k1 and
k2 will be effected only if we have specific prior information about the
probability distribution of k1 and k2 .

12

Cite as: William Green, Jr., course materials for 10.34 Numerical Methods Applied to Chemical Engineering, Fall 2006.

MIT OpenCourseWare ([Link] Massachusetts Institute of Technology. Downloaded on [DD Month YYYY].