Numerical Integration
Download
Report
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