Transcript 슬라이드 1 - POSTECH, Bio-Medical Electrics Laboratory

How to solve PDEs using
MATHEMATIA and MATLAB
2006. 5. 17
G. Y. Park, S. H. Lee and J.K. Lee
Department of Electronic and Electrical Engineering,
POSTECH
Plasma Application
Modeling POSTECH
Contents
• Solving PDEs using MATHEMATICA
- FTCS method
- Lax method
- Crank Nicolson method
- Jacobi’s method
- Simultaneous-over-relaxation (SOR) method
• Solving PDEs using MATLAB
- Examples of PDEs
Plasma Application
Modeling POSTECH
References
• Textbook
- ‘Numerical and Analytical Methods for Scientists and
Engineers Using Mathematica’, Daniel Dubin, Wiley, 2003
- ‘Applied numerical methods in C’, S. Nakamura
Plasma Application
Modeling POSTECH
PDE (Partial Differential Equation)
PDEs are used as mathematical models for phenomena in
all branches of engineering and science.
- Three Types of PDEs:
1) Elliptic:
   (cu)  au  f
 Steady heat transfer, flow and diffusion
u
2) Parabolic: d t    (cu )  au  f
 Transient heat transfer, flow and diffusion
 2u
3) Hyperbolic: d t 2    (cu )  au  f
 Transient wave equation
Plasma Application
Modeling POSTECH
FTCS method for the heat equation
Heat equation in a slab

2
T ( x, t )   2 T ( x, t )  S ( x, t )
t
x
(T (0, t )  T1 (t ),T ( L, t )  T2 (t ),  : therm aldiffusivity)
FTCS ( Forward Euler in Time and Central difference in Space )
Let T ( x j , tn )  Tjn
T jn 1  T jn

T ( x j , tn ) 
t
t
T jn1  2T jn  T jn1
2
T ( x j , tn ) 
2
x
x 2
T jn 1  T jn  tS nj 
t
x 2
(T jn1  2T jn  T jn1 )
Plasma Application
Modeling POSTECH
FTCS method for the heat equation
FTCS
Initial
conditions
Plot
t  0.01
Stability of FTCS and CTCS
FTCS is first-order accuracy in time and second-order accuracy in space.
So small time steps are required to achieve reasonable accuracy.
t  0.05
Courant condition for FTCS
( unstable )
1

2
x
2
t
CTCS method for heat equation
(Both the time and space derivatives
are center-differenced.)
T jn 1  T jn 1  2tS nj 
2 t n
n
n
(
T

2
T

T
j 1
j
j 1 )
x 2
However, CTCS method is unstable
for any time step size.
Plasma Application
Modeling POSTECH
Lax method
Simple modification to the CTCS method
In the differenced time derivative,
Replacement by average value from surrounding grid points
(T jn1  T jn1 )
2t

[T jn1  12 (T jn11  T jn11 )]
2t
The resulting difference equation is
2 t n
n 1
n 1
n 1
n
n
n
1
T j  2 (T j 1  T j 1 )  2tS j 
(
T

2
T

T
j 1
j
j 1 )
2
x
Courant condition for Lax method
t
1
( Second-order accuracy in both time and space )

2
x
4
Plasma Application
Modeling POSTECH
Crank Nicolson Algorithm ( Implicit Method )
BTCS ( Backward time, centered space ) method for heat equation
T jn  T jn1  tS nj  x2t (T jn1  2T jn  T jn1 )
a set of coupled linear equations for
T jn
( This is stable for any choice of time steps,
however it is first-order accurate in time. )
Crank-Nicolson scheme for heat equation
taking the average between time steps n-1 and n,
T jn  T jn1 
1
2
t (S
n
j
 S nj1 )
 x2t (T jn1  2T jn  T jn1  T jn11  2T jn1  T jn11 )

( This is stable for any choice of time steps and second-order accurate in time. )
Plasma Application
Modeling POSTECH
Crank Nicolson Algorithm
Initial
conditions
CrankNicolson
scheme
Plot
Exact solution
T ( x, t )  e
  2t / L2
