ECE695: Special Topics in EM-Plasma Simulations JK LEE (Spring, 2005)

1
ECE695 Special Topics in EM-Plasma
SimulationsJK LEE (Spring, 2005)

• ODE Solvers
• PIC-MCC
• PDE Solvers (FEM and FDM)
• Linear NL Eq. Solvers

2


??- ??   2005-2
????( ??-??-??)    ?? ??? ????? (3-0-3) ????    ???
????-??(?????)     - 01(EECE695) ????-????    ???? - ??
????/???     ?, ? 1615 1730 LG???-0106   ?, ? 1615 1730 LG???-0106   ?, ? 1615 1730 LG???-0106
???    279-2083  279-2083  279-2083























Computational methods for advanced EM including plasmas and charged particles to be described with hands-on experience projects




-ECE261 EM (II) -ECE490D or others



-Mid-Term Quiz -Homeworks -Projects ECE490 Special Topics on Simulation (Spring 2006)




-(A) CK Birdsall and AB Langdon Plasma Phys. via Computer Simul. (Adam Hilger 1991) -(B) S. Nakamura Applied Numerical Methods in FORTRAN (Prentice 1991) -(C) JD Hoffman Numerical Methods for Engineers and Scientists (McGraw Hill 1992)










Week 1-3 ODE solvers with projects in Fortran Week 4-6 ES1 and PIC-MCC in 1d 2d with projects Week 7-9 FEM Week 10-12 PDE solvers via FDM Week 12-15 Nonlinear eq. Solvers and others




ECE490 (Spring 2006) for entry level simulations recommended

3





Week 1-3 ODE solvers with projects in Fortran Week 4-6 ES1 and PIC-MCC in 1d 2d with projects Week 7-9 FEM Week 10-12 PDE solvers via FDM Week 12-15 Nonlinear eq. Solvers and others

