Transcript Numerical Integration

+
Numerical Integration
Techniques
A Brief Introduction
+
2
Objectives
 Start Writing
your OWN Programs
 Make
Numerical Integration accurate
 Make
Numerical Integration fast

CUDA acceleration
3
+
The same Objective
Lord, make me accurate and fast.
- Mel Gibson, Patriot
+
4
Schedule
Methods
• Numerical Differentiation
• Euler Method
Approaches
• The Three-Variable Model
• Heat Diffusion Equation
+
5
Preliminaries
 Basic


Calculus
Derivatives
Taylor series expansion
ò w = ò dw
 Basic

Programming Skills
Octave
+
6
Numerical Differentiation

Definition of Differentiation
df
f (x + h) - f (x)
= f '(x) = lim
h®0
dx
h

Problem: We do not have an infinitesimal h

Solution: Use a small h as an approximation
+
7
Numerical Differentiation
Forward Difference

Approximation Formula
df
f (x + h) - f (x)
» f '(x)approx =
dx
h

Is it accurate?
+
8
Numerical Differentiation
Error Analysis

Taylor Series expansion uses an infinite sum of terms to
represent a function at some point accurately.
1 2
1 3
f (x + h) = f (x) + hf '(x) + h f ''(x) + h f '''(x) +…
2
6

which implies
f (x+ h) - f (x)
1
1 2
f '(x)approx =
= f '(x) + hf ''(x) + h f '''(x)+…
h
2
6
9
+
Truncation Error
1
1 2
et = hf ''(x) + h f '''(x)+… = O(h)
2
6
+
10
Numerical Differentiation
Error Analysis

Roundoff Error
 A computer can not store a real number in its memory
accurately.
 Every number is stored with an inaccuracy proportional to
itself.
 Denoted er

Total Error
e= et + er
Usually we consider Truncation Error more.
+
11
Numerical Differentiation
Backward Difference

Definition
f (x) - f (x - h)
f '(x)approx =
h

Truncation Error
1
1
et = - hf ''(x) + h2 f '''(x) +… = O(h)
2
6

No Improvement!
+
12
Numerical Differentiation
Central Difference

Definition
f (x + h) - f (x - h)
f '(x)approx =
2h