sin(x / L)
Crank Nicolson Algorithm
Plasma Application
Modeling POSTECH
Multiple Spatial Dimensions
FTCS for 2D heat equation
 2

2 
T ( x, y, t )    2  2 T ( x, y, t )  S ( x, y, t )
t
 x y 
 T jkn1  T jkn  tS njk 
t
x 2
(T jn1k  2T jkn  T jn1k ) 
t
y
n
n
n
(
T

2
T

T
jk

1
jk
jk
1 )
2
Courant condition for this scheme
1 x 2 y 2
t 
2 x 2  y 2
( Other schemes such as CTCS and Lax can be easily extended to multiple dimensions. )
Plasma Application
Modeling POSTECH
Wave equation with nonuniform wave speed
2D wave equation
 2
2
2 
2
z ( x, y, t )  c ( x, y) 2  2  z ( x, y, t )
2
t
 x y 
Initial condition :
Boundary
condition :
z ( x, y,0)  z0 ( x, y), z( x, y,0)  v0 ( x, y)
z0  v0  0
z (0, y, t )  z (1, y, t )  z ( x,1, t )  0,
z ( x,0, t )  x(1  x) sin 2t
2
2
Wave speed : c( x, y)  101 1  0.7 exp[2( x  0.5)  5( y  0.5) ]
CTCS method for the wave equation :
c 2jk t 2 n
c 2jk t 2 n
n 1
n
n 1
n
n
n
n
z jk  2 z jk  z jk 
(
z

2
z

z
)

(
z

2
T

T
j 1k
jk
j 1k
jk 1
jk
jk 1 )
x 2
y 2
c 2 t 2 1
 ( for x  y   )
Courant condition :
2

2
Plasma Application
Modeling POSTECH
Wave equation with nonuniform wave speed
Since evaluation of the nth timestep
refers back to the n-2nd step,
for the first step, a trick is employed.
z1jk  z 0jk  v0 jk t  2t a jk
2
a jk 
c 2jk
x
2
(z
0
j 1k
 2z  z
0
jk
0
j 1k
)
c 2jk
y
2
( z njk 1  2T jkn  T jkn 1 )
Since initial velocity and value,
v0 jk  z 0jk  0
 z1jk  0
Plasma Application
Modeling POSTECH
Wave equation with nonuniform wave speed
Plasma Application
Modeling POSTECH
Wave equation with nonuniform wave speed
Plasma Application
Modeling POSTECH
2D Poisson’s equation
Poisson’s equation
(
2
x 2

2
y 2
) ( x, y )   ( x, y )
Centered-difference the spatial derivatives
 j 1k  2 jk   j 1k
x
2

 jk 1  2 jk   jk 1
y
2
Direct Solution for Poisson’s equation
  jk
Jacobi’s method ( Relaxation method )
Direct solution can be difficult to program efficiently.
Relaxation methods are relatively simple to code,
however, they are not as fast as the direct methods.
Idea :
Poisson’s equation can be thought of as the equilibrium solution to the
heat equation with source.
II. Starting with any initial condition, the heat equation solution will
eventually relax to a solution of Poisson’s equation.
I.
(Maximum time step satisfying Courant condition)
2
2

2
1

x

y
  ( x, y )   ( x, y )   ( x, y, t )    ( x, y )   ( x, y ) ( t 
)
2
2
t
2 x  y
2
FTCS

1 x 2 y 2 
1
1
n 1
n
n
n
n
n
n
n

 jk   jk 



(


2



)

(


2



)
jk
j 1k
jk
j 1k
jk 1
jk
jk 1 
2
2 
2
2

2 x  y 
x
y


1 x 2 y 2 
1
1
n
n
n
n





(



)

(



)
jk
j 1k
j 1k
jk 1
jk 1 
2
2 
2
2

2 x  y 
x
y