4
Programs of Initial Value Problem Shooting
Method for BVP ODE
Sung Jin Kim and Jae Koo Lee
5
(No Transcript)
6
(No Transcript)
7
Initial Value Problems of Ordinary Differential
Equations
Program 9-1 Second-order Runge-kutta Method
, y(0)1, y(0)0
Changing 2nd order ODE to 1st order ODE,
y(0)1
(1)
z(0)0
(2)
8
Second-Order Runge-Kutta Method
By 2nd order RK Method
9
Program 9-1
k1 hz l1 -h(bmz kmy) k2 h(z
l1) l2 -h(bm(z l1) km(y k1)) y y
(k1 k2)/2 z z (l1 l2)/2 printf(
" 12.6f 12.5e 12.5e \n", time, y, z )
exit(0)
/ CSL/c9-1.c Second Order Runge-Kutta Scheme
(Solving the problem II of Example 9.6) /
include ltstdio.hgt include ltstdlib.hgt include
ltmath.hgt / time t y,z y,y' kount
number of steps between two lines of printing
k, m, b k, M(mass), B(damping coefficient) in
Example 9.6 int main() int kount, n,
kstep0 float bm, k1, k2, km, l1, l2 static
float time, k 100.0, m 0.5, b 10.0, z
0.0 static float y 1.0, h 0.001 printf(
"CSL/C9-1 Second Order Runge-Kutta Scheme \n" )
printf( " t y z\n" ) printf( " 12.6f 12.5e
12.5e \n", time, y, z ) km k/m bm b/m
for( n 1 n lt 20 n ) for( kount 1
kount lt 50 kount ) kstepkstep1 time
hkstep
2nd order RK
result
CSL/C9-1 Second Order Runge-Kutta Scheme
t y z 0.000000
1.00000e00 0.00000e00 0.100000
5.08312e-01 -6.19085e00 0.200000
6.67480e-02 -2.46111e00 0.300000
-4.22529e-02 -1.40603e-01 0.400000
-2.58300e-02 2.77157e-01 0.500000
-4.55050e-03 1.29208e-01 0.600000
1.68646e-03 1.38587e-02 0.700000
1.28624e-03 -1.19758e-02 0.800000
2.83107e-04 -6.63630e-03 0.900000
-6.15151e-05 -1.01755e-03 1.000000
-6.27664e-05 4.93549e-04
10
Various Numerical Methods
h0.1
Exact Solution
h0.01
h0.001
11
Error Estimation
Error estimation
Fourth order Runge-Kutta
12
Program 9-2 Fourth-order Runge-Kutta Scheme
A first order Ordinary differential equation
y(0)1
Fourth-order RK Method
13
Boundary Value Problems of Ordinary Differential
Equations Scharfetter-Gummel method
Sung Soo Yang and Jae Koo Lee
Homepage http//jkl.postech.ac.kr
14
Boundary value problems
Type of Problems Advantages Disadvantages
Shooting method An existing program for initial value problems may be used Trial-and-error basis. Application is limited to a narrow class of problems. Solution may become unstable.
Finite difference method using the tridiagonal solution No instability problem. No trial and error. Applicable to nonlinear problems with iteration. Problem may have to be developed for each particular problem.
15
Program 10-1
Solve difference equation,
With the boundary conditions,
x 0
1
2
9
10
i 0
1
2
9
10
known y(0)1
Especially for i 1,
16
Program 10-1
For i 10,
Summarizing the difference equations obtained, we
write
Tridiagonal matrix
17
Solution Algorithm for Tridiagonal Equations (1)
R2
Based on Gauss elimination
R3
18
Solution Algorithm for Tridiagonal Equations (2)
19
Solution Algorithm for Tridiagonal Equations (3)
void trdg(float a, float b, float c, float
d, int n) / Tridiagonal solution /
int i float r for ( i 2 i lt n
i ) r ai/bi - 1
bi bi - rci - 1 di di -
rdi - 1 dn dn/bn
for ( i n - 1 i gt 1 i-- )
di (di - cidi 1)/bi
return
Recurrently calculate the equations in increasing
order of i until iN is reached
Calculate the solution for the last unknown by
Calculate the following equation in decreasing
order of i
20
Program 10-1
/ CSL/c10-1.c Linear Boundary Value Problem /
include ltstdio.hgt include ltstdlib.hgt
include ltmath.hgt / ai, bi, ci, di
a(i), b(i), c(i), and d(i) n number of grid
points / int main() int i, n
float a20, b20, c20, d20, x void
trdg(float a, float b, float c, float d,
int n) / Tridiagonal solution / printf(
"\n\nCSL/C10-1 Linear Boundary Value Problem\n"
) n 10 / n Number of grid points /
for( i 1 i lt n i ) x
i ai -2 bi 5 ci -2
di exp( -0.2x )
21
Program 10-1
d1 d1 2 dn dn0.5
bn 4.5 trdg( a, b, c, d, n ) d0
1 / Setting d0 for printing purpose /
printf( "\n Grid point no. Solution\n" )
for ( i 0 i lt n i ) printf(
" 3.1d 12.5e \n", i, di )
exit(0)
CSL/C10-1 Linear Boundary Value Problem Grid
point no. Solution 0 1.00000e00
1 8.46449e-01 2 7.06756e-01
3 5.85282e-01 4 4.82043e-01
5 3.95162e-01 6 3.21921e-01
7 2.59044e-01 8 2.02390e-01
9 1.45983e-01 10 7.99187e-02
22
Program 10-3
An eigenvalue problem of ordinary differential
equation
We assume
23
Scharfetter-Gummel method
Tridiagonal matrix

• 2D discretized continuity eqn. integrated by the
  alternative direction implicit (ADI) method

Scharfetter-Gummel method
24
Mullers method for solving non-linear equations
FORTRAN features
N. Babaeva
25
Muller_4.c
R0.01
The beam is stronger
Real and imaginary roots of the dispersion
relation
R0.01
26
W.H.Press et al Numerical recipes in C The art
of Scientific Computing
Mullers method
Joe D. Hoffman Numerical Methods for Engineers
and scientists McGraw-Hill, Inc. 1993
In Mullers method, three approximations to the
zero are required. The next approximation is the
zero of the parabola that passes through the
three points. Both real zeros and complex zeros
can be found.Mullers method generalizes the
secant method, but uses quadratic interpolation
among three points instead of linear
interpolation between two. Solving for the zeros
of the quadratic allows the method to find
complex pairs of roots. Given three previous
guesses for the root
And the values of polynomial P(x) at those
points, the next approximation is produces by
the following formulas
Figure, illustrating the Mullers method
Mullers method is also used for finding complex
zeros of analytic functions (not just
polynomials) in the complex plane.
27
Plasma Display Panel
Plasma Display Panel
Many Pixels
the flat panel display using phosphor
luminescence by UV photons produced in plasma gas
discharge
PDP structure
28
Simulation domain
Sustain 1
Sustain 2
dielectric layer
dielectric and phosphor layer
address
Finite-Difference Method
?
?
?
j
? Electric field, Density
?
?
?
? Potential, Charge
?
?
?
? Flux of x and y
j1
i
i1
? Light, Luminance, Efficiency, Power, Current
and so on
29
Flow chart
fl2p.c
initial.c
pulse.c
charge.c
field.c
flux.c
continuity.c
source.c
history.c
diagnostics.c
Calculate efficiency
time_step.c
current.c, radiation.c, dump.c, gaspar.c,
mu_n_D.c, gummel.c
30
Basic equations

• Continuity Equation with Drift-Diffusion Approx.

and

• Poissons Equation

? surface charge density on the dielectric
surfaces

• Boundary conditions on dielectric surface

for ion
for electron
Mobility-driven drift term
Isotropic thermal flux term
for secondary electron
for excited species
31
Partial Differential Eqs.
General form of linear second-order PDEs with
two independent variables
In case of elliptic PDEs,
?
Jacobi-Iteration method
?
Gauss-Seidel method
?
Successive over-relaxation (SOR) method
In case of parabolic PDEs,
?
Alternating direction implicit (ADI) method
32
Continuity equation (1)
density nsp
Alternating direction implicit (ADI) method
ADI method uses two time steps in two dimension
to update the quantities between t and t?t.
During first ?t/2, the integration sweeps along
one direction (x direction) and the other
direction (y direction) is fixed. The temporary
quantities are updated at t?t/2. With these
updated quantities, ADI method integrates the
continuity equation along y direction with fixed
x direction between t?t/2 and t?t.
1st step
(k means the value at time t) ( means the
temporal value at time t?t/2 )
Discretized flux can be obtained by
Sharfetter-Gummel method.
Spatially discretized forms are converted to
tridiagonal systems of equations which can be
easily solved.
33
Tridiagonal matrix (1)
Based on Gauss elimination
R2
R3
34
Tridiagonal matrix (2)
35
Tridiagonal matrix (3)
/ Tridiagonal solution / void trdg(float a,
float b, float c, float d, int n)
int i float r for ( i 2 i lt n
i ) r ai/bi - 1
bi bi - rci - 1 di di -
rdi - 1 dn dn/bn
for ( i n - 1 i gt 1 i-- )
di (di - cidi 1)/bi
return
Calculate the equations in increasing order of i
until iN is reached.
Ri
Calculate the solution for the last unknown by
Calculate the following equation in decreasing
order of i
36
Continuity equation (2)
2nd step
From the temporally updated density calculated in
the 1st step, we can calculated flux in
x-direction (?) at time t?t/2. Using these
values, we calculate final updated density with
integration of continuity equation in y-direction.
(k1 means the final value at time t?t) (
means the temporal value at time t?t/2 )
(tridiagonal matrix)
From the final updated density calculated in the
2nd step, we can calculated flux in y-direction
(?k1) at time t.
37
Poissons eq. (1)
Poisson equation is solved with a successive over
relation (SOR) method. The electric field is
taken at time t when the continuity equations are
integrated between t and t?t. Time is integrated
by semi-implicit method in our code. The electric
field in the integration of the continuity
equation between t and t?t is not the field at
time t, but rather a prediction of the electric
field at time t?t. The semi-implicit integration
of Poisson equation is followed as
Density correction by electric field change
between t and t?t
The continuity eq. and flux are coupled with
Poissions eq.
This Poissons eq can be discriminated to x and y
directions, and written in matrix form using the
five-point formula in two dimensions.
38
Poissons eq. (2)
j1
ci, j
ai, j
bi, j
j
is the surface charge density
accumulating on intersection between plasma
region and dielectric.
di, j
j-1
i-1
i
i1
Solved using SOR method
39
Finite Difference Method (II)
Grid 100x50
40
Full rectangle
100V
-100V
100V
-100V
41
Full circle
100V
-100V
100V
-100V
42
What is FEMLAB?

• FEMLAB a powerful interactive environment for
  modeling and
• solving various kinds of scientific and
  engineering problems based
• on partial differential equations (PDEs).

• Overview

• Finite element method
• GUI based on Java
• Unique environments for modeling
• (CAD, Physics, Mesh, Solver, Postprocessing)
• Modeling based on equations (broad application)
• Predefined equations and User-defined
  equations
• No limitation in Multiphysics
• MATLAB interface (Simulink)

• Mathematical application modes and types of
  analysis

• Mathematical application modes
• 1. Coefficient form suitable for linear or
  nearly linear models.
• 2. General form suitable for nonlinear models
• 3. Weak form suitable for models with PDEs on
  boundaries, edges,
• and points, or for models using terms with
  mixed space and time
• derivatives.
• Various types of analysis
• 1. Eigenfrequency and modal analysis
• 2. Stationary and time-dependent analysis
• 3. Linear and nonlinear analysis

Reference Manual of FEMLAB Software
43
Useful Modules in FEMLAB

• Application areas

• Microwave engineering
• Optics
• Photonics
• Porous media flow
• Quantum mechanics
• Radio-frequency components
• Semiconductor devices
• Structural mechanics
• Transport phenomena
• Wave propagation

• Acoustics
• Bioscience
• Chemical reactions
• Diffusion
• Electromagnetics
• Fluid dynamics
• Fuel cells and electrochemistry
• Geophysics
• Heat transfer
• MEMS

• Additional Modules

1. Application of Chemical engineering Module

• Momentum balances
• - Incompressible Navier-Stokes eqs.
• - Darys law
• - Brinkman eqs.
• - Non-Newtonian flow

• Mass balances
• - Diffusion
• - Convection and Conduction
• - Electrokinetic flow
• - Maxwell-stefan diffusion and convection

• Energy balances
• - Heat equation
• - Heat convection and conduction

2. Application of Electromagnetics Module

• - Electrostatics
• - Conductive media DC
• Magnetostatic
• Low-frequency electromagnetics
• - In-plane wave propagation
• Axisymmetric wave propagation
• Full 3D vector wave propagation
• Full vector mode analysis in 2D and 3D

3. Application of the Structural Mechanics Module

• Plane stress
• Plane strain
• 2D, 3D beams, Euler theory
• Shells

44
FEMLAB Environment
Model Navigator
Pre-defined Equations
45
Monoconical RF Antenna

• Introduction of Conical RF Antenna

• Antennas are used to couple guided
  electromagnetic energy to waves
• in free space.
• Antenna engineering design involves minimizing
  reflection losses and
• obtaining desired directional properties for a
  wide range of frequencies.
• - antenna impedance, radiation pattern,
  frequency analysis.
• Modeling can be extremely useful as it reduces
  the need for building
• prototypes and performing measurements.

• Antenna parameters

Antenna
metal
teflon ?r2.07
Ground plane
Coaxial feed
Coaxial feed rinner1.5mm router4.916mm Z50
46
PIC Overview

• PIC Codes Overview

• PIC codes simulate plasma behavior of a large
  number of charges particles using a few
  representative super particles.
• These type of codes solve the Newton-Lorentz
  equation of motion to move particles in
  conjunction with Maxwells equations (or a
  subset).
• Boundary conditions are applied to the particles
  and the fields to solve the set of equations.
• PIC codes are quite successful in simulating
  kinetic and nonlinear plasma phenomenon like ECR,
  stochastic heating, etc.

47
Capacitively Coupled Plasma 1D PIC-MCC

j 1, ? , N


Sheath


Bulk Plasma
Sheath
Substrate

• MCC (Monte-Carlo Collision) Processes
• - Electron-Neutral Collisions
• (Ionization, Scattering, Excitation)
• - Ion-Neutral Collisions
• (Charge-exchange, Scattering)

• 1D Asymmetric Dual-Freq. Voltage-Driven System

48
PIC-MCC Flow Chart

• Particles in continuum space
• Fields at discrete mesh locations in space
• Coupling between particles and fields

I
II
V
IV
III
IV
Fig Flow chart for an explicit PIC-MCC scheme
49
I. Particle Equations of Motion

• Newton-Lorentz equations of motion

• In finite difference form, the leapfrog method

Fig Schematic leapfrog integration
50
III. Electrostatic Field Model

• Possions equation

• Finite difference form in 1D planar geometry

• Boundary condition External circuit

Fig Schematic one-dimensional bounded
plasma with external circuit
51
PDP Structure
AC PDP
Discharge in PDP
52
Striation Profiles in PDP 2D PIC/MCC
Anode
Cathode

• Pressure dependence of striations
• Number of peaks depend on the pressure and
  electrode size.

53
Shooting Method for Boundary Value Problem ODEs
Rewrite the second-order ODE as two first-order
ODEs
We should assume a initial value of z(a). z(a) is
determined for which y(b)B by secant method.
54
Errors in the Solution by the Shooting Method
55
C   CSL/F9-1.FOR    SECOND ORDER RUNGE-KUTTA
SCHEME C                   (SOLVING THE PROBLEM
II OF EXAMPL 9.6)       REAL M,K,K1,K2,L1,L2,KM
      PRINT ,'CSL/F9-1   SECOND ORDER
RUNGE-KUTTA SCHEME'       DATA T,   K,    M,  
B,    Z,   Y,   H          /0.0,100.0, 0.5,
10.0, 0.0, 1.0, 0.001/       PRINT ,'    
T           Y             Z'       PRINT 1,T,Y,Z
1     FORMAT( F10.5,  1P2E13.5)       KMK/M
      BMB/M
56
DO 20 N1,20         DO 10 KOUNT1,50
            TTH             K1HZ
            L1-H(BMZ KMY)
            K2H(ZL1)             L2-H(BM(Z
L1) KM(YK1))             YY(K1K2)/2
            ZZ(L1L2)/2 10      CONTINUE
          PRINT 1,T,Y,Z 20    CONTINUE
      STOP       END
57
C-----CSL/F9-2.FOR   FOURTH-ORDER RUNGE-KUTTA
SCHEME (FORTRAN)       REAL K1,K2,K3,K4
      PRINT       PRINT,'CSL/F9-2.FOR 
FOURTH-ORDER RUNGE KUTTA SCHEME (FORTRAN)' 1    
PRINT       PRINT , 'INTERVAL OF T FOR
PRINTING ?'          READ , XPR       PRINT ,
'NUMBER OF STEPS IN ONE PRINTING INTERVAL ?'
         READ , I       PRINT ,'MAXIMUM T ?'
         READ , XL C  Setting the initial
value of the solution       Y1
58
 HXPR/I       PRINT , 'H', H       XP0
      HHH/2       PRINT       PRINT
,'--------------------  PRINT ,'      
T                Y'       PRINT
,'----------------------       PRINT 82, XP,Y
82    FORMAT( 1X,F10.6, 7X,1PE15.6) 30    DO 40
J1,I          XBXP          XPXPH
         YNY          XMXBHH
         K1HFUN(YN,XB)          K2HFUN(YNK1
/2,XM)          K3HFUN(YNK2/2,XM)
         K4HFUN(YNK3,XP)          YYN
(K1K22K32K4)/6 40    CONTINUE
59
   IF (XP.LE.XL) GO TO 30       PRINT
,'------------------------       PRINT
      PRINT ,' MAXIMUM X LIMIT IS EXCEEDED'
      PRINT 200   PRINT       PRINT,'TYPE
1 TO CONTINUE, 0 TO STOP.'       READ ,K
      IF(K.EQ.1) GOTO 1       PRINT       END
      FUNCTION FUN(Y,X)       FUN -Y/(XXYY
)       RETURN       END
60
Programs of Initial Value Problem Shooting
Method for BVP ODE
Sung Jin Kim and Jae Koo Lee
61
(No Transcript)
62
(No Transcript)
63
Initial Value Problems of Ordinary Differential
Equations
Program 9-1 Second-order Runge-kutta Method
, y(0)1, y(0)0
Changing 2nd order ODE to 1st order ODE,
y(0)1
(1)
z(0)0
(2)
64
Second-Order Runge-Kutta Method
By 2nd order RK Method
65
Program 9-1
k1 hz l1 -h(bmz kmy) k2 h(z
l1) l2 -h(bm(z l1) km(y k1)) y y
(k1 k2)/2 z z (l1 l2)/2 printf(
" 12.6f 12.5e 12.5e \n", time, y, z )
exit(0)
/ CSL/c9-1.c Second Order Runge-Kutta Scheme
(Solving the problem II of Example 9.6) /
include ltstdio.hgt include ltstdlib.hgt include
ltmath.hgt / time t y,z y,y' kount
number of steps between two lines of printing
k, m, b k, M(mass), B(damping coefficient) in
Example 9.6 int main() int kount, n,
kstep0 float bm, k1, k2, km, l1, l2 static
float time, k 100.0, m 0.5, b 10.0, z
0.0 static float y 1.0, h 0.001 printf(
"CSL/C9-1 Second Order Runge-Kutta Scheme \n" )
printf( " t y z\n" ) printf( " 12.6f 12.5e
12.5e \n", time, y, z ) km k/m bm b/m
for( n 1 n lt 20 n ) for( kount 1
kount lt 50 kount ) kstepkstep1 time
hkstep
2nd order RK
result
CSL/C9-1 Second Order Runge-Kutta Scheme
t y z 0.000000
1.00000e00 0.00000e00 0.100000
5.08312e-01 -6.19085e00 0.200000
6.67480e-02 -2.46111e00 0.300000
-4.22529e-02 -1.40603e-01 0.400000
-2.58300e-02 2.77157e-01 0.500000
-4.55050e-03 1.29208e-01 0.600000
1.68646e-03 1.38587e-02 0.700000
1.28624e-03 -1.19758e-02 0.800000
2.83107e-04 -6.63630e-03 0.900000
-6.15151e-05 -1.01755e-03 1.000000
-6.27664e-05 4.93549e-04
66
Various Numerical Methods
h0.1
Exact Solution
h0.01
h0.001
67
Error Estimation
Error estimation
Fourth order Runge-Kutta
68
Program 9-2 Fourth-order Runge-Kutta Scheme
A first order Ordinary differential equation
y(0)1
Fourth-order RK Method
69
Program 9-2 Fourth-order Runge-Kutta Scheme
CSL/C9-2 Fourth-Order Runge-Kutta Scheme
Interval of t for printing ? 1 Number of steps
in one printing interval? 10 Maximum t? 5 h0.1
-------------------------------------- t
y ----------------------------------
---- 0.00000 0.000000e00
1.00000 1.410686e00 2.00000
8.839368e00 3.00000
1.125059e02 4.00000 3.734231e03
5.00000 3.357971e05
6.00000 8.194354e07 --------------------
------------------ Maximum t limit exceeded
Main algorithm for 4th order RK method at program
9-2
do for( j 1 j lt nstep_pr j ) t_old
t_new t_new t_new h yn y t_mid
t_old hh yn y k1 hfun(
yn, t_old ) ya yn k1/2 k2 hfun( ya,
t_mid ) ya yn k2/2 k3 hfun( ya, t_mid
) ya yn k3 k4 hfun( ya, t_new )
y yn (k1 k22 k32 k4)/6 double
fun(float y, float t) float fun_v fun_v
ty 1 / Definition of f(y,t) / return( fun_v
)
4th order RK
Interval of t for printing ? 1 Number of steps in
one printing interval? 100 Maximum t? 5 h0.01
-------------------------------------- t
y -----------------------------------
--- 0.00000 0.000000e00
1.00000 1.410686e00 2.00000
8.839429e00 3.00000
1.125148e02 4.00000 3.735819e03
5.00002 3.363115e05
--------------------------------------
70
Program 9-3 4th order RK Method for a Set of ODEs
A second order Ordinary differential equation
y1(0)1,
y1(0)1, y1(0)0
The 4th order RK method for the set of two
equations(Nakamuras book p332)
do for( n 1 n lt ns n ) t_old t_new
/ Old time / t_new t_new h / New time /
t_mid t_old hh / Midpoint time / for( i
1 i lt No_of_eqs i ) yai yi f(
k1, ya, t_old, h ) for( i 1 i lt
No_of_eqs i ) yai yi k1i/2 f(
k2, ya, t_mid, h ) for( i 1 i lt
No_of_eqs i ) yai yi k2i/2 f(
k3, ya, t_mid, h ) for( i 1 i lt
No_of_eqs i ) yai yi k3i f(
k4, ya, t_new, h ) for( i 1 i lt
No_of_eqs i ) yi yi (k1i k2i2
k3i2 k4i)/6
void f(float k, float y, float t, float h)
k1 y2h k2 -y1h / More
equations come here if the number of equations
are greater./ return
71
Results of Program 9-3
CSL/C9-3 Fourth-Order Runge-Kutta Scheme
for a Set of Equations Interval of
t for printing ? 1 Number of steps in one print
interval ? 10 Maximum t to stop calculations ?
5.0 h 0.1 t y(1),
y(2), ..... 0.0000 1.00000e00
0.00000e00 1.0000 5.40303e-01
-8.41470e-01 2.0000 -4.16145e-01
-9.09298e-01 3.0000 -9.89992e-01
-1.41123e-01 4.0000 -6.53646e-01
7.56800e-01 5.0000 2.83658e-01
9.58925e-01 6.0000 9.60168e-01
2.79420e-01 Type 1 to continue, or 0 to stop. 1
Interval of t for printing ? 1 Number of steps
in one print interval ? 1 Maximum t to stop
calculations ? 5.0 h 1 t y(1),
y(2), ..... 0.0000 1.00000e00
0.00000e00 1.0000 5.41667e-01
-8.33333e-01 2.0000 -4.01042e-01
-9.02778e-01 3.0000 -9.69546e-01
-1.54803e-01 4.0000 -6.54173e-01
7.24103e-01 5.0000 2.49075e-01
9.37367e-01 Type 1 to continue, or 0 to stop.
Interval of t for printing ? 1 Number of steps
in one print interval ? 100 Maximum t to stop
calculations ? 5.0 h 0.01 t
y(1), y(2), ..... 0.0000
1.00000e00 0.00000e00 1.0000
5.40302e-01 -8.41471e-01 2.0000
-4.16147e-01 -9.09298e-01 3.0000
-9.89992e-01 -1.41120e-01 4.0000
-6.53644e-01 7.56802e-01 5.0000
2.83662e-01 9.58924e-01 Type 1 to continue, or
0 to stop. 1
72
Shooting Method for Boundary Value Problem ODEs
Definition a time stepping algorithm along with
a root finding method for choosing
the appropriate initial conditions which
solve the boundary value problem.
Second-order Boundary-Value Problem
y(a)A and y(b)B
73
Computational Algorithm of Shooting Method
Computational Algorithm
1. Solve the differential equation using a
time-stepping scheme with the initial
conditions y(a)A and y(a)A.
2. Evaluate the solution y(b) at xb and compare
this value with the target value of y(b)B.
3. Adjust the value of A (either bigger or
smaller) until a desired level of tolerance
and accuracy is achieved. A bisection or secant
method for determining values of A, for
instance, may be appropriate.
4. Once the specified accuracy has been achieved,
the numerical solution is complete and is
accurate to the level of the tolerance chosen and
the discretization scheme used in the
time-stepping.
74
Shooting Method for Boundary Value Problem ODEs
Rewrite the second-order ODE as two first-order
ODEs
We should assume a initial value of z(a). z(a) is
determined for which y(b)B by secant method.
75
Example of Shooting Method
T(0)0 and T(1.0)100
Ex)
Sol) Rewrite the second-order ODE as two
first-order ODEs
T(0)0.0
S(0)T(0)
By modified Euler method
Let ?x0.25, S(0.0)(1)50, S(0.0)(2)100
76
Solution by the Shooting Method
77
Errors in the Solution by the Shooting Method
78
Equilibrium Method
y(a)A and y(b)B
Ex)
T(0)0 and T(1.0)100
Let ?x0.25, Ta0, ?2.0
x0.25
0 0 -100
-2.25 1.0 0 1.0 -2.25 1.0 0
1.0 -2.25
T1 T2 T3

