BOT/ECOL 5756 Ecological Systems Modeling
Lecture 8
Lecture 8: Introduction to differential equations: Analytical solutions
Motivation • • • An important part of the analysis phase is finding the solution(s) to your model. Whenever possible, you should try to calculate exact (analytical) solutions; i.e., solve the differential (or difference) equations that describe the system "by hand". For large or complex models, finding analytical solutions may be impossible. In this case, you will need to employ numerical methods to solve your model (e.g., use STELLA; we will discuss numerical methods next time). Additionally, when analytical solutions are not possible, you should at least solve for the steadystate behavior (this is good practice regardless of whether or not closed-form solutions exist). For large models that can be decomposed into sets of smaller (sub-)models, you should try to find solutions to the sub-models.
• •
Steady-state vs. equilibrium Let’s begin with analytical methods for finding the steady-state solutions. But, first we should clarify what it means to be in a “steady-state” and how this differs from being “at equilibrium.” • • • • The terms equilibrium and steady-state are often used interchangeably in ecology, although many disciplines (physyology, chemistry, physics, engineering) have specific meanings for each. In general, equilibrium describes the average condition of a system (e.g., a state variable over some specific period of time) and, steady-state connotes no change. Examples of usage include: o Equilibrium (or “steady-state equilibrium”) is an average condition of a system where the trajectory remains unchanged in time, e.g.:
(1)
BOT/ECOL 5756 Ecological Systems Modeling
Ogle, Fall 2006
o Thermodynamic equilibrium describes a condition in a system where the distribution of mass and energy moves towards maximum entropy, e.g.:
o Dynamic equilibrium occurs when there are unrepeated average states through time: e.g.:
Examples also include “limit cycles” or “periodic cycles”– recall the “classical” lynx-hare problem. Here’s a plot of the predicted hare pop’n density from King & Schaffer (2001) Ecology 82:814-830
(2)
BOT/ECOL 5756 Ecological Systems Modeling
Ogle, Fall 2006
o Steady-state (or “static equilibrium”) occurs where force and reaction (or input and output) are balanced and the properties of the system remain unchanged over time, e.g.:
o Transient state. Note that there is almost always a transient period before the system reaches an equilibrium or steady-state whereby the system properties are dynamic and changing over time, e.g. consider the 3-compartment (closed-system) model from Lecture 7:
120
Transient state(s)
100
x1(t) x2(t) x3(t)
x1(t), x2(t), x3(t)
80 60 40 20 0 0 20 40 60 80 100
Steady-state(s)
Time (t)
Computing steady-state solutions • • Using the definition above: steady-state (or “static equilibrium”) occurs when the inputs and outputs are balanced and the properties of the system remain unchanged over time. That is, an individual compartment i is in a steady-state if: d x i (t ) = 0 or x i (t ) = 0 ɺ dt
(7.1)
(3)
BOT/ECOL 5756 Ecological Systems Modeling
Ogle, Fall 2006
And the ENTIRE system is in a steady-state if all of the time derivatives are zero for ALL compartments; that is, if Eqn (7.1) is true for all i = 1, 2, .., n states (or compartments). • • Thus, we find the steady-state solutions, call these x i* , by solving Eqn (7.1) for xi and the solution gives us x i* . Some examples: (1) Consider the model for soil water content (W): dW = P − λ ⋅ B ⋅W − ε ⋅W (7.2) dt Where P is precipitation input, B is root mass density, λ is a root water extraction parameter, and ε is an evaporation rate constant. Find the steady-state water content (W*), assuming that P, B, λ, and ε do not change with time. The solution is arrived at by solving the following equation for W*: P − λ ⋅ B ⋅ W * −ε ⋅ W * = 0 (7.3)
W* = P λ ⋅B +ε (7.4)
(2) Consider the logistic growth model (see Lecture 4, pg 1) for population size (N):
d N (t ) N (t ) = r ⋅ 1 − ⋅ N (t ) dt K What is the steady-state solution N (i.e., find N*)?
(7.5)
(3) Recall the Mediterranean chaparral scrubland example from Lab 2 and the model for biomass dynamics of live biomass (LB) and dead biomass (DB): d LB( t ) = p ⋅ LB − s ⋅ LB dt d DB( t ) = s ⋅ LB − d ⋅ DB dt What are the steady-state solutions for LB and DB (call these LB* and DB*)?
(7.6)
(4)
BOT/ECOL 5756 Ecological Systems Modeling
Ogle, Fall 2006
(4) How about the nitrogen cycle model for the humid savanna (Lecture 7):
o What are the steady-state solutions for G*, H*, D*, L*, O*, M*, and N*? Since all of the rate coefficients (parameters) are constant (i.e., do not vary with time and do not depend on the state of system), we could solve for the steady-states by hand. Or, we could also use Maple to simultaneously solve a set of 7 equations. But, this becomes a little tricky. o For example, consider first solving for steady-state nitrogen in grass biomass (G*), this means solving the following equation for G*: 0 = u ⋅ N * ⋅G * + a ⋅ G * −bs ⋅ G * −d s ⋅ G * −d r ⋅ G * −c ⋅ G *
0 = G * ⋅ ( u ⋅ N * + a − bs − d s − d r − c )
(7.7)
o So, what are the possible steady-state values for G*? o Can try to solve for the steady-states in Maple... o The Maple solutions look messy, but they could be simplified, which would probably have to be done “by hand.” Mazancourt et al. (1999) give the simplified steady-state solution for nitrogen in grass biomass (G*) as a function of steady-state inorganic nitrogen (N*):
where:
(5)
BOT/ECOL 5756 Ecological Systems Modeling
Ogle, Fall 2006
o Graphical representation of grass biomass steady-state solutions:
Percentage of consumption (x)
Steady-state nitrogen in grass biomass as a function of the percentage of primary production consumed by the herbivores, at five fractions of nitrogen losses by herbivores: βh = 0.01, 0.10, 0.25, and 0.50; solid lines are short-term trends (~decades) and dashed lines are longterm trends (~centuries).
(6)
BOT/ECOL 5756 Ecological Systems Modeling
Lecture 8
Stability of steady-states
• • Once you have computed the steady-states for the system variables, you may be interested in evaluating the “stability” of these steady-state conditions. There is detailed literature and extensive theory on how to evaluate the stability of steady-states (and equilibrium points in general), and these details are described in many mathematical or theoretical biology texts. We will not go into detail about stability theory and analysis. The general ideas are: o A steady-state condition (or equilibrium point) can be unstable or stable o Stability is a property of the model o Consider the situation where the system is at a steady-state (as predicted by the model), if we perturb the system slightly (e.g., slightly increase or decrease the level of one or more of the state variables), then the steady-state condition is: unstable if the system does not return to the steady-state condition following the perturbation stable if it returns to the steady-state condition following the perturbation
•
Solving differential equations
• • • We are often interested in the transient states (e.g., responses to disturbance, environmental change, etc.) Thus, we need to know how the system and the state variables change over time I.e., if we have a system of equations that describe the time dynamics of our system such that: x 1 ( t ) = f 1( x , t , θ ) ɺ x 2 (t ) = f 2 ( x , t ,θ ) ɺ
⋮
x n (t ) = f n ( x , t ,θ ) ɺ where xi represents compartment i for i = 1, 2, ..., n different states. • • • •
(7.8)
The goal is to solve the system of equations in (7.8) to obtain solutions for x 1( t ), x 2 ( t ),..., x n ( t ) We’ve already discussed that there are two approaches to solving a system of differential equations: Analytical vs. Numerical approaches. But, we haven’t discussed how you actually compute analytical solutions. Before delving into analytical solution methods, we first need to understand the different types of differential equation models that we might work with.
(7)
BOT/ECOL 5756 Ecological Systems Modeling
Lecture 8
Difference vs. differential equations
•
Difference equations: time step is discrete, recursive relationship, Markov chain. Examples:
x n +1 = ρ ⋅ x n
x x n +1 = ρ ⋅ x n ⋅ 1 − n k •
1st order linear difference equation 1st order nonlinear difference equation
(7.9) (7.10)
Differential equations: time step is continuous. Examples:
d x (t ) = ρ ⋅ x (t ) dt d x (t ) x(t ) = ρ ⋅ x(t ) ⋅ 1 − dt k 1st order linear differential equation 1st order nonlinear differential equation (7.11) (7.12)
Ordinary vs. partial differential equations
• • •
ODE's: Ordinary differential equation. The unknown function, e.g., x(t) in Eqns (7.11) and (7.12), depends on a SINGLE independent variable (e.g., t ). Both of the above examples (exponential growth, Eqn (7.11); logistic growth, Eqn (7.12)) are ODE's (i.e., ρ and k are constants; x(t) only depends on the variable, t ). PDE's: Partial differential equations. The unknown function depends on TWO or MORE independent variables (e.g., time and “space”). For example: ∂ ∂2 u( t , z ) = k ⋅ 2 u( t , z ) + Q ∂t ∂z ∂ ∂2 w ( t , z ) = Dw ⋅ 2 w ( t , z ) ∂t ∂z Heat equation; u = u(t,z); Q = source/sink Simple diffusion equation; e.g., water vapor gradient along stomatal cavity, w = w(t,z) (7.13)
(7.14)
Systems of differential equations: Some terminology
•
ORDER: order of highest derivative Examples given in Eqns (7.11) and (7.12) are first order Examples given in Eqns (7.13) and (7.14) are second order LINEAR vs. NONLINEAR: if x is the unknown function, then the differential equation is linear if it is a linear function of x, x', x'', ... Examples given in Eqns (7.11), (7.13) and (7.14) are linear differential equations Example given in Eqn (7.12) is a nonlinear differential equation (notice that we have an x2).
•
(8)
BOT/ECOL 5756 Ecological Systems Modeling
Ogle, Fall 2006
d2 y( t ) = − g ⋅ sin ( y( t )) dt 2 Why is this a nonlinear differential equation? Another nonlinear example:
•
HOMOGENOUS vs. NON-HOMOGENOUS (the mathematical definition): for the ODE written in the form: d2 y dy + p( t ) ⋅ + q( t ) ⋅ y = g ( t ) 2 dt dt (First, note that this is a 2nd order, linear ODE.) Eqn (7.15) is HOMOGENOUS if g(t) = 0 Eqn (7.15) is NON-HOMOGENOUS if g(t) ≠ 0. E.g., g(t) = 4, g(t) = t2, etc. How would you classify the exponential growth model in Eqn (7.11), the logistic growth model in Eqn (7.12), and the heat equation in Eqn (7.13)? (7.15)
•
How would you classify the equations in the nitrogen cycle model in de Mazancourt et al. (1999)? ODE vs. PDE? Order? Linear vs. nonlinear? Homogenous vs. non-homogenous?
(9)
BOT/ECOL 5756 Ecological Systems Modeling
Ogle, Fall 2006
Analytical solutions
• •
The goal is to solve for the unknown function such as x(t), or the state variables at time t, to get an exact solution. Example, Eqn (7.11) for exponential growth; 1st order, linear, homogenous ODE, whereby dx = ρ ⋅ x . How do we go about solving this ODE for x(t)? dt
o
Need initial conditions: x(0) = x0. More generally, we write the initial conditions as x(t0) = x0; i.e., we don't have to start at t = 0, we may start at some other time, t = t0. We also note that Eqn (7.11) is separable. What is a separable differential equation?
o
Definition. Separable equations. For a 1st order differential equation:
Rewrite in form: M ( t , x ) + N ( t , x )
dx = f(t, x) dt
dx =0 dt If M is a function of only t (or is constant) and N is a function of only x (or is constant), then equation can be written as: M(t)dt = – N(x)dx. Each side of the equation depends on only one of the variables, and we solve for x by integrating the left-hand-side with respect to t and the right-hand-side with respect to x.
o
Now note that Eqn (7.11) is separable and thus can be written as: 1 dx = ρ ⋅ dt x (7.16)
o
Now we integrate the left-hand-side with respect to (wrt) x and integrate the right-hand-side wrt t:
∫ x dx = ∫ ρ ⋅ dt
o
1
(7.17)
Which gives us: ln ( x ( t )) = ρ ⋅ t + C where C is some constant. Exponentiation gives a general solution for x: x ( t ) = c ⋅ e ρ ⋅t where c = exp(C) is also a constant. (7.19) (7.18)
o
Use the initial conditions to solve for c to get specific solution for x(t): x (0) = c ⋅ e ρ ⋅0 = c , which implies that c = x (0) = x 0
Thus, the specific solution for x(t) is: (7.20)
(10)
BOT/ECOL 5756 Ecological Systems Modeling
Ogle, Fall 2006
x ( t ) = x 0 ⋅ e ρ ⋅t
o
(7.21)
Of course, you could solve this equation in Maple. You can use the dsolve function:
dsolve({ODE, ICs}, y(x))
Where ODE is the differential equation, ICs are the initial conditons, and y(x) is the function that you want to solve for. • 1st order, linear, non-homogenous ODE's.
o General procedure:
Write equation as:
dy + p( x ) ⋅ y = g ( x) dx
Find an integrating factor, m(x), where m(x) =
(7.22)
e ∫ p ( x ) dx
m( x) g ( x)dx + c Solution is: y(x) = ∫ , where c is some constant; use initial conditions to solve m( x ) for c. o Specific example:
dy − 2⋅ x⋅ y = x dx
and initial conditions: y(0) = 0
The proceedure is:
m(x) = y(x) =
e ∫ p ( x ) dx = e − ∫ 2 xdx = e − x
−x2
2
∫ (e
) ⋅ ( x)dx + c e−x
2
2
y(x) =
− 1 e−x + c 2 e
− x2
= − 1 + c ⋅e x 2
2
Employ initial conditions to solve for c:
y(0) = −
1 2
+ c ⋅ e 0 = − 1 + c ⋅1 = 0 2
2
thus, c = ½ Hence, the specific solution for y(x) is:
2 y(x) = − 1 1 − e x 2
Again, we could solve this equation in Maple.
(11)
BOT/ECOL 5756 Ecological Systems Modeling
Lecture 8
o Another example: Simple diffusion model for water vapor concentration in stomata. See the PDE in Eqn (7.14) the PDE:
∂ ∂2 w ( t , z ) = Dw ⋅ 2 w ( t , z ) ∂t ∂z
Where Dw is a diffusion coefficient (cm/s), and w = w(t,z) is the water vapor density (e.g., mol H20/cm2) at distance time t and distance z, e.g.:
Outside leaf H2O vapor L
z
0 Inside leaf
Initial conditions:
w(0,z) = wi (water vapor density inside leaf)
Boundary conditions: w(t,0) = wi (water vapor density inside leaf) w(t,L) = ws (water vapor density at the leaf surface)
Even though we won't solve the PDE, we'll learn about its steady-state behavior.
ˆ ˆ Use w to denote w(z,t) at steady-state, where w satisfies ˆ ˆ ∂w ∂2w = 0 , thus: Dw 2 = 0 . ∂t ∂z
Note that Dw is a constant (Dw ≠ 0), so we need to find the solution to:
ˆ ˆ We want to solve for w ; note that w does not depend on t
ˆ ∂ 2w ∂z 2
=0
At steady-state this problem reduces to a 2nd order, linear, homogenous ODE. Having the second derivative be zero implies that the first derivative is a constant, which ˆ means that w vs. z is a straight line at steady-state; i.e.,
ˆ w (z) = a + b⋅z
Use boundary conditions to solve for a and b.
w(0) = a = wi w(L) = a + b⋅L = wi +b⋅L = ws
which gives b = (ws-wi)/L
ˆ Thus, the SS-solution is: w (z) = wi + z⋅(ws-wi)/L
Note: transpiration rate at steady-state (Tss, mol⋅cm–2⋅s–1) is given by
Tss = − Dw
ˆ dw w − wi = − Dw s dz L
(12)
BOT/ECOL 5756 Ecological Systems Modeling
Ogle, Fall 2006
Systems of 1st order linear ODE’s
• •
Encountered in situations where there are several dependent variables (e.g., in our chaparral example we have 2 dependent (state) variables: live biomass and dead biomass). Example: Find analytical solution for the chaparral problem. Our system of equations is: From Lab 1 and the chaparral model:
d LB = p ⋅ LB − s ⋅ LB dt
Where, LB = live biomass (kg / ha) t = time (yr) s = average senescence rate (yr –1)
d DB = s ⋅ LB − d ⋅ DB dt
DB = dead (kg / ha) p = average production rate (yr –1) d = average decomposition rate (yr –1)
o Note that the equation for dLB/dt is a 1st order linear homogenous equation for LB. Following the methods above methods, we could solve for LB(t). o Note that the equation for dDB/dt is a 1st order linear non-homogenous equation for DB (the non-homogenous part results from the s⋅LB term). Once we've found the solution for LB and plugged it into the equation for dDB/dt, we can follow the methods described above to solve for DB. o Again, could solve these equations in Maple. o How about the 1st order, nonlinear, system of equations for the de Mazancourt et al. nitrogen cycle model? We can try this in Maple...
•
For compartment models that are 1st order linear constant-coefficient systems of n state variables and n differential equations, in general, for systems with 2+ state variables we can use matrix methods to find solutions (i.e., involves finding eigenvalues and eigenvectors).
Summary
•
Now we know how to find analytical solutions for: o 1st order linear homogenous ODE's o 1st order linear non-homogenous ODE's o Steady-state solutions for linear homogenous PDE's o Systems of 1st order linear ODE's Other problems that we did not discuss: o 2nd order and higher order linear equations o Nonlinear ODE's o Systems of higher order equations o PDE's
•
(13)