Numerical Methods for Partial Differential Equations
Download
Report
Transcript Numerical Methods for Partial Differential Equations
Numerical Methods for Partial
Differential Equations
CAAM 452
Spring 2005
Lecture 9
Instructor: Tim Warburton
Today
• Introduction to a very basic finite volume method
CAAM 452 Spring 2005
Mass Conservation
• For the following we are going to consider a 1-dimensional
domain which we parameterize with the variable x.
x
• Now imagine that this line is a figurative representation of
a pipe which contains fluid.
• At every point on the line the fluid has a density measured
in mass per meter ( with units kgm-1 )
• We define a new function : [0, ) which is a
non-negative, real valued function defined on the space-time
domain.
CAAM 452 Spring 2005
Derivation of Mass Conservation Law
• Next we consider an arbitrary section of the pipe, say [a,b]
• We now assume that the fluid is not created or destroyed at any point
inside the section and is traveling with velocity u (which is a function of
space and time). For the moment we will assume that u is positive (i.e.
the fluid is flowing in the direction of positive x)
• This allows us to state the following:
– The time rate of change of the total fluid inside the section [a,b]
changes only due to the flux of fluid into and out of the pipe at the
ends x=a and x=b.
• A simple formula relating these two quantities is:
b
d
x, t dx u b, t b, t u a, t a, t
dt a
CAAM 452 Spring 2005
• In detail:
b
d
x, t dx u b, t b, t u a, t a, t
dt a
The time rate of
change of total
mass in the
section of pipe
[a,b]
The flux out of
the section at
the right end of
the section of
pipe per unit
time
The flux into the
section at the
left end of the
section of pipe
per unit time
CAAM 452 Spring 2005
Use the Fundamental Theorem of
Calculus
• Look carefully at the right hand side:
b
d
x, t dx u b, t b, t u a, t a, t
dt a
• Clearly we may rewrite this as:
d
x
,
t
dx
u x, t x, t dx
dt a
x
a
b
b
• From which
we may deduce:
b
a t x, t x u x, t x, t dx 0
CAAM 452 Spring 2005
Finally…
• Assuming that the integrand of:
x
,
t
a t x u x, t x, t dx 0
b
is continuous and noting that this relation holds for
all choices of a,b then we may deduce:
x, t u x, t x , t 0
t
x
• In short hand:
u
0
t
x
CAAM 452 Spring 2005
Advection Equation
• Let’s choose a simple, constant, fluid velocity
u x, t u
• Then the pde reduces to the advection equation:
u
0
t
x
• This is a pretty easy equation to solve . Consider
the change of variables: t t
x x ut
CAAM 452 Spring 2005
• Note that the units of each variable are
consistent. t t
x x ut
• Basic calculus:
t x
u
t t t t x t
x
t x
x x t x x x
• From which we obtain:
u
t
x
u
u
x
t
x
t
0
CAAM 452 Spring 2005
Solution and Interpretation
• So we know: 0
t
• Which we can instantly solve: x, t 0 x
0 x ut
where:
0 x : x, t 0
• So an interesting property of the advection
equation is the way that the profile of the
solution does not change shape but it does shift
in the positive x direction with constant velocity
CAAM 452 Spring 2005
Space Time Diagram
• Let’s track the information:
t
1
Slope =
u
x
• The dashed lines are x ut const which are known as
characteristics of the equation.
• If we choose a point on one of these dashed lines and track
back down to t=0 and we will find the value of the density which
applies at all points on the dashed line
CAAM 452 Spring 2005
Recall
• Assuming that the integrand of:
b
a t x, t x u x, t x, t dx 0
is continuous and noting that this relation holds for all
choices of a,b then we may deduce:
x, t u x, t x , t 0
t
x
• Well – this does not hold if the density is discontinuous
and the integral equation is the appropriate representation:
b
d
x, t dx u b, t b, t u a, t a, t
dt a
CAAM 452 Spring 2005
Building a Finite Volume Solver
1) Let’s consider the advection equation:
b
d
x, t dx u b, t u a, t
dt a
2) Next we take a finite portion of the real line from
x1 to xN divided into N-1 equal length sections
x1
dx
xN
3) In each section we will approximate the density by a
constant value
i
i 1,..., N 1
CAAM 452 Spring 2005
Piecewise Constant Approximation
x1
xN
In the n’th section the density will be approximated by the
constant:
n
CAAM 452 Spring 2005
• Choosing a = xi, b=xi+1,
b
d
x, t dx u b, t u a, t
dt a
• Is approximated (to first order in time) by:
n1
n
i i dx
u xi1 , t u xi , t
dt
where the time axis has been divided into sections of length
dt and the i’th cell average at time n*dt is represented by
in
xi 1
x, t ndt dx
xi
• Outstanding question: Now we have to figure out how to
evaluate the density at the interval end points given the cell
averages.
CAAM 452 Spring 2005
Upwind Treatment for Flux Terms
t
1
Slope =
u
• Recall that the solution shifts from left to right
as time increases.
• Idea: use the upwind values
u xi1 , t u xi , t u in u in1
CAAM 452 Spring 2005
Basic Upwind Finite Volume Method
n1
i
in dx
dt
u in u in1
simplify
n 1
i
dt n dt n
1 u i u i 1
dx
dx
u
dt
dx
in1 1 in in1
Note we must supply a value for the left most average at each time step:
n
0
CAAM 452 Spring 2005
Recall
Basic Upwind Finite Volume Method
d
qi dx uq xi1 , t uq xi , t
dt
Approximate fluxes with upwind flux
d
qi dx uqi uqi1
dt
n1
i
in dx
dt
Approximate time derivative and look for
solution in qin
u in u in1
u
dt
dx
in1 1 in in1
Note we must supply a value for the left most average at each time step:
n
0
CAAM 452 Spring 2005
Convergence
• We have constructed a physically reasonable numerical scheme to
approximate the advection equation.
• However, we need to do some extra analysis to determine how
good at approximating the true PDE the discrete scheme is.
• Let us suppose that the i’th subinterval cell average of the actual
solution to the PDE at time T=n*dt is denoted by
qin
1
xi 1 xi
xi 1
q x, ndt
xi
xi 1
1
q x, ndt
dx xi
where q satisfies:
d
dt
xi 1
xi
d
q x, t dx dxqi uq xi 1 , t uq xi , t
dt
CAAM 452 Spring 2005
Error Equation
• The goal is to estimate the difference of the exact solution
and the numerically obtained solution at some time
T=n*dt.
• So we are interested in the error:
E q
n
i
n
i
n
i
T
,n
dt
• For the given finite volume scheme dt and dx will be
related in a fixed manner (i.e. dt = Cdx for some C,
independent of dx).
• Suppose we let
dt 0
and
E n O dt s
then
the scheme is said to be of order s.
CAAM 452 Spring 2005
Norms and Definitions
• We define the discrete p-norms:
E
p
1/ p
p
dx Ei
i
i
• We say that the scheme is convergent at time T in
the norm ||.|| if:
n
lim E 0
dt 0
ndt T
• It is said to be accurate of order s if:
E n O dt s as dt 0
CAAM 452 Spring 2005
Local Truncation Error
• Suppose at the beginning of a time step we actually have the exact
solution -- one question we can ask is how large is the error
committed in the evaluation of the approximate solution at the end
of the time step.
n
n
i qi
dt n dt n
1 u qi u qi 1
dx
dx
1
Rin : in1 qin1
dt
n 1
i
• i.e. choose
• Then
• We expand
n
i 1
n1
about
i
q ,q
xi,tn with Taylor series:
____
2
q
dx
q O (dx 3 )
qin1 qin dx
2
x 2 x
____
____
2 2
q dt q
qin 1 qin dt
O (dt 3 )
2
t 2 t
____
2
CAAM 452 Spring 2005
Estimating Truncation Error
• Inserting the formulas for the expanded q’s:
Rin :
1 dt n dt n
n 1
1
u
q
u
q
q
i
i 1 i
dt dx
dx
1 dt u q n
i
dx
____
____
2
2
1
dt
q dx q
3
u qin dx
O
(
dx
)
2
dt dx
x 2 x
____
____
2
2
n
q dt q
3
O (dt )
qi dt
2
t 2 t
CAAM 452 Spring 2005
Estimating Truncation Error
• Removing canceling terms:
1 dt u q n
i
dx
____
____
2
2
1
dt
q dx q
3
Rin : u qin dx
O
(
dx
)
2
dt dx
x 2 x
____
____
2
2
n
q dt q
3
O (dt )
qi dt
2
t 2 t
CAAM 452 Spring 2005
Estimating Truncation Error
• Simplifying:
____
____
2
2
q dx q
dt
3
dx u dx x 2 x 2 O(dx )
1
Rin :
____
____
dt
2
2
q
dt
q
3
dt
O
(
dt
)
2
t 2 t
____
____
q u q
t x
n
Ri :
____
____
dt 2 q
dx 2 q
2
2
2 O(dt ) u 2 O(dx )
2 x
2 t
CAAM 452 Spring 2005
Final Form
____
____
definition
of q
• Using the
q
q
u
t x
n
Ri :
____
____
2
dt 2 q
dx
q
2
2
2 O(dt ) u 2 O(dx )
2 x
2 t
____
____
2
2
dt q
dx q
Rin : 2 O(dt 2 ) u 2 O(dx 2 )
2 t
2 x
____
____
2
2
q
Using that: q
2
u
t 2
x 2
____
2
udx udt q
2
Rin :
1
O
(
dx
)
2
2
dx x
CAAM 452 Spring 2005
Interpretation of Consistency
____
2
udx
udt
q
n
2
Ri :
1
O
(
dx
)
2
2
dx x
So the truncation error is O(dx) under the assumption that
dt/dx is a constant..
This essentially implies that the numerical solution diverges
from the actual solution by an error of O(dx) every time step.
If we assume that the solution q is smooth enough then the
truncation error converges to zero with decreasing dx. This
property is known as consistency.
CAAM 452 Spring 2005
Error Equation
in qin Ein
• We define the error variable:
n 1
n
N
• We next define the numerical iterator N:
• Then:
E
n 1
N q E
q
N q E N q N q q
N q E N q dtR
n
n
n
n
n
n 1
n
n
n
n
n 1
n
• So the new error consists of the action of the numerical
scheme on the previous error and the error commited in the
approximation of the derivatives.
CAAM 452 Spring 2005
Abstract Scheme
• Without considering the specific construction of the scheme
suppose that the numerical N operator satisfies:
N P N Q P Q
• i.e. N is a contraction operator in some norm then…
CAAM 452 Spring 2005
Estimating Error in Terms of Initial Error
and Cumulative Truncation Error
E n1 N q n E n N q n dtR n
E n1 N q n E n N q n dt R n
triangle inequality
E n1 q n E n q n dt R n
contraction property of N
E n1 E n dt R n
E n1 dt R n1 dt R n
.....
E
n 1
mn
E dt R m
0
by induction
m 1
CAAM 452 Spring 2005
Error at a Time T (independent of
dx,dt)
E
n 1
mn
E dt R m
0
m 1
E n1 E 0 T max R m
m 1,..n
• If the method is consistent (and the actual solution
is smooth enough) then:
E n1 E 0 T O dx
• As dx -> 0 the initial error ->0 and consequently
the numerical error at time T tends to zero with
decreasing dx (and dt).
CAAM 452 Spring 2005
Specific Case: Stability and Consistency
for the Upwind Finite Volume Scheme
• We already proved that the upwind FV scheme is
consistent.
• We still need to prove stability of.
in1 1 in in1
CAAM 452 Spring 2005
Stability in the Discrete 1-norm
in1 1 in in1
N 1
n 1
1
dx in1
i 1
N 1
dx 1 in in1
triangle inequality
i 1
N 1
N 1
dx 1 dx in1 assuming 0 1
n
i
i 1
i 1
N 1
dx dx in
n
0
i 1
dx 0n n
1
So here’s the interesting story. In the case of a zero boundary
condition then we automatically observe that the operator is a
contraction operator.
CAAM 452 Spring 2005
Boundary Condition
• Suppose , are two numerical solutions with
n
0
• Then:
n
0
n1 n1 1 dx 0n 0n n n
n n
1
1
• i.e. if we are spot on with the left boundary
condition the N iterator is indeed a contraction.
CAAM 452 Spring 2005
Relaxation on Stability Condition
• Previous contraction condition on the numerical
iterator N
N P N Q P Q
• A less stringent condition is:
N P N Q 1 dt P Q
• Where alpha is a constant independent of dt as dt>0
CAAM 452 Spring 2005
Relaxation on Stability Condition
• In this case the stability analysis yields:
E n 1 N q n E n N q n dtR n
E n 1 q n E n q n dt R n
E n 1 1 dt E n dt R n
1 dt 1 dt E n 1 dt R n 1 dt R n
.....
E
n 1
1 dt
eT
n 1
mn
E dt 1 dt
0
m 1
E 0 T max R m
m 1,..n
nm
Rm
Ndt=T
CAAM 452 Spring 2005
Interpretation
E n1 eT E 0 T max R m
m1,..n
, Ndt=T
• Relaxing the stability yields a possible exponential
growth – but this growth is independent of T so if
we reduce dt (and dx) then the error will decay to
zero for fixed T.
CAAM 452 Spring 2005
Homework 2
• Due on Friday
• Recall report presentation policy
• Time will be given in class on Wednesday to go
over questions.
CAAM 452 Spring 2005
Homework 2
Q1) Solve the pde
2
0
t
x
( x, t 0) e
Q2) Solve the pde
analytically.
Q3) Solve the pde
analytically
analytically on the domain
x2
, [0, )
2
x 2t
t
x
( x, t 0) e
x2
3
t
x
( x, t 0) e
x2
and explain what happens to the density along the
characteristics.
PTO
CAAM 452 Spring 2005
Homework 2 cont
Q4) Implement the numerical approximation of:
3
0
t
x
( x, t 0) e
x2
Geometry:
By:
N i
i 1
x1 4, xN 4, xi
x
1
xN
N 1
N 1
Scheme:
in1 1 in in1
dt
3
dx
Initial Condition :
xi xi 1
,t 0
2
Boundary Condition :
0
i
00 0
CAAM 452 Spring 2005
Homework 2 cont
Q4 cont)
For N=10,40,160,320,640,1280 run to t=100, with:
dx = (8/(N-1))
dt = dx/6
On the same graph, plot t on the horizontal axis and error
on the vertical axis. The graph should consist of a sequence
of 6 curves – one for each choice of dx.
Comment on the curves.
NOTE: For the purposes of this test we define error as:
x x
error n max in i i 1 , ndt
1i N
2
CAAM 452 Spring 2005
Lecture 5
• We will define stability for a numerical scheme and
investigate stability for the upwind scheme.
• We will compare this scheme with a finite difference
scheme.
• We will consider alternative ways to approximate
the flux functions.
CAAM 452 Spring 2005
Homework 2
correction
Q4) Implement the numerical approximation of:
3
0
t
x
( x, t 0) e
x2
Geometry:
By:
N i
i 1
x1 4, xN 4, xi
x
1
xN
N 1
N 1
Scheme:
in1 1 in in1
dt
3
dx
Initial Condition :
xi xi 1
,t 0
2
Boundary Condition :
0
i
00 0
CAAM 452 Spring 2005
Homework 2 cont
Q4 cont)
For N=10,40,160,320,640,1280 run to t=10, with:
dx = (8/(N-1))
dt = dx/6
Use
0n1 1 u
dt n dt n
0 u N 1
dx
dx
On the same graph, plot t on the horizontal axis and error on the vertical
axis. The graph should consist of a sequence of 6 curves – one for each
choice of dx.
Comment on the curves.
NOTE: For the purposes of thistest we define
error as:
xi xi 1
n
n
error max i
, ndt
1i N
2
CAAM 452 Spring 2005