x0.50
x0.75
79
Executing Programs for 490D class
An user ID will be made as team name at pam10
computer.
A password is the same as an user ID
You can execute your programs to connect your PC
to pam10 computer by x-manager or telnet.

• A C program is compiled by a gcc compiler at the
  linux PC.
• gcc o a a.c lm
• where, a is a changeable execution sentence.

80
Boundary Value Problems of Ordinary Differential
Equations Scharfetter-Gummel method
Sung Soo Yang and Jae Koo Lee
Homepage http//jkl.postech.ac.kr
81
Boundary value problems
Type of Problems Advantages Disadvantages
Shooting method An existing program for initial value problems may be used Trial-and-error basis. Application is limited to a narrow class of problems. Solution may become unstable.
Finite difference method using the tridiagonal solution No instability problem. No trial and error. Applicable to nonlinear problems with iteration. Problem may have to be developed for each particular problem.
82
Program 10-1
Solve difference equation,
With the boundary conditions,
x 0
1
2
9
10
i 0
1
2
9
10
known y(0)1
Especially for i 1,
83
Program 10-1
For i 10,
Summarizing the difference equations obtained, we
write
Tridiagonal matrix
84
Solution Algorithm for Tridiagonal Equations (1)
R2
Based on Gauss elimination
R3
85
Solution Algorithm for Tridiagonal Equations (2)
86
Solution Algorithm for Tridiagonal Equations (3)
void trdg(float a, float b, float c, float
d, int n) / Tridiagonal solution /
int i float r for ( i 2 i lt n
i ) r ai/bi - 1
bi bi - rci - 1 di di -
rdi - 1 dn dn/bn
for ( i n - 1 i gt 1 i-- )
di (di - cidi 1)/bi
return
Recurrently calculate the equations in increasing
order of i until iN is reached
Calculate the solution for the last unknown by
Calculate the following equation in decreasing
order of i
87
Program 10-1
/ CSL/c10-1.c Linear Boundary Value Problem /
include ltstdio.hgt include ltstdlib.hgt
include ltmath.hgt / ai, bi, ci, di
a(i), b(i), c(i), and d(i) n number of grid
points / int main() int i, n
float a20, b20, c20, d20, x void
trdg(float a, float b, float c, float d,
int n) / Tridiagonal solution / printf(
"\n\nCSL/C10-1 Linear Boundary Value Problem\n"
) n 10 / n Number of grid points /
for( i 1 i lt n i ) x
i ai -2 bi 5 ci -2
di exp( -0.2x )
88
Program 10-1
d1 d1 2 dn dn0.5
bn 4.5 trdg( a, b, c, d, n ) d0
1 / Setting d0 for printing purpose /
printf( "\n Grid point no. Solution\n" )
for ( i 0 i lt n i ) printf(
" 3.1d 12.5e \n", i, di )
exit(0)
CSL/C10-1 Linear Boundary Value Problem Grid
point no. Solution 0 1.00000e00
1 8.46449e-01 2 7.06756e-01
3 5.85282e-01 4 4.82043e-01
5 3.95162e-01 6 3.21921e-01
7 2.59044e-01 8 2.02390e-01
9 1.45983e-01 10 7.99187e-02
89
Program 10-3
An eigenvalue problem of ordinary differential
equation
We assume
90
Program 10-3 (Inverse Power method)
Step1
for all i and
are set to an arbitrary initial guess.
is solved by the tridiagonal solution for
Step2
Step3 The next estimate for is calculated
by the equation
is solved by the tridiagonal solution for
Step4
Step5 The operations similar to Step 3 and 4
are repeated as the iteration cycle
t increases.
Step6 The iteration is stopped when the
convergence test is satisfied.
Criterion for convergence
91
Program 10-3
while (TRUE) k k1 for
(i1 iltn i) fbi fiei for
(i1 iltn i) ai asi
bi bsi ci csi di
dsifbi trdg(a, b, c, d,
n) sb 0 s 0 for (i1 i
lt n i) fi di
s s fifi sb sb
fifbi eb ei ei
sb/s printf("3.1d 12.6e \n", k,
ei) if (fabs(1.0-ei/eb) lt ep) break
if (k gtit) printf("Iteration
limit exceeded.\n") break

main() while (TRUE) k 0
n 10 ei 1.5 it 30 ep
0.0001 for (i1 i lt n i)
asi -1.0 bsi 2.0 csi
-1.0 fi 1.0 dsi 1.0
printf("\n It. No. Eigenvalue\n")
Step 2,4
Step 1
Step 3
Step 6
92
Program 10-3
It. No. Eigenvalue 1 8.196722e-02 2
8.102574e-02 3 8.101422e-02 4
8.101406e-02 Eigenvalue 0.0810141
Eigenfunction i f(i) 1 2.84692e-01 2
5.46290e-01 3 7.63595e-01 4
9.19018e-01 5 1.00000e00 6
1.00000e00 7 9.19018e-01 8
7.63594e-01 9 5.46290e-01 10
2.84692e-01 ------------------------------------
Type 1 to continue, or0 to stop.
z 0 for (i1 i lt n i) if
(fabs(z) lt fabs(fi)) z fi for (i1 i
lt n i) fi fi/z eigen ei
printf("Eigenvalue g \n", eigen)
printf("\n Eigenfunction\n") printf("i
f(i)\n") for (i1 iltn i)
printf("3.1d 12.5e \n", i, fi)
printf("------------------------------------\n")
printf("Type 1 to continue, or0 to stop.
\n") scanf("d", kstop) if (kstop !
1) exit (0)
93
Scharfetter-Gummel method
Tridiagonal matrix

• 2D discretized continuity eqn. integrated by the
  alternative direction implicit (ADI) method

Scharfetter-Gummel method
94
Scharfetter-Gummel method
Flux is constant between half grid points and
calculated at the grid point
?
?
Flux is a linear combination of and
Analytic integration between i and i1 leads to
95
Scharfetter-Gummel method
where,
The main advantage of SG scheme is that is
provides numerically stable estimates of the
particle flux under all conditions
(Big potential difference)
(small potential difference)
96
Scharfetter-Gummel method
In case of
Drift flux
97
Scharfetter-Gummel method
In case of
Diffusion flux