1 / 97

ECE695 Special Topics in EM-Plasma

SimulationsJK LEE (Spring, 2005)

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

??- ?? 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

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

Programs of Initial Value Problem Shooting

Method for BVP ODE

Sung Jin Kim and Jae Koo Lee

(No Transcript)

(No Transcript)

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)

Second-Order Runge-Kutta Method

By 2nd order RK Method

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

Various Numerical Methods

h0.1

Exact Solution

h0.01

h0.001

Error Estimation

Error estimation

Fourth order Runge-Kutta

Program 9-2 Fourth-order Runge-Kutta Scheme

A first order Ordinary differential equation

y(0)1

Fourth-order RK Method

Boundary Value Problems of Ordinary Differential

Equations Scharfetter-Gummel method

Sung Soo Yang and Jae Koo Lee

Homepage http//jkl.postech.ac.kr

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.

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,

Program 10-1

For i 10,

Summarizing the difference equations obtained, we

write

Tridiagonal matrix

Solution Algorithm for Tridiagonal Equations (1)

R2

Based on Gauss elimination

R3

Solution Algorithm for Tridiagonal Equations (2)

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

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 )

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

Program 10-3

An eigenvalue problem of ordinary differential

equation

We assume

Scharfetter-Gummel method

Tridiagonal matrix

- 2D discretized continuity eqn. integrated by the

alternative direction implicit (ADI) method

Scharfetter-Gummel method

Mullers method for solving non-linear equations

FORTRAN features

N. Babaeva

Muller_4.c

R0.01

The beam is stronger

Real and imaginary roots of the dispersion

relation

R0.01

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.

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

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

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

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

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

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.

Tridiagonal matrix (1)

Based on Gauss elimination

R2

R3

Tridiagonal matrix (2)

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

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.

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.

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

Finite Difference Method (II)

Grid 100x50

Full rectangle

100V

-100V

100V

-100V

Full circle

100V

-100V

100V

-100V

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

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

FEMLAB Environment

Model Navigator

Pre-defined Equations

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

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.

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

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

I. Particle Equations of Motion

- Newton-Lorentz equations of motion

- In finite difference form, the leapfrog method

Fig Schematic leapfrog integration

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

PDP Structure

AC PDP

Discharge in PDP

Striation Profiles in PDP 2D PIC/MCC

Anode

Cathode

- Pressure dependence of striations
- Number of peaks depend on the pressure and

electrode size.

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.

Errors in the Solution by the Shooting Method

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

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

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

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

IF (XP.LE.XL) GO TO 30 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

Programs of Initial Value Problem Shooting

Method for BVP ODE

Sung Jin Kim and Jae Koo Lee

(No Transcript)

(No Transcript)

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)

Second-Order Runge-Kutta Method

By 2nd order RK Method

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

Various Numerical Methods

h0.1

Exact Solution

h0.01

h0.001

Error Estimation

Error estimation

Fourth order Runge-Kutta

Program 9-2 Fourth-order Runge-Kutta Scheme

A first order Ordinary differential equation

y(0)1

Fourth-order RK Method

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

--------------------------------------

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

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

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

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.

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.

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

Solution by the Shooting Method

Errors in the Solution by the Shooting Method

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

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.

Boundary Value Problems of Ordinary Differential

Equations Scharfetter-Gummel method

Sung Soo Yang and Jae Koo Lee

Homepage http//jkl.postech.ac.kr

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.

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,

Program 10-1

For i 10,

Summarizing the difference equations obtained, we

write

Tridiagonal matrix

Solution Algorithm for Tridiagonal Equations (1)

R2

Based on Gauss elimination

R3

Solution Algorithm for Tridiagonal Equations (2)

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

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 )

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

Program 10-3

An eigenvalue problem of ordinary differential

equation

We assume

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

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

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)

Scharfetter-Gummel method

Tridiagonal matrix

- 2D discretized continuity eqn. integrated by the

alternative direction implicit (ADI) method

Scharfetter-Gummel method

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

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)

Scharfetter-Gummel method

In case of

Drift flux

Scharfetter-Gummel method

In case of

Diffusion flux