Introduction to Mathematical Modelling: Ordinary Differential Equations and Heat Equations Knut–Andreas Lie SINTEF ICT, Dept. Applied Mathematics
INF2340 / Spring 2005 – p.
Overview 1. Ordinary differential equations • Example of ODE models − radioactive decay, Newton’s second law, population models • Numerical solution of ODEs − forward and backward Euler, Runge-Kutta methods 2. Heat equations • from physical problem to simulator code • steady heat conduction − finite differences, linear algebra • unsteady heat conduction − finite differences, boundary conditions
INF2340 / Spring 2005 – p.
Ordinary Differential Equations (ODEs) y
(m)
= f t, y, y, ˙ y¨, . . . , y
(n−1)
0001
Simple example: radioactive decay Given a quantity q of a radioactive matter, which decays at a certain constant rate k. The model reads dq(t) = r(t) ˙ = −k · q. dt Solution: q(t) = q(t0 )e−k(t−t0 ) .
In general: finding a solution is not so easy, although there are approaches in certain special cases −→ numerical approach!
INF2340 / Spring 2005 – p.
ODEs cont’d Very simple example from high school physics: Consider the motion of a body with mass m under constant force f , which is initially at rest at position x0 . Newton’s second law reads f = ma = m¨ x(t) where x(t) is the position of the body at time t. The corresponding ODE then reads f = m¨ x(t),
˙ =0 x(0) = x0 , x(0)
This equation is easily integrated x(t) = x0 +
f 2 2m t . INF2340 / Spring 2005 – p.
Population models Consider the dynamics of a single species (isolated or with no predators) • constant birth rate b per time and individual • constant death rate d per time and individual • hence, constant growth rate λ = b − d
The model (Maltus, 1798): p(t) ˙ = λ · p(t)
Solution p(t) = p0 eλt predicts exponential growth or decay.
INF2340 / Spring 2005 – p.
Population models cont’d Is there any realism in this? • between 1700 and 1960: growth rate of about 0.02,
population doubles in 34.67 years
• generally: limited resources on earth slows down growth
More realism (Verhulst et al., 19th century): • linear rates: b(t) = b0 − b1 p(t),
d(t) = d0 − d1 p(t)
−→ new model: p(t) ˙ = −k(p(t) − p∞ )
Solution: p(t) = p∞ + (po − p∞ )e−kt
INF2340 / Spring 2005 – p.
Population models cont’d What about realism now? • Populations tend to follow a S-shape (logistic model) • New model: • New solution:
p(t) ˙ = a · p(t) − b · p2 (t) p(t) =
a·p0 b·p0 +(a−b·p0 )e−at
And so on .. adding more than one species (predator–pray) .. Do we have equilibrium, is it attractive or repellent, ..?
INF2340 / Spring 2005 – p.
Population models cont’d An ODE model from modern research, describing the dynamics of HIV-1 infection in vivo (Perelson& Nelson, SIAM Review 41/1, 1999) : The rate of change of uninfected cells T , productively infected cells T ∗ , and virus V : ´ ` dT = s + pT 1 − T /Tmax − dT T − kV T dt dT ∗ = kV T − δT ∗ dt dV = N δT ∗ − cV. dt Here: dT – death rate of uninfected cells δ – death rate of infected cells p – rate of proliferation (continuous development of cells in tissue) N – virus production per infected cell c – clearance rate INF2340 / Spring 2005 – p.
Numerical solution of ODEs – Euler’s method We wish to solve the equation: y 0 = f (y, t),
y(a) = α,
a≤t≤b
Obvious solution – use finite differences: • generate a mesh: ti = a + ih, for each i = 0, . . . , N , where
h = (b − a)/N is called stepsize
• apply forward differences to the equation
y(ti+1 ) = y(ti ) + hf (y(ti ), ti ),
i = 0, . . . , N − 1
This gives an explicit formula for each y(ti+1 ) once y0 is known.
INF2340 / Spring 2005 – p.
Example: Euler’s method for u0 = 3u 20
Exact h=1/4 h=1/16 h=1/64
18 16 14 12 10 8 6 4 2 0
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
INF2340 / Spring 2005 – p. 1
Another Euler method – backward Euler Once again we consider: y 0 = f (y, t),
y(a) = α,
a≤t≤b
and introduce a mesh: ti = a + ih, for each i = 0, . . . , N This time we apply backward differences y(ti+1 ) = y(ti ) + hf (y(ti+1 ), ti+1 ),
i = 0, . . . , N − 1
This gives an equation for each y(ti+1 ) once y0 is known.
INF2340 / Spring 2005 – p. 1
Two different methods • forward Euler: explicit method
y(ti+1 ) = y(ti ) + hf (y(ti ), ti )
the new value is given by a formula • backward Euler: implicit method
y(ti+1 ) = y(ti ) + hf (y(ti+1 ), ti+1 )
the new value is given by an algebraic equation
INF2340 / Spring 2005 – p. 1
Example: u0 = −4π sin(4πx) u(xi ) = u(xi−1 ) − 4π sin(4πxi−1 )
u(xi ) = u(xi−1 ) − 4π sin(4πxi )
Exact Forward Euler Backward Euler
1 0.8 0.6 0.4 0.2 0 −0.2 −0.4 −0.6 −0.8 −1 0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
INF2340 / Spring 2005 – p. 1
Discretisation errors When using finite differences in forward Euler, we make a local discretisation error at each point y(ti+1 ) − y(ti ) = f (y(ti ), ti ) + error h Consider y˙ = f (y). Expanding y(ti+1 ) by a Taylor polynomial: ˙ i ) + y¨(τ ) 12 h2 − y(ti ) = hf (y(ti )) + h · error y(ti ) + hy(t ˙ i ) = f (y(ti )), for ti ≤ τ ≤ ti+1 . Now since y(t error = y¨(τ ) ·
1 2h
d = f (y(τ )) · 12 h dt
We say that forward Euler is a first-order method.
INF2340 / Spring 2005 – p. 1
Discretisation errors cont’d We say that the scheme is consistent if y(t + h) − y(t) − f (y(t), t) → 0 as h → 0 l(h) = max h t∈[a,b] Similarly, we define the global discretisation error as e(h) = max yi − y(ti ) , ti ∈[a,b]
where y(ti ) is the exact solution and yi is computed by our scheme. The scheme is convergent if e(h) → 0 for t → 0
INF2340 / Spring 2005 – p. 1
Higher order methods – Runge–Kutta Consider the forward Euler method and try to evaluate f (·) at some other point yi+1 = yi + hf (yi + β). Let us repeat the error analysis y(ti+1 ) − y(ti ) = f (yi )h + f˙(yi ) 12 h2 + f¨(τ ) 16 h3 = hf (yi + β) + h · error Assume now that the error equals 16 f¨(τ )h2 . Now f˙(yi ) = f 0 (yi )f (yi ) and we have f (yi ) + f 0 (yi )f (yi ) 12 h = f (yi + β) = f (yi ) + βf 0 (yi ) + β 2 + . . . This means that β = 12 h f (yi ).
INF2340 / Spring 2005 – p. 1
Runge-Kutta methods cont’d Thus, we have a new second-order method yi+1 = yi + hf yi +
1 2 h f (yi )
0001
There are other alternatives also, e.g., yi+1 = yi +
1 2h
0002
f (yi ) + f yi + hf (yi )
00010003
The Runge-Kutta methods are all on the form:
•
Approximate the solution at a point ti ≤ τ ≤ ti+1 by the intermediate step w = yi + (τ − ti )f (yi )
•
Combine yi , w, f (yi ), f (w) to get a more accurate solution yi+1 .
For higher-order methods we use more intermediate steps.
INF2340 / Spring 2005 – p. 1
From ODEs to PDEs So far, we have used ODEs: • involve derivatives with respect to only one variable • if the unknown function depends on more variables, these
have been treated as constants
Now, partial differential equations (PDEs): • involve (partial) derivatives with respect to more than one
variable
• the unknown function depends on more than one variable • stationary problems: no time-dependence • unsteady problems. time-dependence, but maybe steady
limit
INF2340 / Spring 2005 – p. 1
Mathematical model: ∂t u = ∆u + f • Models propagation of heat within a given object • Examples − a heated rod (1D) − a heated plate or smoothing of images (2D) − heat distribution in a furnace (3D) • Functions of interest − u(x,t), u(x,y), u(x,y,t), u(x,y,z), u(x,y,z,t), .. • More generally, heat propagation depends on the
conductivity of the material
ut = ∇(K(x, u)∇u) + f (x, t, u)
INF2340 / Spring 2005 – p. 1
Heat conduction in the continental crust x=0
Ts
x x=b -Q
• Knowing the temperature at the earth’s surface and the heat flow from the mantle, what is the temperature distribution through the continental crust? • Interesting question in geology and geophysics – and for those nations exploring oil resources..
INF2340 / Spring 2005 – p. 2
Physical and mathematical model • Our prototype differential equation (an ODE in 1D): ∂2u − 2 = f (x) ∂x Here u is the temperature. • Needs to be equipped with boundary conditions • Very simple equation, but it has applications to − fluid flow in channels − deflection of electric cables − strength analysis of beams − .. • In multidimensions −∆u = f is the elliptic Poission equation, which is fequently occuring in mathematical models
INF2340 / Spring 2005 – p. 2
Heat conduction in the continental crust cont’d x=0
Ts
x x=b -Q
Physical assumptions:
• • •
Crust of infinite area Steady state heat flow Heat generated by radioactive decay
Physical quantities:
• • •
u(x) : temperature q(x) : heat flux (velocity of heat) s(x) : heat release per unit time and mass
INF2340 / Spring 2005 – p. 2
Derivation of the model x=0
inflow
s(x)=R exp(-x/L)
x
outflow
x=b
Physical principles:
•
First law of thermodynamics: net outflow of heat = total generated heat
•
Fourier’s law: heat flows from hot to cold regions (i.e. heat velocity is poportional to changes in temperature) q(x) = −λu0 (x)
•
Heat generation due to radioactive decay: s(x) = R exp(−x/L)
INF2340 / Spring 2005 – p. 2
Derivation of the model cont’d x=0 q(x-h/2) h
s(x) q(x+h/2)
x x=b
From the first law of thermodynamics: q(x + h/2) − q(x − h/2) = s(x) h Using a Taylor expansion: q(x + h/2) − q(x − h/2) 1 000 = q 0 (x) + q (x)h2 + . . . h 24 Hence, as h → 0 we have q 0 (x) = s(x)
INF2340 / Spring 2005 – p. 2
Derivation of the model cont’d Ts
x=0
x x=b -Q
Combining the 1st law of thermodynamics (q 0 = s) with Fourier’s law (q = −λu), we get d − dx
„
du λ dx
«
= s(x)
Boundary conditions:
• •
u(0) = Ts (at the surface of the earth) q(b) = −Q (at the bottom of the crust)
INF2340 / Spring 2005 – p. 2
Mathematical model d 0010 du 0011 λ = Re−x/L , − dx dx
u(0) = Ts ,
λ(b)u0 (b) = −Q
Observe that u depends upon seven parameters: u = u(x; λ, R, L, b, Ts , Q)! Suppose that we want to investigate the influence of the different parameters. Assume (modestly) three values of each parameter: −→ Number of possible combinations: 36 = 729. Using scaling we can reduce the six physical parameters λ, R, L, b, Ts , Q to only two!
INF2340 / Spring 2005 – p. 2
Scaling We introduce dimensionless quantities (and assume that λ is constant): x=x ¯b,
u/λ, u = Ts + Qb¯
s(b¯ x) = R¯ s(¯ x)
This gives ¯ d2 u − 2 = γe−¯x/β , d¯ x
u ¯ = 0,
d¯ u (1) = 1 d¯ x
where we have two dimensionless quantities β = b/L,
γ = bR/Q
Dropping the bars, we get an equation on the form −u00 (x) = f (x), x ∈ (0, 1),
u(0) = 0,
u0 (1) = 1
INF2340 / Spring 2005 – p. 2
Discretisation • Introduce a grid xi = (i − 1)h and compute the unknown at grid points ui = u(xi )
u2 u1
x=0
u3 u4
u5
x=1
x
• Differential equation fulfilled at each node −u00 (xi ) = f (xi ) • Approximate by standard finite differences ui+1 − 2ui + ui−1 = −h2 fi ,
i = 1, . . . , n − 1
As opposed to the ODEs we have seen earlier, this is a linear system of unknowns
INF2340 / Spring 2005 – p. 2
Discretising boundary conditions • u(0) = 0 simply becomes u1 = 0 • u0 (1) = 1 can be approximated as un+1 − un−1 =1 2h • Problem: un+1 is not in the mesh! • Solution: Use the discrete differential equation for i = n: un−1 − 2un + ui+1 = −h2 fn and the discrete boundary condition to eliminate un+1 • The result is 2un−1 − 2un = −2h − h2 fn
INF2340 / Spring 2005 – p. 2
Linear system of equations u1 u1
=0 −2u2
+u3
u2
−2u3 . .
= −h2 f2 +u4 . . . .
. . .
= −h2 f3 . . . . . . . .
. .
.
.
.
.
.
.
.
.
.
.
.
un−2
.
.
−2un−1
+un
2un−1
−2un
= −h2 fn−1 = −2h − h2 fn
INF2340 / Spring 2005 – p. 3
Using linear algebra We write the system as Au = b: 2 1 6 6 61 6 6 6 60 6 6 6. 6 . 6 6 6 . 6. 6 6 6 . 6. 6 6. 6. 6. 6 6. 6. 6. 4 0
0
0
0
−2
1
0
1 . .
−2 . .
.
.
.
.
.
.
.
.
.
.
···
···
···
··· . .
··· . .
··· . .
.
.
.
.
.
.
.
1 . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
···
···
···
.
1
−2
0
2
32
3
2
3
0 u1 7 6 7 76 7 6 7 76 6 7 6 7 2 0 7 6 u2 7 6 −h f2 7 7 7 6 7 76 7 6 7 76 6 7 6 7 2 0 7 6 u3 7 6 −h f3 7 7 7 6 7 76 7 6 7 . 7 6 6 7 7 6 7 . . . 76 . 7 6 7 . . 7 76 . 7 6 . 76 6 7 7 . 7 6 . 7 = 6 7 . 7 76 . 7 6 . . 7 76 . 7 6 . 76 6 7 7 . 76 . 7 6 7 . . 76 . 7 6 . 7 . 7 6 . 7 6 7 . 76 6 7 7 . 76 . 7 6 7 . . 76 . 7 6 . 7 . 7 76 . 7 6 7 6 7 76 2 6 7 6 7 1 5 4un−1 5 4 −h fn−1 7 5 −2 un −2h − h2 fn 0
This system can be solved by Gaussian elimination.
INF2340 / Spring 2005 – p. 3
Implementation Assembly of matrix in pseudocode: Matrix( real ) A(n,n) Vector(real ) b(n ), u(n) // Assemble matrix and left−hand side for i =1:n if i 1 // Left boundary A (1,1) = 1; b(1) = 0 else if i n // Right boundary A(i , i −1) = 2; A(i , i ) = −2; b( i ) = −2∗h − h∗h∗f(x(i )); else // Interior A(i , i −1) = 1; A(i , i ) = −2; A(i , i +1) = 1; b( i ) = −h∗h∗f(x(i )); end
INF2340 / Spring 2005 – p. 3
In C++ int main(int argc, char ∗∗argv) { : double ∗b, ∗u; double ∗∗A; : b = new double[n]; u = new double[n]; A = new double∗[n]; A [0] = new double[n∗n]; for ( i =1; i