- 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
©2001 by Nam Sun Wang
- Department of Chemical Engineering
dr/dt = a*r - b*r*f ... r=rabbits (prey) df/dt =-c*f + d*r*f ... f=foxes (predator) where a, b, c, and d are model parameters
Parameters: a = 2 b = 2 c = 1 d = 1 Initial conditions: r(0)=2 f(0)=1
Solution:
Euler's method for numerical integration of ODEs.
With a=2, b=2, c=1, and d=1, the dynamic equations are:
dr/dt = 2*r - 2*r*f = 2*r*(1-f)
df/dt = -f + r*f = f*(r-1)
Start with r=2 and f=1, and the slopes at this starting point are:
dr/dt = 2*2*(1-1) = 0
df/dt = 1*(2-1) = 1
With h=0.1
r (at t=0.1) = r(0) + dr/dt*h = 2 + 0*0.1 = 2.0
f (at t=0.1) = f(0) + df/dt*h = 1 + 1*0.1 = 1.1
Solution:
Runge-Kutta's method for numerical integration of ODEs.
1st evaluation is the same as Euler's Method.
r1 = 2
f1 = 1
1st slope for r and f (from Problem 1) ...
kr1= fr(r1,f1) = 0
kf1= ff(r1,f1) = 1
2nd evaluation is at:
r2 = 2 + kr1*h/2 = 2 + 0*0.05 = 2.0
f2 = 1 + kf1*h/2 = 1 + 1*0.05 = 1.05
2nd slope for r and f ...
kr2= fr(r2,f2) = 2*2.0*(1-1.05) = -0.2
kf2= ff(r2,f2) = 1.05*(2.0-1) = 1.05
3rd evaluation is at:
r3 = 2 + kr2*h/2 = 2 + -0.2 *0.05 = 1.99
f3 = 1 + kf2*h/2 = 1 + 1.05*0.05 = 1.0525
3rd slope for r and f ...
kr3= fr(r3,f3) = 2*1.99*(1-1.0525) = -0.20895
kf3= ff(r3,f3) = 1.0525*(1.99-1) = 1.041975
4th evaluation is at:
r4 = 2 + kr3*h = 2 + -0.20895 *0.1 = 1.979105
f4 = 1 + kf3*h = 1 + 1.041975*0.1 = 1.1041975
u = 9.8116/(1+9.8116) = 0.9075
4th slope for r and f ...
kr4= fr(r4,f4) = 2*1.979105*(1-1.1041975) = -0.41244
kf4= ff(r4,f4) = 1.1041975*(1.979105-1) = 1.08113
Weighted average slope for r and f ...
krave = (kr1+2*kr2+2*kr3+kr4)/6 = ( 0-2*0.2 -2*0.20895-0.41244)/6 = -0.20506
kfave = (kf1+2*kf2+2*kf3+kf4)/6 = ( 1+2*1.05+2*1.04198+1.08113)/6 = 1.04418
Update the values of r and f at t=h=0.1
r = r + krave*h = 2 + -0.20506*0.1 = 1.97949
f = f + kfave*h = 1 + 1.04418*0.1 = 1.10442
dx [a*x1 - b*x1*x2]
-- = B(x) =
dt [-c*x2 + d*x1*x2]
where x is a column of two functions
x = [x1]
[x2]
Is transformation B defined above linear? Why?
Solution:
No, because B fails both linearity tests:
B(x+y)¹B(x)+B(y).
B(a*x)¹B(x) where a is a scalar.
Specifically, the product term x1*x2 makes the transformation nonlinear.
Solution:
Do not even bother with eigenvalues and eigenvectors of B.
The eigenvalue-eigenvector idea applies to linear
transformations, and B is not linear. (Thus, in the following
parts, we linearize the system of ODEs around steady-state
points first before we apply eigenvalue-eigenvector to find
stability. If B were already linear to start with, we do not
need to examine the linearized approximation of it.)
Solution:
Set d/dt=0:
0 = a*r - b*r*f = r*(a-b*f)
0 = -c*f + d*r*f = f*(-c+d*r)
There are two solutions:
Steady-state I: rss=0, fss=0
Steady-state II: rss=c/d, fss=a/b
Solution:
Introduce deviation variables R and F.
R=r-rss
F=f-fss
Linearize around steady-state points (rss,fss)
dR/dt = (a-b*fss)*R - b*rss*F
dF/dt = d*fss*R + (-c+d*rss)*F
Steady-state I: rss=0, fss=0
dR/dt = a*R
dF/dt = -c*F
Characteristic equation: (l-a)*(l+c)=0
Eigenvalues are l=a and l=-c; saddle point; unstable.
(Note that the real part of all eigenvalues, not just some, must be negative for a steady-state to be stable.)
Steady-state II: rss=c/d, fss=a/b
dR/dt = - b*c/d*F
dF/dt = d*a/b*R
Characteristic equation: l2+a*c=0
Eigenvalues are l=+/-i*sqrt(a*c); purely imaginary; center.
Thus, the trajectories circle around in a periodical fashion.
Solution:
The trajectories are given by a series of functions
g(r,f,r0,f0)=0
Or, equivalently,
h = h(r,f) = h(r(t),f(t)) = h(r(0),f(0)) = h(r0,f0) = constant
which means we seek
d(h(r,f))/dt=0.
by suitably combining the given dr/dt and df/dt ODEs.
multiply by 2*a (r*(dr/dt)=(dr2/dt)/2=(a-b*f)*r2)
multiply by 2*b (f*(df/dt)=(df2/dt)/2=(-c+d*r)*f2)
------------------------------------------------------
Add: a*(dr2/dt) + b*(df2/dt) + c*(d(r*f)/dt) + d*(dr/dt) + e*(df/dt) = 0
d/dt(a*r2+b*f2+c*r*f+d*r+e*f) = 0
a*r2+b*f2+c*r*f+d*r+e*f = constant
The above equation gives a series of trajectories for different initial conditions.
Here is another approach that eliminates the "t" component (trajectories residing in the phase diagram are functions of (r, f) but not explicit function in t)
dr/df = (dr/dt)/(df/dt) = (r*(a-b*f))/(f*(-c+d*r))
(-c+d*r)/r*dr = (a-b*f)/f*df
Integrate: -c*ln(r)+d*r = a*ln(f)-b*f + constant
Rearrange: ln(fa*rc) - d*r - b*f = constant = ln(f0a*r0c) - d*r0 - b*f0
where r0 and f0 are initial conditions.
Solution:
^
r | + /----\ -
| / \ Direction of rotation: clockwise
| | | (because fox follows rabbit)
c/d +- | * |
| \ /
| \------/
| - +
+-------+-------->
0 a/b f
slope of trajectories: dr/df = (dr/dt)/(df/dt) = (r*(a-b*f))/(f*(-c+d*r))
Along the vertical line that passes through the steady state point (i.e., f=a/b), slope=0
Along the horizontal line that passes through the steady state point (i.e., r=c/d), slope=±¥
Consider only the quadrant where 0£r & 0£f (becaue negative population densities make no physical sense)
Within this quadrant, the slope of the trajectories are identified in the phase diagram with "+" or "-".
time rabbits foxes ---------------------- : : : : : : : : :Set up the equation(s) needed to estimate the parameters a, b, c, and d in the model from the above population data. Note that field data have measurement errors.
Solution:
Because of measurement errors, avoid differentiation; regress based
on integrated measurements.
Given dynamic equations are:
dr/dt = a*r - b*r*f
df/dt = -c*f + d*r*f
Integrate (or approximation with summation) the above dynamic equations:
r-r0 = a*Ri - b*RFi
f-f0 = -c*Fi + d*RFi
For continuous measurement, the symbols Ri, Fi, and RFi are:
Ri = integral of r(t) from t=0 to t=t.
Fi = integral of f(t) from t=0 to t=t.
RFi = integral of r(t)*f(t) from t=0 to t=t.
For discrete measurement, the symbols Ri, Fi, and RFi are:
Ri = summation of ri from i=0 to i.
Fi = summation of fi from i=0 to i.
RFi = summation of ri*fi from i=0 to i.
Regress r & f against Ri, Fi, RFi to find the model parameters a, b, c, and d
[a]
[r-r0] = [ Ri -RFi 0 0 ][b]
[f-f0] [ 0 0 -Fi RFi ][c]
[d]
The solution is given by the normal equation for linear regression
a = (XT*X)-1*XT*y
where a is a vector of the model coefficients, X is a matrix of the
integrated state variables, and y is a vector of the measured state
variables.
[a]
y=[r-r0] X=[ Ri -RFi 0 0 ] a = [b]
[f-f0] [ 0 0 -Fi RFi ] [c]
[d]
A(x)=[cos(f)*x1-sin(f)*x2]
[sin(f)*x1+cos(f)*x2]
where f is a scalar, and x is a column of two numbers
x=[x1]
[x2]
Solution:
An orthogonal transformation is one that does not change the length.
See lecture note on linear transformation.
Let y=A*x
Orthogonal transformation: (y,y)=(A*x,A*x)=(x,x) ... no change in length
Thus, we need to show that the above relationship is true for the linear
transformation defined in question.
For a vector x=(x1,x2)T,
First, define scalar product: (x,y)=x1*y1+x2*y2
|A*x|2 = (cos(f)*x1-sin(f)*x2)2 + (sin(f)*x1+cos(f)*x2)2
= (x1)2 + (x2)2 = |x|2
Do not confuse the definition of an orthogonal transformation,
which is |A*x|=|x|, with its consequences, one of which is
(Ax,Ay)=(x,y).
Do not confuse an orthogonal linear transformation with
orthogonal vectors. (Note the plural; we need to have two
vectors to even speak of the angle between them).
L(y)=(d2/dx2+a)*y = d2y/dx2+a*y B.C. y(0)=y(1)=0 where a is a scalar.
Solution:
Eigenvalue-eigenvector equation:
L(y)=d2y/dx2+a*y=l"y B.C. y(0)=y(1)=0
where
l" is the eigenvalue(s) for the linear operator L.
With minor rearrangement, we have,
d2y/dx2=(l"-a)y
d2y/dx2=l'y
Thus,
l'=(l"-a) is the eigenvalue of the second-derivative operator d2/dx2
l'=(l"-a) is the eigenvalue of the second-derivative operator d2/dx2
Eigenvalues of d2/dx2: -l2=l'
lj=p*j for j=0,1,2,..
The eigenvalues for this problem are simply shifted from those for the the plain second-derivative operator.
Eigenvectors of d2/dx2:
yj=cos(lj*x) and yj=sin(lj*x)
Of the above, keep the sine series to satisfy the boundary conditions.
Since sin(0*x)=0, we lose this function. The eigenvectors are:
sin(p*x), sin(2p*x), sin(3p*x),...
Solution:
infinite
Solution:
Yes, because symmetrical linear transforms have orthogonal eigenvectors.
ó1
Define scalar product: (f,g)= ô f(x)*g(x)*dx
õ0
integral (from 0 to 1) f(x)*g(x)*dx (if you don't see the integral symbol above)
(Note that orthogonality depends on how we define a scalar product.
Had we defined it differently, the eigenvectors would not be orthogonal.)
Solution:
Yes. Given n linearly independent
basis vectors (functions), we
can exactly represent any given vector (function) f taken from an
n-dimensional LVS as
a linear combination of other linearly independent basis vectors
(functions) yj taken from the same LVS.
Eigenvectors are linearly independent; thus, they conveniently
serve as the basis vectors.
f=a1*y1+a2*y2+a3*y3+...
Since the eigenvectors in this problem are orthogonal, we apply
the projection formula to find the coefficients aj
aj=(f,yj)/(yj,yj) j=1,2,...
where yj=sin(lj)=sin(jpx)
Had the eigenvectors being non-orthogonal, we would not be able to apply the
above formula.
The dimension of a LVS composed of continuous functions is
infinite; thus, we need an infinite number of basis vectors to
represent a given function exactly. Of course, infinity is a
mathematical abstraction; in practice, we only approximate a
given vector (function) with a finite number of terms. Thus, we
expect a linear combination of the basis vectors to match the
given vector (function) only at finite number of points. Not
being able to match the two end points should not bother us.
Lastly, it makes no sense to speak of the boundary
conditions of a given function f. Here is a question to tickle
your brain: what if the interval was x=[-1 1] instead of x=[0 1]?
Solution:
Simply map x=[a,b] into z=[0,1] with the function
z=(x-a)/(b-a).
d2y/dx2=-l*y B.C. y(0)=y(1)=0 where l is a scalar.
Solution:
(This was the way the question was originally proposed to me;
it was not a valid question.)
Strictly speaking, no. Eigenvalue and eigenvectors applies to
linear transformation. We can speak of a second-derivative
operator D2 = d2/dx2; we can
speak of an operator that multiplies a scalar constant l* (which is not to be confused with the
eigenvalue); we can speak of another operator
which is a linear combination of the above two operators
d2/dx2+l*. We
can speak of an eigenvalue-eigenvector equation given by the
above ODE (which really means eigenvalue-eigenvector of a
second-derivative operator, but not eigenvalue-eigenvector the given ODE itself.
The ODE itself, the LHS, or the RHS, is a vector (i.e., function),
not a linear transformation. We cannot even speak of the eigenvector-eigenvalue
of a vector, only a linear transformation.
|