# PDE Mathematica

Document Sample

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 ),  : thermal diffusivit )
y

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
2                    T jn1  2T jn  T jn1
T ( x j , tn ) 
x 2
x 2
t
T jn1  T jn  tS n            (T jn1  2T jn  T jn1 )
x 2
j

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 )                                               t 1

x   2
2
CTCS method for heat equation
(Both the time and space derivatives
are center-differenced.)
2 t n
T jn1  T jn1  2tS n          (T j 1  2T jn  T jn1 )
x
j         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 )       [T jn 1  1 (T jn1  T jn1 )]
              2       1         1

2t                               2t
The resulting difference equation is
n 1     n 1     n 1             2 t n
T j  2 (T j 1  T j 1 )  2tS j 
1                         n
(T j 1  2T jn  T jn1 )
x 2
Courant condition for Lax method
t    1
   ( Second-order accuracy in both time and space )
x 2
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 n  x2t (T jn1  2T jn  T jn1 )
j

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 n 1 )
j

 x 2t (T jn1  2T jn  T jn1  T jn1  2T jn 1  T jn1 )
1                    1      
( This is stable for any choice of time steps and second-order accurate in time. )

Plasma Application
Modeling POSTECH
Crank Nicolson Algorithm

Initial                 Exact solution
  2t / L2
conditions              T ( x, t )  e                  sin(x / L)

Crank-
Nicolson
scheme

Plot
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 )
 x y 
t                          
t                                     t
 T jk1  T jk  tS n 
n        n
(T jn1k  2T jk  T jn1k ) 
n
(T jk 1  2T jk  T jk 1 )
n          n      n

x 2                                    y 2
jk

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 
z ( x, y, t )  c ( x, y) 2  2  z ( x, y, t )
2
 x y 
t 2
       
z ( x, y,0)  z0 ( x, y), z ( x, y,0)  v0 ( x, y)

Initial condition :
z0  v0  0
Boundary            z (0, y, t )  z (1, y, t )  z ( x,1, t )  0,
condition :         z ( x,0, t )  x(1  x) sin 2t

Wave speed : c( x, y )  10   0.7 exp[ 2( x  0.5)  5( y  0.5) ]
1                          2             2
1

CTCS method for the wave equation :
n 1            n 1
c 2 t 2 n                               c 2 t 2 n
z jk  2 z jk  z jk          ( z j 1k  2 z jk  z j 1k )          ( z jk 1  2T jk  T jk 1 )
n             jk                    n      n           jk                    n      n

x 2                                    y 2
c 2 t 2 1
Courant condition :          ( for x  y   )
  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 0  v0 jk t  2 a jk
2
t
jk

c2                                            c2
a jk                           2z  z           )            ( z n 1  2T jk  T jk 1 )
jk           0          0       0             jk                     n      n
(z   j 1k              j 1k
x                                            y
2                   jk                        2       jk

Since initial velocity and value,
v0 jk  z 0  0
jk

 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              jk 1  2 jk   jk 1
                                 jk
x     2
y   2

Direct Solution for Poisson’s equation
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 :
I.  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.
(Maximum time step satisfying Courant condition)
                                               1 x 2 y 2
  ( x, y)   ( x, y)   ( x, y, t )    ( x, y)   ( x, y) ( t 
2                                        2
)
t                                              2 x  y
2       2

FTCS
1 x 2 y 2             1                                1                              
 jk   jk 
n 1   n
   jk  2 ( jn1k  2 jk   jn1k )  2 ( jk 1  2 jk   jk 1 ) 
n                     n          n      n

2 x 2  y 2 
         x                               y                              

1 x 2 y 2             1                        1                      
                  jk  2 ( jn1k   jn1k )  2 ( jk 1   jk 1 ) 
n         n

2 x 2  y 2 
         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,
 x 2 y 2
   n 1
 (1   ) 
n

2 x 2  y 2
jk                 jk

          1                        1              n 1 
    jk  2 ( j 1k   j 1k )  2 ( jk 1   jk 1 ) 

n         n 1           n

         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                                  MATLAB PDE

Run: dftcs.m           >> dftcs
dftcs - Program to solve the diffusion equation using the Forward Time Centered Space scheme.
Enter time step: 0.0001                      T ( x, t )     2T ( x, t )

Enter the number of grid points: 51             t             x 2
Solution is expected to be stable

Plasma Application
Modeling Group
O.V. Manuilenko                                                                  POSTECH
MATLAB The Language of Technical Computing                             MATLAB PDE

Run: dftcs.m          >> 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

Plasma Application
Modeling Group
O.V. Manuilenko                                                                POSTECH
MATLAB The Language of Technical Computing                               MATLAB PDE

Run: neutrn.m >> neutrn        Program to solve the neutron diffusion equation using the FTCS.
Enter time step: 0.0005
n( x, t )     2 n ( x, t )
Enter the number of grid points: 61                                                           n ( x, t )
Enter system length: 2 => System length is subcritical             t              x 2
Solution is expected to be stable
Enter number of time steps: 12000

Plasma Application
Modeling Group
O.V. Manuilenko                                                                       POSTECH
MATLAB The Language of Technical Computing                            MATLAB PDE

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: 4 => System length is supercritical
Solution is expected to be stable
Enter number of time steps: 12000

Plasma Application
Modeling Group
O.V. Manuilenko                                                                POSTECH
MATLAB The Language of Technical Computing                             MATLAB PDE

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                                              FTCS

Plasma Application
Modeling Group
O.V. Manuilenko                                                                POSTECH
MATLAB The Language of Technical Computing                             MATLAB PDE

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                                             Lax

Plasma Application
Modeling Group
O.V. Manuilenko                                                                POSTECH
MATLAB The Language of Technical Computing                            MATLAB PDE

Run: relax.m      >> 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

Plasma Application
Modeling Group
O.V. Manuilenko                                                               POSTECH

DOCUMENT INFO
Shared By:
Categories:
Tags:
Stats:
 views: 461 posted: 2/9/2012 language: English pages: 28
How are you planning on using Docstoc?