- Nam Sun Wang
- Department of Chemical Engineering
- University of Maryland
- College Park, MD 20742-2111
- 301-405-1910 (voice)
- 301-314-9126 (FAX)
e-mail: nsw@eng.umd.edu
©2002 by Nam Sun Wang
- Department of Chemical Engineering
Solution:
We apply Newton's method to find the root to f(x)=0
x=exp(a)
ln(x)=a
Let f(x) = 0 = ln(x)-a
Initial guess: x0=1 (or any positive value)
Iterate until convergence to the specified accuracy.
xk+1=xk-f(xk)/f'(xk)
=xk-(ln(xk)-a)*xk
Note that x0£0 is a
bad initial guess because the logarithmic function is undefined
for negative arguments.
óp
A(f) = ô f(x)*cos(j*x)*dx
õ0
Solution:
No, it is not a transformation.
Reason: Input to A is a function (i.e., vector); output from A
is a series of discrete numbers (a different type of
vector). A transform turns one vector into another vector
defined in the same vector space; it does not make a human out of
a monkey. Since A is not a transform, we cannot speak of a
linear transform. A is called a Fourier cosine series (except
for the first term j=0), not a Fourier cosine transform.
Solution:
Strictly speaking, no. Reason: Input to A is a function of x
(i.e., vector) defined in x=[0, p];
output from A is a function of j defined in a different interval j=real
number=[-¥, +¥] (i.e., another type of
vector). The answer would be "yes", had we restricted j to just [0, p]. A is linear.
Reason:
Addition: A(f+g)=A(f)+A(g)
Multiplication by a scalar: A(af)=aA(f)
A is called a Fourier cosine transform. (There are many
different definitions of Fourier transforms.) If we wish
to extend f(x) to x=[-¥, +¥], we implicitly assume
that f(x) repeats with a period of 2p
(because the cosine basis functions have a common period of
2p).
Furthermore, f(x)=f(2p-x)
(equivalently, f(p-x)=f(p+x) symmetric around x=p).
Solution:
As commented in the last part, if j=[-¥, +¥],
A is not a transformation; thus, we cannot even speak of
eigenvalues. However, had j been restricted to [0, p], we would find eigenvalues and eigenvectors by
solving the following equation.
óp
A(f) = ô f(x)*cos(j*x)*dx = l*f(j)
õ0
There is no function f(x) that satisfies the above relationship; thus,
A has no eigenvalue-eigenvector. If there were eigenvalues and
eigenvectors for Fouirer transform, Fourier transform would have
been much easier to evaluate.
f(x) = Sj=0n ajcos(j*x) f(x) = a0cos(0x) + a1cos(x) + a2cos(2x) + ... + ancos(n*x)Find the coefficients for the Fourier cosine series. You may find the following definite integrals useful.
For any integers j and k óp ô cos(j*x)*cos(j*x)*dx = p/2 for j=1,2... õ0 óp ô cos(j*x)*cos(k*x)*dx = 0 for j=0,1,2,... k=0,1,2,... j¹k õ0
Solution:
The projection formula for orthogonal basis vectors gives the
coefficients.
aj = (f,cos(j*x)) / (cos(j*x),cos(j*x))
1 óp
a0 = - ô f(x)*dx
p õ0
2 óp
aj = - ô f(x)*cos(j*x)*dx for j=1,2...
p õ0
The first term j=0 has a slightly different formula because
óp
ô 1*1*dx = p
õ0
Solution:
Yes, because we are not dealing with an interval
symmetrically centered at x=0. For example, had the interval been
x=[-p/2,+p/2] instead, we would not have been able to capture
f(x)=x. On the other hand, this cosine series cannot capture a constant,
say, f(x)=1.
Solution:
Normally, a series of sine functions complete the Fourier
series, if x covers any intervals of 2p.
sin(j*x) for j=1,2,...
With a combination of both sine and cosine basis functions, we
can approximate any continuous function f, given sufficient
number of terms. (Sidebar: Caution! In x=[0, p], sine and cosine functions are NOT
mutually orthogonal, because the periodicity of these functions
is 2*p, not p. Orthogonality depends on
the definition of the scalar product, which, in turn, depends on
the limits of integration. Incidentally, sine and cosine
function are orthogonal over any 2p intervals, not just [0, 2p] or [-p, p].
For any real a and for any integers j and k,
óa+2p
ô cos(j*x)*cos(k*x)*dx = 0 for j=0,1,2,... k=0,1,2,... j¹k
õa
óa+2p
ô sin(j*x)*sin(k*x)*dx = 0 for j=0,1,2,... k=0,1,2,... j¹k
õa
óa+2p
ô cos(j*x)*sin(k*x)*dx = 0 for j=0,1,2,... k=0,1,2,... j¹k
õa
The given set of cosine basis functions over x=[0, p] already constitutes a complete set, and
there is no need to add more. Had there been a need to include
more basis functions, the additional basis functions do not need
to be orthogonal; they only have to be linearly independent.
Orthogonality condition is not required, although it greatly
simplifies evaluating the coefficients.
Solution:
First, map the given interval x=[a, b] to
x'=[0,1].
x'=p*(x-a)/(b-a)
Or
x=(b-a)/p*x' + a
Then, we apply the Fourier cosine formula derived earlier.
1 óp
a0 = - ô f((b-a)/p*x'+a)*dx'
p õ0
2 óp
aj = - ô f((b-a)/p*x'+a)*cos(j*x')*dx' for j=1,2...
p õ0
Solution:
This problem is similar to the material presented in
"Eigenvalue-Eigenvector of Differential Operator & Fourier Series"
of the lecture notes posted on 11/21.
Here are some examples.
dx1/dt = 2*x1 - x2 + 9 I.C.: at t=0, x1(0)=1 dx2/dt = 9*x1 - 4*x2 + 38 x2(0)=2
Solution:
Find steady-state(s)
dx1/dt = 0 = 2*x1 - x2 + 9
dx2/dt = 0 = 9*x1 - 4*x2 + 38
--> x1ss=-2 x2ss=5
Introduce deviation variables around the steady-state
X1=x1-x1ss=x1+2
X2=x2-x2ss=x2-5
The ODEs become,
dx1/dt = 2*(x1+2) - (x2-5)
dx2/dt = 9*(x1+2) - 4*(x2-5)
The ODEs expressed in deviation variables in the "standard" linear/homogeneous form,
dX1/dt = 2*X1 - 1*X2 I.C.: at t=0, X1(0)=3
dX2/dt = 9*X1 - 4*X2 X2(0)=-3
In vector notation:
dX/dt = A*X
where A=[2 -1]
[9 -4]
Eigenvalue of A
A*X=l*X
det|A-lI|=0=-l2-2*l-1
--> l1=-1 l2=-1 (repeated)
Eigenvector of A
v1=[1]
[3] (There is only 1 linearly independent eigenvector.)
Jordan normal form.
(lI-A)*v2=v1
v2=[1]
[2]
Similarity transform A*V=V*J
[2 -1][1 1] = [1 1][-1 1]
[9 -4][3 2] [3 2][ 0 -1]
Solution:
X = exp(A*t)*X0 = V*exp(J*t)*V-1*X0
X = [1 1][exp(-t) t*exp(-t)][-2 1]*[ 3]
[3 2][ 0 exp(-t) ][ 3 -1] [-3]
X = [ 3*exp(-t) + 12*t*exp(-t)]
[-3*exp(-t) + 36*t*exp(-t)]
Finally,
x1 = X1+x1ss = X1-2 = 3*exp(-t)+12*t*exp(-t)-2
x2 = X2+x2ss = X2+5 = -3*exp(-t)+36*t*exp(-t)+5
A*v=l*vOr, equivalently,
(A-lI)*v=0Pre-multiplying the inverse of (A-lI) to both sides of the last equation yields,
v=(A-lI)-1*0Can we apply the last equation to find an eigenvector? Briefly, why?
Solution:
No, because eigenvalues satisfy det|A-lI|=0.
In other words, (A-lI) is singular;
(A-lI)-1 does not exist.
(A-lI)*v2=v1Pre-multiplying the inverse of (A-lI) to both sides of the last equation yields,
v2=(A-lI)-1*v1Can we generate v2 from v1 based on the last equation? Briefly, why?
Solution:
No, because eigenvalues satisfy det|A-lI|=0.
In other words, (A-lI) is singular;
(A-lI)-1 does not exist.
dx
-- = (m-D)*x
dt
ds m*x
-- = D*(sf-s) - ---
dt Y
where m(s)=specific growth rate, which is a function of s.
Y=yield coefficient=constant
sf=substrate concentration in the feed stream=constant
D=F/V=dilution rate (inverse of residence time)=constant
--------->+ +------->
F, sf | | F, x, s
| | | |
| |
| V, x, s |
| |
+---------+
bioreactor
Solution:
Set d/dt=0:
0 = dx/dt = (m-D)*x --> xss=0 or m(sss)=D
0 = ds/dt = D*(sf-s) - m/Y*x --> sss=sf or xss=Y*(sf-sss)
The number of steady-states is the number of times the function
m(sss) crosses D plus 1.
Solution:
Introduce deviation variables X and S.
X=x-xss
S=s-sss
Linearize around steady-state points (xss,sss)
dX/dt = (m(sss)-D)*X + m'(sss)*xss*S
dS/dt = -m(sss)/Y*X - [D+m'(sss)/Y*xss]*S
Steady-state I (washout): xss=0, sss=sf
dX/dt = (m(sf)-D)*X
dS/dt = -m(sf)/Y*X - D*S
Characteristic equation: [l-(m(sf)-D)]*[l-(-D)] = 0
Eigenvalues are l1=m(sf)-D
and l2=-D
Stable if l1<0 & l2<0 --> m(sf)<D
Steady-state(s) II (nonwashout): xss=Y*(sf-sss)¹0, m(sss)=D
dX/dt = m'(sss)*xss*S
dS/dt = -D/Y*X - [D+m'(sss)/Y*xss]*S
Characteristic equation:
0=det[l -m'(sss)*xss ]
[D/Y l+D+m'(sss)/Y*xss]
l[l+D+m'(sss)/Y*xss] + D*m'(sss)*xss/Y = 0
l2 + (D+m'(sss)/Y*xss]*l + D*m'(sss)*xss/Y = 0
The characteristic equation is a quadratic equation of the form: l2 + b*l + c = 0
l<0 if -b+sqrt(b2-4*a*c)<0
--> 4*a*c>0 --> c=D*m'(sss)*xss/Y>0
Stable if m'(sss)>0 ... slope condition for stability
D=1, m(s)=s, sf=1, Y=0.5, x(0)=1, s(0)=1
Solution:
Euler's method for numerical integration of ODEs.
With the given parameter values, the dynamic equations are:
dx/dt = (s-1)*x
ds/dt = (1-s) - 2*s*x
Start with x=1 and s=1, and the slopes at this starting point are:
dx/dt = (1-1)*1 = 0
ds/dt = (1-1) - 2*1*1 = -2
With h=0.1, advance to t=0.1
x(0.1) = x(0) + dx/dt*h = 1 + 0*0.1 = 1.0
s(0.1) = s(0) + ds/dt*h = 1 + -2*0.1 = 0.8
y=a1*X<1> + a2*X<2> +... + an*X<n> y=X*aWe choose the regression coefficient a by minimizing the sum of squared errors (sse).
Min sse=eTesubject to equality constraint
e=y-X*a where y and X are given
Solution:
Rewrite the problem in the "standard" form. Specifically,
express the equality constraint in the "f=0" format.
Min L=sse=eTe
subject to equality constraint: f=0=e-y+X*a
In the above two equations, y and X are given; the remaining
variables a and e are to be solved for. Thus, the variables that
corresponds to the generic unknown vector variable "u" in our
lecture note are "a" and "e" in this problem. We couple the
equality constraint to the quantity to be minimized.
Min H(e,a,l)=L+lTf=eTe+lT(e-y+X*a)
The necessary condition for a minimum is dH=0:
Eqn 1. d H(e,a)/de = 0 = 2*eT + lT
Eqn 2. d H(e,a)/da = 0 = lT*X
Eqn 3. d H/dl=0=f (This recovers the equality constraint "f=0")
Solve the above three equations for the three variables: a, e, and
l.
Solution:
0 = 2*eT + lT --> lT=-2eT
0 = lT*X --> eT*X=0 (The error vector is orthogonal to the basis vectors.)
Finally, substituting into the equality constraint, we get,
--> 0=XT*e=XT*(y-X*a)
--> a=(XT*X)-1*XT*y
Thus, we recover the normal equation in linear regression.
0£aDerive the necessary condition(s) for a minimum. (An example of this inequality constraint: the counterfeit example from the last exam where we try to approximate the green color on a $100 bill with 5 different types of ink -- the regression coefficients must be all non-negative because we cannot physically mix a negative quantity of ink. Concentrations must not be negative.)
Solution:
In lecture, we derived the necessary condition for inequality
constraints with the notation f(u)£0.
Min L=sse=eTe
subject to equality constraint: f=0=e-y+X*a
plus inequality constraint: g=-a£0
Couple each constrain to H with a Lagrangian multiplier.
Min H(e,a)=eTe+lT(e-y+X*a)+mT(-a)
Necessary condition:
d H(e,a)/de = 0 = 2*eT + lT
d H(e,a)/da = 0 = lT*X - mT
d H/dlT=0=f (just the equality constraint)
d H/dmT=0=g (just the active inequality constraint)
At the boundary, m³0, g=-a=0
Inside he boundary, m=0, g=-a<0
ótf
Min J=f(x(tf),tf) + ô L(x,u,t)*dt
õt0
subject to
dx/dt=f(x,u,t) x(t0)=x0=given
y(x(tf),tf)=0 ... equality constraint at tf
where the subscripts 0 and f denote the initial value at
t0 and the final value at tf respectively
of the corresponding variable.
Solution:
We fold the equality constraint y=0
into the original J with a Lagrangian multiplier.
J'new = Joriginal + mT*y
In addition, we allow tf to vary, since it is
unspecified.
dJ' = ...(same old stuff from the lecture notes)... + dJ'/dtf*dtf
We then set to zero the factor preceeding dtf.
dJ'/dtf=0
The rest of the necessary conditions remain identical to the case
where the final time tf is specified and where the
equality condition y=0 is absent.
dx
-- = (m-D)*x
dt
ds m*x
-- = D*(sf-s) - ---
dt Y
where m(s)=specific growth rate, which is a function of s.
Y=yield coefficient=constant
sf=substrate concentration in the feed stream=constant
D=F(t)/V(t)=dilution rate, which is a function of t.
--------->+
F, sf |
| | |
| |
| V, x, s |
| |
+---------+
bioreactor
Set up all the relevant equations that are needed to find the
feeding profile F(t) and harvest time t1 such that
we maximize biomass production and minimize the time
required. One possible objective function is,
Max J= V1*x1 - w*t1
where w=weighting factor
(ratio of operation cost to the price of biomass product)
The subscript 1 indicates the various values at harvest time.
Additional information: In a sustained cyclic operation, the
initial concentration is identical to the final concentration
during harvest: x0=x1,
s0=s1. For the organism to grow well, we
need a 5% inoculum, i.e., V0=0.05*V1.
Clearly identify the state variables and mark which equations are
to be solved along with relevant boundary conditions. Provide
any missing information that you may need to solve the problem
(e.g., are there anything else that needs to be specified to
close the problem?). Do not actually solve the equations because
of limited time. However, do describe briefly how to proceed
after equations have been set up.
Solution:
This problem tests whether we can apply the necessary
condition for optimality to a specific example.
d ds --[r2*(--)] = r2*e*v(s) dr dr B.C. at r=1, s(1)=1; at r=0, ds(0)/dr=0where e is the ratio of the characteristic reaction rate to the characteristic diffusion rate. Suggested by HN
Solution:
This is a regular perturbation problem (because the small
parameter is associated with the highest derivative term). It is
similar to the "Regular Perturbation Method: Diffusion-Reaction"
example of the lecture notes posted on 11/25, except the equation
is in spherical coordinate rather than cartesian coordinate.
In O(1) approximation, we ignore e.
d ds
--[r2*(--)] = 0 --> s0(r)=1
dr dr
To find the existence/thickness/location of the boundary layer,
we resealing the given ODE with either R=r/d or R=(1-r)/d.
d ds
--[R2*(--)] = R2*d2*e*v(s)
dR dR
Consider the following cases:
LHS dominates ... d=1
--> We recover the outer solution.
RHS dominates ... d>>1
--> We extent r way beyond r=1. --> No.
LHS and RHS balance ... d=sqrt(1/e)
--> We extent r way beyond r=1. --> No.
We cannot rescale other than d=1.
Thus, there is no boundary layer.
Solution:
This is a singular perturbation problem. We solve the outer
region and inner region separately. In terms of e'=1/e<<1,
we have,
d ds
e'--[r2*(--)] = r2*v(s)
dr dr
Outer solution: In O(1) approximation, we ignore e' (the LHS).
s(r)=s0(r)+e's1(r)+...
v(s0)=0 --> s0(r)=constant=0
where s0 is a root of the algebraic equation v(s)=0.
For a nonlinear v(s), it is in general possible to have multiple
roots. In practice, s=0 is normally a root (because there is no
reaction without reactant). We now rescale the independent
variable r to rewrite the ODE valid in the inner region (i.e.,
within the boundary layer).
Re-scale: R=(1-r)/d dR=-dr
d dS
e'(---)[(1-dR)2*(---)] = (1-dR)2*v(S)
ddR ddR
Matching the LHS with the RHS, we obtain the thickness of the boundary layer.
e'/d2=1
--> d=e'1/2
d dS
--[(1-e'1/2R)2*(--)] = (1-e'1/2R)2*v(S)
dR dR
We should not linearize v(S) around the value of S at r=1 or R=0
because S is not a constant but is rapidly changing in the
boundary layer region. We expand S in power series of e'1/2 to capture the e'1/2 terms in the rescaled ODE.
Let S(R)=S0(R)+e'1/2*S1(R)+e'*S2(R)+...
d2S0
--- = v(S0) B.C. S(0)=1 limR->¥S0(R)=limr->1s0(r)=0
dR2
Inner Solution: The above is a nonlinear ODE. Multiplying by dS0/dR yields,
(dS0/dR)(d2S0/dR2)=d(dS0/dR)2/dR/2=(dS0/dR)*v(S0)
óS0
(dS0/dR)2 = ô 2*v(S0)*dS0
õ0
Let us assume certain forms of v(s) so that we can integrate the RHS, although
it is not necessary that we express the RHS in a closed form.
For v(s)=1 (0th order reaction kinetics)
(dS0/dR)2 = 1
dS0 = -dR
--> S0 = 1-R (for R<1)
S0 = 0 (for R³1) (Concentration cannot be negative.)
This matches with v(s0)=0 in the outer solution.
For v(s)=s (1st order reaction kinetics)
(dS0/dR)2 = S02
1/S0*dS0 = -dR
--> S0 = exp(-R)
This matches with v(s0)=0 in the outer solution.
For v(s)=s/(b+s) (saturation kinetics)
:
:
Uniformly valid solution:
soverall = sout + Sin - smatching
= 0 + exp(-R) - 0
= exp(-(1-r)/e'1/2) = exp(-e1/2(1-r))
Solution:
We can "flatten" out the spherical geometry by introducing
z=s*r. Doing so changes the diffusion equation from a spherical
coordinate to a cartesian coordinate.
d2z
--- = e*r*v(z/r) = e*z
dr2
B.C.: at r=0, ds(0)/dr=0, z=r*s=0
at r=1, s(1)=1, z(1)=r*s=1
Let z=exp(lr)
--> characteristic equation for eigenvalue: l2-e=0
--> l1=e1/2, l2=-e1/2
z(r)=s(r)*r = c1*exp(l1r) + c2*exp(l2r)
Apply B.C. to find the integration constants c1 and c2.
At r=0, z(0)=0=c1+c2
At r=1, z(1)=1=c1*exp(l1) + c2*exp(l2)
--> c1=-c2=1/(exp(l1)-exp(l2))
Finally,
z(r) = (exp(l1r)-exp(l2r))/(exp(l1)-exp(l2))
s(r) = z(r)/r = (exp(l1r)-exp(l2r))/r/(exp(l1)-exp(l2))
|