Truncation Error
1 2
1 4 (5)
et = h f '''(x) +
h f (x) +… = O(h2 )
6
120

More accurate than Forward Difference and Backward
Difference
+
13
Numerical Differentiation
Example

Compute the derivative of function
f (x) = ex

At point x=1.15
f '(1.15) = f (1.15) » 3.1582
14
+
Use Octave to compare these methods
Blue – Error of Forward Difference
Green – Error of Backward Difference
Red – Error of Central Difference
+
15
Numerical Differentiation
Multi-dimensional & High-Order
 Multi-dimensional

Apply Central Difference for different parameters
d2 f d df f (x+ 2h)+ f (x - 2h) - 2 f (x)
=
»
2
dx
dx dx
4h2
 High-Order

Apply Central Difference several times for the same
parameter
¶2 f
¶ ¶f f (x + h, y+ h) + f (x - h, y- h) - f (x+ h, y- h) - f (x - h, y+ h)
=
»
¶x¶y ¶y ¶x
4h2
+
16
Euler Method
IVP
 The


Initial Value Problem
Differential Equations
Initial Conditions
ì y'(t) = f (t, y(t))
í
î y(t0 ) = y0
 Problem

What is the value of y at time t?
+
17
Euler Method
Explicit Euler Method
 Consider
Forward Difference
y(t + Dt) - y(t)
y'(t) »
Dt
 Which
implies
y(t + Dt) » y(t)+ Dt ×y'(t)
+
18
Euler Method
Explicit Euler Method
 Split
time t into n slices of equal length Δt
×t0 = 0
×
×ti = i ×Dt
×t = t
×n
 The
Explicit Euler Method Formula
y(ti+1 ) = y(ti )+ Dt ×y'(ti )
+
19
Euler Method
Explicit Euler Method - Algorithm
+
20
Euler Method
Explicit Euler Method - Error

Using Taylor series expansion, we can compute the
truncation error at each step
y(t + Dt) = y(t)+ Dt ×y'(t)+O(Dt )
2
et-step = O(Dt )
2

We assume that the total truncation error is the sum of
truncation error of each step
t
et = ×O(Dt 2 ) = O(Dt)
Dt

This assumption does not always hold.
+
21
Euler Method
Implicit Euler Method
 Consider
Backward Difference
y(t) - y(t - Dt)
y'(t) »
Dt
 Which
implies
y(t) » y(t - Dt)+ Dt ×y'(t)
+
22
Euler Method
Implicit Euler Method
 Split
the time into slices of equal length
y(ti+1 ) » y(ti )+ Dt ×y'(ti+1 )
 The
above differential equation should be solved
to get the value of y(ti+1)
Extra computation
 Sometimes worth because implicit method is more
accurate

+
23
Euler Method
A Simple Example
 Try
to solve IVP
ì y'(t) = e-t + t
í
î y(0) = 1
 What
 The
is the value of y when t=0.5?
analytical solution is
1 2
y = -e + t + 2
2
-t
+
24
Euler Method
A Simple Example
 Using
explicit Euler method
-ti
yi+1 = yi + dt × (e + ti )
 We
choose different dts to compare the accuracy
ì dt1 = 0.1 Þ t = 0, 0.1, 0.2,..., 0.5
ï
í dt2 = 0.05 Þ t = 0, 0.05, 0.1,..., 0.5
ï dt = 0.01 Þ t = 0, 0.01, 0.02,..., 0.5
î 3
+
25
Euler Method
A Simple Example
t
exact
dt=0.05
error
dt=0.025
error
dt=0.012
5
error
0.1
0.2
0.3
0.4
0.5
1.10016
1.20126
1.30418
1.40968
1.51846
1.10030
1.20177
1.30525
1.41150
1.52121
0.00014
0.00050
0.00107
0.00182
0.00274
1.10022
1.20151
1.30470
1.41057
1.51982
0.00006
0.00024
0.00052
0.00089
0.00135
1.10019
1.20138
1.30444
1.41012
1.51914
0.00003
0.00011
0.00025
0.00044
0.00067
At some given time t, error is proportional to dt.
+
26
Euler Method
Instability
 For
some equations called Stiff Equations, Euler
method requires an extremely small dt to make the
result accurate
y'(t) = -k× y(t), k > 0
 The
Explicit Euler Method Formula
yi+1 = yi - Dt ×k×yi = (1- Dt ×k)yi
 The
choice of Δt matters!
+
27
Euler Method
Instability

Assume k=5
ì y'(t) = -5y(t)
í
î y(0) =1

Analytical Solution is
-5t
y(t) = e

Try Explicit Euler Method with different dts
28
+ Choose dt=0.002, s.t.
Works!
1
0 <1- Dt ×k <1 × Dt <
k
29
+ Choose dt=0.25, s.t.
1
2
-1 <1- Dt ×k < 0 × < Dt <
k
k
Oscillates, but works.
30
+ Choose dt=0.5, s.t.
Instability!
2
Dt > × 1- Dt ×k < -1
k
+
31
Euler Method
Stiff Equation – Explicit Euler Method

For large dt, explicit Euler Method does not guarantee an
accurate result
t
exact
dt=0.5
error
dt=0.25
error
dt=0.002
error
0.4
0.135335
1
6.389056
-0.25
2.847264
0.13398
0.010017
0.8
0.018316
-1.5
-0.015625 1.853096
0.017951
0.019933
1.2
0.002479
2.25
-0.000977 1.393973
0.002405
0.02975
1.6
0.000335
-3.375
-0.000061 1.181943
0.000322
0.039469
2
0.000045
5.0625
82.897225
906.71478
5
10061.733
21
111507.98
31
0.000015
0.000043
0.04909
0.663903
+
32
Euler Method
Stiff Equation – Implicit Euler Method
 Implicit
Euler Method Formula
yi+1 = y i - dt × 5× y i+1
 Which
implies
yi
yi+1 =
1+ 5dt
33
+
Choose dt=0.5,
Oscillation eliminated!
Not elegant, but works.
+
34
The Three-Variable Model
Single Cell
 The
Differential Equations
ѶtV = Ñ Ñ(DÑV) - (I fi + I so + I si ) + I ext
Ñ
Ѷt v = (1- p)(1- v) / t v- - pv / t v+
Ñ
+
¶
w
=
(1p)(1w)
/
t
pw
/
t
w
w
Ñt
I fi = -vp(V -Vc )(V -Vm ) / t d
I so = (V -Vo )(1- p) / t o + p / t r
I si = -w(1+ tanh[k(V -Vcsi )]) / (2t si )
ì1 if V ³ Vc
p=í
î0 if V < Vc
ì1 if V ³ VV
q=í
î0 if V < VV
+
35
The Three-Variable Model
Single Cell
 Simplify
the model
ì¶tV = -(I fi + I so + I si ) + I ext
ï
í¶t v = (1- p)(1- v) / [1000q+19.2(1- q)]- 0.3pv
ï¶ w = (1- p)(1- w) /11- pw / 667
î t
I fi = 0.04vp(V + 72)(15 -V)
I so = (V + 85)(1- p) / 8.3+ 2 p
I si = -w(1+ tanh[0.1V]) / 0.897
ì1 if V ³ -72
ì1 if V ³ -79.5
p=í
q=í
î0 if V < -72
î0 if V < -79.5
+
36
The Three-Variable Model
Single Cell
 Using




explicit Euler method
Select simulation time T
Select time slice length dt
Number of time slices is nt=T/dt
Initialize arrays vv(0, …, nt), v(0, …, nt), w(0, …, nt) to store
the values of V, v, w at each time step
+
37
The Three-Variable Model
Single Cell
 At
each time step
Compute p,q from value of vv of last time step
 Compute Ifi, Iso, Ifi from values of vv, v, w of previous
time step

I fi = 0.04 × v(i -1)× p× (vv(i -1) + 72)× (15 - vv(i -1))
I so = (vv(i -1) + 85)× (1- p) / 8.3+ 2 p
I si = -w(i -1)× (1+ tanh[0.1× vv(i -1)]) / 0.897
+
38
The Three-Variable Model
Single Cell
 At

each time step
Use explicit Euler method formula to compute new vv,
v, and w
vv(i) = vv(i -1) + dt × (-(I fi + I so + I si ) + I ext )
v(i) = v(i -1) + dt × (1- p)× (1- v(i -1)) / (1000q+19.2(1- q)) - 0.3p× v(i -1)
w(i) = w(i -1) + dt × (1- p)× (1- w(i -1)) /11- p× w(i -1) / 667
39
+
The Three-Variable Model
+
40
Heat Diffusion Equations

The Model

The first equation describes the heat conduction



׶u
× = aDu
× ¶t
×u = f
× t=0
Function u is temperature distribution function
Constant α is called thermal diffusivity
The second equation initial temperature distribution
+
41
Heat Diffusion Equation
Laplace Operator

Laplace Operator Δ (Laplacian)

Definition: Divergence of the gradient of a function
Df = × 2 f = × ×(×f )

Divergence measures the magnitude of a vector field’s
source or sink at some point
div( f ) = Ñ Ñf

Gradient is a vector point to the direction of greatest
rate of increase, the magnitude of the gradient is the
greatest rate
Ñf
+
42
Heat Diffusion Equation
Laplace Operator
 Meaning
of the Laplace Operator
Du = × ×(×u)
= × ×(thegradient of temperaturedistribution)
= themagnitudeof temperaturechangeover space
 Meaning

of Heat Diffusion Operator
At some point, the temperature change over time equals
the thermal diffusivity times the magnitude of the greatest
temperature change over space
+
43
Heat Diffusion Equation
Laplace Operator
 Cartesian
coordinates
x1, x2 ,… , xn
¶2 f
Df = × 2
i=1 ¶xi
n
 1D
space (a cable)
¶2 f
Df = 2
¶x
+
44
Heat Diffusion Equation
Laplace Operator
 Compute
 Similar
Laplacian Numerically (1d)
to Numerical Differentiation
¶2 f f (x+ dx) + f (x - dx) - 2 f (x)
Df = 2 =
¶x
dx2
+
45
Heat Diffusion Equation
Laplace Operator
 Boundaries
 Assume
(1d), take x=0 for example
the derivative at a boundary is 0
¶f
=0
¶x x=0
f (Dx) - f (-Dx)
×
=0
2Dx
× f (Dx) = f (-Dx)
 Laplacian
at boundaries
Df
=
x=0
f (Dx) + f (-Dx) - 2 f (0) 2 f (Dx) - 2 f (0)
=
Dx2
Dx2
+
46
Heat Diffusion Equation

Exercise

Write a program in Octave
to solve the following heat
diffusion equation on 1d
space:
׶u
× ¶t = aDu
×
×
px
u(x,
0)
=
10sin(
)
×
10
×
×0 × x ×10
×a = 0.9
×
+
47
Heat Diffusion Equation

Exercise

TIPS:

Write a program in Octave
to solve the following heat
diffusion equation on 1d
space:

Store the values of u in a 2d
array, one dimension is the
time, the other is the
cable(space)

Use explicit Euler Method

Choose dt, dx carefully to
avoid instability

You can use mesh()
function to draw the 3d
graph
׶u
× ¶t = aDu
×
×
px
u(x,
0)
=
10sin(
)
×
10
×
×0 × x ×10
×a = 0.9
×
48
+
Store u in a two-dimensional array
u(i,j) stores value of u at point x=xi, time t=tj.
u(i,j) is computed from u(i-1,j-1), u(i,j-1), and u(i+1,j+1).
+
49
Heat Diffusion Equation

Explicit Euler Method Formula
u(x, t) = u(x,t - dt)+ dt ×aDu

Discrete Form
ì
u(i -1, j -1) + u(i +1, j -1) - 2u(i, j -1)
ïu(i, j ) = u(i, j -1) + dt × a
2
dx
ï
ï
2u(1, j -1) - 2u(0, j -1)
íu(0, j ) = u(0, j -1) + dt × a
2
dx
ï
ï
2u(n-1, j -1) - 2u(n, j -1)
u(n,
j
)
=
u(n,
j
-1)
+
dt
×
a
ïî
dx2
+
50
Heat Diffusion Equation

Stability
2dt × a
0 < 1<1
2
dx
Þ 2dt × a < dx2

dt must be small enough to avoid instability
+
The END
Thank You!
51