x=1/a
Solution:
We apply Newton's method to find the root to f(x)=0
Let f(x) = 0 = x-1/a
xk+1 = xk - f(xk)/f'(xk)
= xk - (xk-1/a)/1 = 1/a
We ended up where we had started -- bad luck. When we devise an
iterative scheme blindly, we have a 50% chance
convergence. Thus, we express f(x) a bit differently and hope we
get lucky this time.
f(x) = 0 = 1/x - a where a is given
xk+1 = xk - f(xk)/f'(xk)
= xk - (1/xk-a)/(-xk2)
= xk + xk - a*xk2
= xk(2-a*xk)
Start with x=1 or some other reasonable approximation.
Iterate the above formula until convergence to the specified accuracy.
pH indicator ------Absorbance----- value conc. A400 A401 ... A700 -------------------------------------- 1.5 0.1 : : : 2.3 0.5 : : : : : : : : : : : : :Find the regression model from the above data. Given the spectrum of an unknown sample, the regression model should return an estimate of its pH value. If you use vector-matrix notation, clearly indicate the dimensions and what the rows or columns are.
Solution:
See Homework #5, Problem 2.
Let X be mean-centered variance-scaled absorbance spectra
of n calibration samples (dimension of X: n´301); in essense, X is the data under
"Absorbance" in the above table. Let Y be a mean-centered
variance-scaled output matrix containing pH and concentration of
pH indicator for n calibration samples (dimension of Y: n´2). We apply regression to relate X
to Y.
Note that the following naive normal equation from linear
regression does not work because XT*X is singular (for
n<301) or nearly singular (for n>301) and cannot be
inverted.
Step 1. Find the eigenvalues for the square input covariance
matrix XT*X (which is the matrix that gives us trouble
with the naive regression approach) and list them in decreasing
order (li, i=1,2,...).
Find the correponding orthogonal eigenvectors (Vi).
The first eigenvector corresponding to the largest eigenvalue of
XT*X is the direction along which there is the most
variance in the calibration spectra; and so on.
Step 2. Represent X in terms of the new coordinate system given
by the eigenvectors.
Y=X*a
where a=(XT*X)-1*XT*Y ... No good!
Instead, we have to trim the number of factors. The following is
one such approach utilizing concepts from eigenvalue-eigenvector.
vi=X*Vi i=1,2,...,m
Step 3. Express Y as a linear combination of X, which is
now expressed in the new eigenvector coordinate.
Y = a1*v1 + a2*v2 + a3*v3 + ... + am*vm
where the coefficients are given by the projection formula:
ai = (Y,vi)/(vi,vi)
ai = (Y,vi) if vi is normalized, i.e., (vi,vi)=1
Step 4. Finally, the regression equation Y(x) predicts Y (which
is composed of pH and pH indicator concentration, dimension of
Y=1´2) as a function of an input
spectrum x (dimension of x: 1´301).
Y=x*a
The number of terms to keep in the model (m) is that when the sum
of squared error (sse) on the calibration sample turns a corner
as the number of terms i is increased. In general, m<n,
m<301. Since we are interested in only the pH value, not the
pH indicator concentration, we perform regression as decribed
above but utilize only the pH portion of the output while
disgarding the pH indicator concentration portion.
Alternatively, we can perform regression on just the pH portion
without the pH indicator concentration.
[4*x1 ]
A*x = [x1+4*x2+4*x3 ]
[2*x1+4*x3 ]
[x1]
where x = [x2]
[x3]
Solution:
See Homework #6, Problem 3.
Note that matrix A is just a collection of numbers; thus, strictly, it is
inappropriate to speak of the eivenvalue of a matrix A, which is different from
the eigenvalue of a linear transformation A defined by the operator in the
problem statement.
[4 0 0 ][x1]
A*x = [1 4 4 ][x2]
[2 0 4 ][x3]
Eigenvalues-Eigenvectors
A*x=l*x
(l*I-A)*x=0
eigenvalues
det|l*I-A)|=0
(l-4)3=0
l1=4, l2=4, l3=4 (repeated 3 times)
eigenvectors
[ 0 0 0 ][x1]
(l*I-A)*x=[ -1 0 -4 ][x2] = 0
[ -2 0 0 ][x3]
-x1-4*x3=0 and -2*x1=0
Thus, x1=x3=0 & x2=arbitrary -->normalize--> x2=1
In this problem, there is only one eigenvector. x=[ 0 1 0
]T
Solution:
With only one eigenvector. We cannot even speak of
orthogonality (unless we have more than one vector). In
general, the eigenvectors are not orthogonal (unless the linear
transformation is a special one, e.g., symmetric).
Solution:
To find exp(A*t), we go through the Jordan normal form J.
Because there is only one eigenvector, the Jordan normal form J
has only one Jordan block.
The problem reduces to finding a similarity transform T such that
A*T=T*J.
A*T = T*J
T has three columns. T=[T1 T2 T3]
A*[T1 T2 T3] = [T1 T2 T3]*J
[ l 1 0 ]
A*[T1 T2 T3] = [T1 T2 T3]*[ 0 l 1 ]
[ 0 0 l ]
[A*T1 A*T2 A*T3] = [lT1 T1+lT2 T2+lT3 ]
[ 4 0 0 ][ 0 T12 T13 ] [ 0 T12 T13 ][ 4 1 0 ]
A*T=[ 1 4 4 ][ 1 T22 T23 ]=[ 1 T22 T23 ][ 0 4 1 ]
[ 2 0 4 ][ 0 T32 T33 ] [ 0 T32 T33 ][ 0 0 4 ]
1st column: A*T1=lT1 --> T1 is the eigenvector.
T1=v1=[0 1 0]T.
2nd column: A*T2=lT2+T1
4*T12 =0+4*T12
1*T12+4*T22+4*T32=1+4*T22 --> 1*T12+4*T32=1 --> T32=0.25
2*T12 +4*T32=0+4*T32 --> 2*T12=0 --> T12=0
Thus, T2=[0 0 0.25]T
Note that T2 is not necessarily orthogonal to T1.
T2 should not be forced to be orthogonal to T1; nor should T2 be forced to be normalized.
3rd column: A*T3=lT3+T2
4*T13 =0 +4*T13
1*T13+4*T23+4*T33=0 +4*T23 --> 1*T13+4*T33=0 --> T33=-T13/4=-0.03125
2*T13 +4*T33=0.25+4*T33 --> 2*T13=0.25 --> T13=0.125
Thus, T3=[0.125 0 0.03125]T
Note that T3 is not (and should not be forced to be) orthogonal to T1 and T2.
Summary (combine all 3 columns of T)
[ 0 0 0.125 ]
T = [ 1 0 0 ]
[ 0 0.25 -0.03125 ]
Note that T is not orthogonal.
Form exp(J*t)
[ exp(4*t) t*exp(4*t) t*t*exp(4*t)/2 ]
exp(J*t)=[ 0 exp(4*t) t*exp(4*t) ]
[ 0 0 0 ]
Form exp(A*t)
A=T*J*T-1
exp(A*t)=T*exp(J*t)*T-1
[ 0.03125*exp(4*t) 0 0 ]
=[ 0.03125*t*exp(4*t)+0.125*t*t*exp(4*t) 0.03125*exp(4*t) 0.125*t*exp(4*t) ]
[ 0.0625*t*exp(4*t) 0 0.03125*exp(4*t) ]
dx/dt = A*x x(0)=x0
Solution:
x(t)=exp(A*t)*x0
=T*exp(J*t)*T-1*x0
ót
g(t) = ô K(t-t)*f(t)*dt
õ0
Solution:
This is a linear transform that changes an input vector f into
an output vector g. Show that the following two conditions
of being a linear transform is satisfied.
ót
g(t) = A(f) = ô K(t-t)*f(t)*dt
õ0
Condition 1 for a linear transform (addition of two vectors)
A(f+h)=A(f)+A(h)
ót
A(f+h) = ô K(t-t)*[f(t)+h(t)]*dt
õ0
ót
= ô K(t-t)*f(t)*dt
õ0
ót
+ ô K(t-t)*h(t)*dt
õ0
= A(f) + A(h)
Condition 2 for a linear transform (multiplication by a scalar a)
A(a*f)=a*A(f)
ót
A(a*f) = ô K(t-t)*[a*f(t)]*dt
õ0
ót
= a*ô K(t-t)*f(t)*dt
õ0
= a*A(f)
Do not confuse a linear transformation A with a linear function f.
Condition 1 for a linear function (addition of two arguments)
f(s+t)=f(s)+f(t)
Condition 2 for a linear function (multiplication by a constant)
f(a*t)=a*f(t)
To simplify the math we can extend the lower limit from 0 to
-¥ if we set f(t)=0 for t£0 (i.e., outside the original range).
Eigenvalues-Eigenvectors
A(f) = l*f
ót
A(f) = ô exp(-(t-t)/T)/T*f(t)*dt = l*f(t)
õ-¥
ót
ô exp(t/T)*f(t)*dt = l*T*exp(t/T)*f(t)
õ-¥
ót
ô ef(t)*dt = l*T*ef(t) where ef(t)=exp(t/T)*f(t)
õ-¥
As always, f(t)=0 is a trivial case. We seek a function
ef(t) whose integral (LHS) is the same as the function
itself (RHS). Equivalently, we seek a function whose
derivative is identical to to the function itself. An
exponential function satisfies this criterion.
ef(t)=exp(a*t) a>0 when the lower integration limit is -¥
1/a=l*T --> l=1/(a*T)
Eigenvalue: l=1 a=1/T
Eigenvector: any constant. f=1
Since any multiple of an eigenvector is also an eigenvector, we
simply set f=1 to be the eigenvector.
In fact, there are an infinite number of eigenvalues and eigenvectors, as shown below.
Eigenvalue: l=2 a=1/(2T)
Eigenvector: ef(t)=exp(a*t)=exp(t/2/T)=exp(t/T)*f(t) -->
f(t)=exp(-t/2/T)
Eigenvalue: l=3
Eigenvector: f(t)=exp(-2*t/3/T)
In general:
Eigenvalue: l=any positive real number
Eigenvector: ef(t)=exp(a*t)=exp(t/l/T)=exp(t/T)*f(t) -->
f(t)=exp((1/l-1)*t/T)
Solution:
A symmetric transformation satisfies the following: (A*x,y)=(x,A*y).
Since we deal with a scalar product, we must define it first.
Thus, the answer depends on how we define it. Let us consider a
scalar product that is defined as an integral, such as the one
below.
ó¥
(f,h) = ô exp(-t/T)*f(t)*h(t)*dt
õ0
Since the eigenvectors are exponential functions of different
decay constants, the scalar product between two different
exponential functions is basically the integral of exponential
functions. The integral is never 0; the scalar product of
eigenvectors is not 0. Thus, the eigenvectors are not
orthogonal, and we conclude that this is not a symmetric
transformation.
Solution:
Eigenvalue-eigenvector #1: A*x1=l1*x1
Eigenvalue-eigenvector #2: A*x2=l2*x2
Symmetric transformation: (A*x1,x2)=(x1,A*x2)
Proof:
(A*x1,x2)=(l1*x1,x2)=l1*(x1,x2)
(x1,A*x2)=(x1,21*x2)=l2*(x1,x2)
Subtract the above two equations one from the other:
(A*x1,x2)-(x1,A*x2) = 0 = (l1-l2)*(x1,x2)
For distinct eigenvalues, (l1-l2)¹0,
This forces (x1,x2)=0. Thus, x1 and x2 are orthogonal.
ds/dt=J(s)-v(s)
where s is the reactant concentration
J(s) is the mass transfer rate: J(s)=kL(sb-s)
v(s) is the reaction rate: v(s)=vm*s/(K+s+Ki*s2)
Solution:
At steady-state, ds/dt=0. We have J(sss)=v(sss).
The above is a cubic equation in s. In general, there are
three roots for a polynomial of degree 3. Since we deal with a
real system, only real roots are of interest to us: either 1, 2,
or 3 distinct real roots are possible. --> 1, 2, or 3
steady-states.
J=v
kL(sb-s)=vm*s/(K+s+Ki*s2)
kL(sb-s)*(K+s+Ki*s2)-vm*s=0
Solution:
Perform linearlized stability analysis around each
steady-state point.
Deviation variable: S=s-sss
dS/dt=(J'(sss)-v'(sss))*S
Eigenvalue: l=J'(sss)-v'(sss)
Stable if the above eigenvalue is negative: J'(sss)-v'(sss)<0
We do not need to worry about the "real" part of the eigenvalue
being negative, because complex eigenvalues must appear in
conjugate pairs -- and we have only one eigenvalue for this
one-dimensional system. Consequently, there is only
exponential approach (no oscillatory approach) to a steady-state.
Even without rigorous stability analysis, when there is only one
steady-state for a one-dimensional real system (not an
artificially made-up mathematical one), the lone steady-state
must be stable. Otherwise, the state diverges to either +¥ or -¥, which is strictly a mathematical
concept, as all real physical systems are bound. Likewise, when
there are three steady-states for a one-dimensional real system,
the middle steady-state must be unstable and the two end
steady-states must be stable.
Solution:
Stability criterion: (slope of J) < (slope of v)
Solution:
The plot is S-shaped.
sss
^
| **************
| **
| *
| *
| *
| **
| **
| *
| *
| *
| **
| **
|**********
+--------------------------->
kL
d2y --- = f2*f(y) Boundary Conditions: at x=0 dy/dx=0; at x=1 y=1 dx2where f is the Thiele modulus. Unless specified otherwise, use the following saturation kinetics:
f(y)=a*y/(b+y)
Solution:
Let y=exp(l*x)
Characteristic equation: l2-f2*k = 0
Roots: l1=-f*sqrt(k), l2=f*sqrt(k)
Solution: y=A*exp(-f*sqrt(k)*x) + B*exp(f*sqrt(k)*x)
Apply BC at x=0: y'=0=-A*f*sqrt(k) + B*f*sqrt(k) --> A=B
Apply BC at x=1: y =1= A*exp(-f*sqrt(k)) + B*exp(f*sqrt(k))
--> A=B=1/(exp(-f*sqrt(k))+exp(f*sqrt(k)))=1/2/cosh(f*sqrt(k))
y=cosh(f*sqrt(k)*x)/cosh(f*sqrt(k))
Solution:
Express the given 2nd-order ODE as two 1st-order ODEs.
dy/dx=z
dz/dx=f2*f(y)
Integrate from x=0 to x=0.1 with Euler's method.
IC: y(0)=0.5 (... guess this)
z(0)=0
Euler's method
y(0.1)=y(0)+z(0)*step =0.5+0*0.1
z(0.1)=z(0)+f(y(0))*step=0 +f(0.5)*0.1
Guess different values of y(0) so as to satisfy y(1)=1.
dy/dx ==> 0.5(yi+1-yi-1)/Dx d2y/dx2 ==> (yi+1-2yi+yi-1)/(Dx)2where Dx is the step size.
Solution:
Similar to Homework #3, Problem #2.
Express the above as a set of algebraic equations of the form
g(y)=0 and solve numerically with Newton's method.
Solution:
See Homework #13, Problem #3.
e=f2
Solution:
Similar to Homework #14, Problem #2.
e=1/f2
Ellipse: x2/a2 + y2/b2 = 1
Solution:
Solution:
drag ~ v2*sin2y*(cross-section area)Set up all the relevant equations, clearly mark which ones are to be solved, and briefly comment on how to proceed. Provide any missing information that you may need to solve the problem (e.g., are there anything else that needs to be specified?). Solve the equations only when time permits. Suggested by Anonymous
Solution:
The problem statement does not specify the final position xf.
We either specify it or impose a condition to determine it.
|