Transcript Numerical Quadrature (Integration)
Numerical Differentiation and Quadrature (Integration)
d
x
,
b a
Daniel Baur
ETH Zurich, Institut für Chemie- und Bioingenieurwissenschaften ETH Hönggerberg / HCI F128 – Zürich E-Mail: [email protected]
http://www.morbidelli-group.ethz.ch/education/index
Daniel Baur / Numerical Methods for Chemical Engineers / Numerical Quadrature 1
Numerical Differentiation
Problem
: Though an analytical derivative can be found for all differentiable functions, it is often impractical to calculate it
Solution
: Approximate the derivative numerically
Method of finite differences
: Remember that: d
x
0 lim
h
0 ( 0
h h
( 0 ) Therefore: d
x
0 ( 0
h h
( 0 ) for small
h
Daniel Baur / Numerical Methods for Chemical Engineers / Numerical Quadrature 2
Numerical Quadrature (Integration)
Problem
: Generally, it is not possible to find the antiderivative (Stammfunktion) of a function f(x) in order to solve a definite integral in the interval [a,b]
Solution
: Approximate the area under the curve numerically
Trapezoidal rule
: Divide the interval into sub-intervals and approximate the integral by the sum of the areas of the resulting trapezoids Daniel Baur / Numerical Methods for Chemical Engineers / Numerical Quadrature 3
Trapezoidal rule
b a
h
n
h A
h
2
f x
1
h
2
b a
f x
1 a
h
2 x 1
h
2 x 2
f x
1 x n-1 b ( 2 ))
h
2
n i
1 1
f x i
Daniel Baur / Numerical Methods for Chemical Engineers / Numerical Quadrature
n
1 ) 4
How does Matlab do it?
quad : Low accuracy, non-smooth integrands, uses adaptive recursive Simpson rule quadl : High accuracy, smooth integrands, uses adaptive Gauss/Lobatto rule (degree of integration routine related to number of points) quadgk : High accuracy, oscillatory integrands, can handle infinite intervals and singularities at the end points, uses Gauss/Konrod rule (re-uses lower degree results for higher degree approximations) Degree
p
of an integration rule = Polynomials up to order
p
can be integrated accurately (assuming there is no numerical error) Daniel Baur / Numerical Methods for Chemical Engineers / Numerical Quadrature 5
Matlab Syntax Hints
All the integrators use the syntax result = quadl(int_fun, a, b, ...);
int_fun
is a function handle to a function that takes one input
x
returns the
function value at x
; it must be vectorized and Use parametrizing functions to pass more arguments to
int_fun
if needed f_new = @(x)int_fun(x, K);
f_new
is now a function that takes only
x
value that
int_fun
would return when
K
as input and returns the is used as second input Note that
K
must be defined in this command This can be done directly in the integrator call: result = quadl(@(x)int_fun(x, K), a, b); Matlab 2012a and newer: integral(...) Daniel Baur / Numerical Methods for Chemical Engineers / Numerical Quadrature 6
Assignment 1
1.
2.
Consider the function log d
x
0 1
x
0 Use the method of finite differences to approximate the derivative of f(x) at x = 1 . Vary h between 10 -15 and 10 -1 using logspace(-15, -1, 200) , and calculate the error of the finite differences approximation compared to the analytical solution for each h .
Plot the error vs. h using loglog . What do you observe? What could be the cause for this behavior?
Daniel Baur / Numerical Methods for Chemical Engineers / Numerical Quadrature 7
Exercise
Mass Transfer into Semi-Infinite Slab Consider a liquid diffusing into a solid material The liquid concentration at the interface is constant The material block is considered to be infinitely long, the concentration at infinity is therefore constant and equal to the starting concentration inside the block c 0 = const.
c(z, t) 0 Daniel Baur / Numerical Methods for Chemical Engineers / Numerical Quadrature z c ∞ = const.
8
Exercise (continued)
Using a local mass balance, we can formulate an ODE Accumulation Mass In Mass Out j in z j out z+ Δz d
m
d d d
c t t
V
d
c j z
d
t
j
z z
z z
d
c
d
t A j z A j z
z
where
j
is the diffusive flux in [kg m -2 s -1 ] With Δz
c
t
j z
0, we arrive at a PDE in two variables Daniel Baur / Numerical Methods for Chemical Engineers / Numerical Quadrature 9
Exercise (continued)
By combining this local mass balance with Fick’s law, a PDE in one variable is found:
j
D
d
c
d
z
c
t
D
2
c
z
2 The analytical solution of this equation (found by combination of variables) reads:
c
0
c
0 erf 2 0 exp d
s
z
4
Dt
Daniel Baur / Numerical Methods for Chemical Engineers / Numerical Quadrature 10
Assignment 2
1.
Write a Matlab program to calculate and plot the concentration profile in the slab Use the following values: c ∞ = 0; c 0 = 1 Create a vector zeta = linspace(1e-6, 3) , calculate the value of c for each zeta , then plot c vs. zeta Use integral or quadl for the integration 2.
Create a function which calculates an integral with the trapezoidal rule Use the form: function F = trapInt(f, n, a, b) Where
f n
is a function handle to the function that is to be integrated, is the number of points and
a
and
b
denote the interval Use your function to solve the diffusion problem and compare it to the result you obtained with quadl or integral .
Daniel Baur / Numerical Methods for Chemical Engineers / Numerical Quadrature 11
Assignment 2 (continued)
3.
Improve your function by adding a convergence check: In addition to computing the integral with
n
calculate it with 2
n
points points, simultaneously while the results differ by more than 10 -6 , double
n
calculation and iterate the Terminate the calculation and issue a warning if
n
exceeds 10 6 Daniel Baur / Numerical Methods for Chemical Engineers / Numerical Quadrature 12