Jacobi method
Simultaneous OverRelaxation (SOR)
The convergence of the Jacobi method is quite slow.
Furthermore, the larger the system, the slower the convergence.
Simultaneous OverRelaxation (SOR) :
the Jacobi method is modified in two ways,

n 1
jk
 (1   ) 
n
jk
 x 2 y 2
2 x 2  y 2


1
1
    jk  2 ( jn1k   jn11k )  2 ( jkn 1   jkn 11 ) 
x
y


1. Improved values are used as soon as they become available.
2. Relaxation parameter ω tries to overshoot for going to the final result.
( 1<ω<2)
Plasma Application
Modeling POSTECH
Simultaneous OverRelaxation (SOR)
MATLAB The Language of Technical Computing
Run: dftcs.m
MATLAB PDE
>> dftcs
dftcs - Program to solve the diffusion equation using the Forward Time Centered Space scheme.
Enter time step: 0.0001
Enter the number of grid points: 51
Solution is expected to be stable
O.V. Manuilenko
T ( x, t )
 2T ( x, t )

t
x 2
Plasma Application
Modeling Group
POSTECH
MATLAB The Language of Technical Computing
Run: dftcs.m
MATLAB PDE
>> dftcs
dftcs - Program to solve the diffusion equation using the Forward Time Centered Space scheme.
Enter time step: 0.00015
Enter the number of grid points: 61
WARNING: Solution is expected to be unstable
O.V. Manuilenko
Plasma Application
Modeling Group
POSTECH
MATLAB The Language of Technical Computing
Run: neutrn.m >> neutrn
Program to solve the neutron diffusion equation using the FTCS.
Enter time step: 0.0005
Enter the number of grid points: 61
Enter system length: 2 => System length is subcritical
Solution is expected to be stable
Enter number of time steps: 12000
O.V. Manuilenko
MATLAB PDE
n( x, t )
 2 n( x, t )

 n( x, t )
t
x 2
Plasma Application
Modeling Group
POSTECH
MATLAB The Language of Technical Computing
Run: neutrn.m >> neutrn
MATLAB PDE
Program to solve the neutron diffusion equation using the FTCS.
Enter time step: 0.0005
Enter the number of grid points: 61
Enter system length: 4 => System length is supercritical
Solution is expected to be stable
Enter number of time steps: 12000
O.V. Manuilenko
Plasma Application
Modeling Group
POSTECH
MATLAB The Language of Technical Computing
MATLAB PDE
Run: advect.m >> advect
advect - Program to solve the advection equation using the various hyperbolic PDE schemes:
FTCS, Lax, Lax-Wendorf
u ( x, t )
u ( x, t )

0
Enter number of grid points: 50
t
x
Time for wave to move one grid spacing is 0.02
Enter time step: 0.002
Wave circles system in 500 steps
Enter number of steps: 500
FTCS
O.V. Manuilenko
FTCS
Plasma Application
Modeling Group
POSTECH
MATLAB The Language of Technical Computing
MATLAB PDE
Run: advect.m >> advect
advect - Program to solve the advection equation using the various hyperbolic PDE schemes:
FTCS, Lax, Lax-Wendorf
Enter number of grid points: 50
Time for wave to move one grid spacing is 0.02
Enter time step: 0.02
Wave circles system in 50 steps
Enter number of steps: 50
Lax
O.V. Manuilenko
Lax
Plasma Application
Modeling Group
POSTECH
MATLAB The Language of Technical Computing
Run: relax.m
MATLAB PDE
>> relax
relax - Program to solve the Laplace equation using Jacobi, Gauss-Seidel and SOR methods on a
square grid
 2 ( x, y)  2 ( x, y )

0
Enter number of grid points on a side: 50
x 2
y 2
Theoretical optimum omega = 1.88184
Enter desired omega: 1.8
Potential at y=L equals 1
Potential is zero on all other boundaries
Desired fractional change = 0.0001
O.V. Manuilenko
Plasma Application
Modeling Group
POSTECH