Transcript Computational E&M: 490D

Computational Eng./Sci.
CES
EE
Phys/Chem/Bio
Mech/Ch.Eng
ECE490O: Special Topics
in EM-Plasma Simulations
JK LEE (Spring, 2006)
•
•
•
•
ODE Solvers
PIC-MCC
PDE Solvers (FEM and FDM)
Linear & NL Eq. Solvers
ECE490O: ODE
JK LEE (Spring, 2006)
Gonsalves’ lecture notes (Fall 2005)
Programs of Initial Value Problem
& Shooting Method for BVP ODE
Sung Jin Kim and Jae Koo Lee
Plasma Application
Modeling Group
POSTECH
Gonsalves’ lecture notes (Fall 2005)
ODE  IV : y  f ( y, t )
 Euler ( Explicit ) : yn 1  yn  h f n
 Modif .Euler ( Implicit ) :
yn 1  yn  h2  f n 1  (1   ) f n 
iter. for NL
 R  K (2nd Order ) :
yn 1  yn  h2  f n  f ( yn  k1 , tn 1 ) 
 R  K (3rd ) : k1  h f n
k2  h f ( yn  k21 , tn  h2 )
k3  h f ( y n  k n  2 k 2 , t n  h )
yn 1  yn  16 h (k1  4 k2  k3 )
 R  K (4th) : k2  h f ( yn  k21 , tn  h2 )
k3  h f ( yn 
k2
2
, tn  h2 )
k4  h f ( yn  k3 , tn  h)
yn 1  yn  16  k1  2 k2  2 k3  k4 
Leaptrog : yn 1  yn 1  2 h f n
Two  step :
yn  1  yn  12 h f n
2
yn 1  yn  h f n  1
2
Initial Value Problems of Ordinary Differential Equations
Program 9-1 Second-order Runge-kutta Method
d2
d
M 2 y (t )  B y (t )  ky(t )  0 ,
dt
dt
y(0)=1, y’(0)=0
Changing 2nd order ODE to 1st order ODE,
d
y (t )  f ( y, z , t )  z (t )
dt
M
(1)
y(0)=1
d
z (t )  Bz (t )  ky (t )  0
dt
d
B
k
z (t )  g ( y, z , t )  
z (t ) 
y (t )
dt
M
M
(2) z(0)=0
Plasma Application
Modeling, POSTECH
Second-Order Runge-Kutta Method
By 2nd order RK Method
k1  hf ( yn , zn , t )  h  zn
l1  hg ( yn , zn , t )  h(
B
k
zn 
yn )
M
M
k2  hf ( yn  k1, zn  l1, t )  h( zn  l1 )
k
 B

l2  hg( yn  k1 , zn  l1 , t )  h  ( zn  l1 )  ( yn  k1 ) 
M
 M

1
k1  k2 
2
1
z n 1  z n  l1  l2 
2
yn 1  yn 
Plasma Application
Modeling, POSTECH
Program 9-1
/* CSL/c9-1.c Second Order Runge-Kutta Scheme
(Solving the problem II of Example 9.6) */
#include <stdio.h> #include <stdlib.h>
#include <math.h>
/* 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, kstep=0;
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 <= 20; n++ ){
for( kount = 1; kount <= 50; kount++ ){
kstep=kstep+1;
time = h*kstep ;
k1 = h*z;
2nd order RK
l1 = -h*(bm*z + km*y);
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);
result
}
CSL/C9-1 Second Order Runge-Kutta Scheme
t
y
z
0.000000 1.00000e+00 0.00000e+00
0.100000 5.08312e-01 -6.19085e+00
0.200000 6.67480e-02 -2.46111e+00
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
Plasma Application
Modeling, POSTECH
Various Numerical Methods
C1  1, C2 

h=0.1
1.0

Exact solution
Modified Euler
Second order RK
Fourth orde Rk
0.8
 
0.6
2
y  et C1 cos   2 t  C2 sin   2 t

Y(m)
Exact Solution:
0.4
0.2
0.0
-0.2
0.0
0.2
0.4
0.6
0.8
1.0
time
h=0.01
1.0
Exact solution
Euler
Modified Euler
Second order RK
Fourth order Rk
0.8
0.8
0.6
Y(m)
Y(m)
0.6
0.4
0.2
Exact solution
Euler
Modified
Second order RK
Fourth order RK
h=0.001
1.0
0.4
0.2
0.0
0.0
-0.2
0.0
0.2
0.4
0.6
Time
0.8
1.0
-0.2
0.0
0.2
0.4
0.6
0.8
1.0
Time
Plasma Application
Modeling, POSTECH
Error Estimation
Error estimation
Fourth order Runge-Kutta
Exact solution
h=0.001
h=0.005
h=0.01
h=0.05
h=0.1
1.0
0.8
0.0025
0.0020
Error
Y(m)
0.6
0.4
0.2
Euler
Modified Euler
Second order RK
Fourth order RK
2
reference(h )
3
reference(h )
5
reference(h )
0.0030
0.0015
0.0010
0.0005
0.0
0.0000
-0.2
0.0
0.2
0.4
0.6
time
0.8
1.0
0.00
0.01
0.02
0.03
0.04
0.05
h
Plasma Application
Modeling Group
POSTECH
Program 9-2 Fourth-order Runge-Kutta Scheme
A first order Ordinary differential equation
y ' (t )  
y (t )
t 2  y 2 (t )
y(0)=1
Fourth-order RK Method
k1  hf ( yn , tn )
k
h

k2  hf  yn  1 , tn  
2
2

yn 1  yn 
1
k1  2k2  2k3  k4 
6
k
h

k3  hf  yn  2 , tn  
2
2

k4  hf  yn  k3 , tn  h
Plasma Application
Modeling, POSTECH
Program 9-2 Fourth-order Runge-Kutta Scheme
Main algorithm for 4th order RK method at program 9-2
do{ for( j = 1; j <= nstep_pr; j++ ){
t_old = t_new;
t_new = t_new + h;
yn = y;
t_mid = t_old + hh;
yn = y;
k1 = h*fun( yn, t_old );
ya = yn + k1/2; k2 = h*fun( ya, t_mid );
ya = yn + k2/2; k3 = h*fun( ya, t_mid );
ya = yn + k3 ; k4 = h*fun( ya, t_new );
y = yn + (k1 + k2*2 + k3*2 + k4)/6;
}
double fun(float y, float t) {
float fun_v;
fun_v = t*y +1; /* Definition of f(y,t) */
return( fun_v );
}
4th order RK
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
h=0.1
-------------------------------------t
y
-------------------------------------0.00000
0.000000e+00
1.00000
1.410686e+00
2.00000
8.839368e+00
3.00000
1.125059e+02
4.00000
3.734231e+03
5.00000
3.357971e+05
6.00000
8.194354e+07
-------------------------------------Maximum t limit exceeded
Interval of t for printing ?
1
-------------------------------------Number of steps in one printing interval? t
y
100
-------------------------------------Maximum t?
0.00000
0.000000e+00
5
1.00000
1.410686e+00
h=0.01
2.00000
8.839429e+00
3.00000
1.125148e+02
4.00000
3.735819e+03
5.00002
3.363115e+05
--------------------------------------
Plasma Application
Modeling, POSTECH
Program 9-3 4th order RK Method for a Set of ODEs
A second order Ordinary differential equation
y1 ' ' (t )  y1 (t )  0
d
y1 (t )  y2 (t )
dt
y1(0)=1, y1’(0)=0
d
y2 (t )   y1 (t )
dt
* The 4th order RK method for the set of two equations(Nakamura’s book p332)
do{
void f(float k[], float y[], float *t, float *h)
for( n = 1; n <= ns; n++ ){
{
t_old = t_new; /* Old time */
k[1] = y[2]**h;
t_new = t_new + h; /* New time */
k[2] = -y[1]**h;
t_mid = t_old + hh; /* Midpoint time */
/* More equations come here if the number of equations
for( i = 1; i <= No_of_eqs; i++ ) ya[i] = y[i];
are greater.*/
f( k1, ya, &t_old, &h );
return;
for( i = 1; i <= No_of_eqs; i++ ) ya[i] = y[i] + k1[i]/2;
}
f( k2, ya, &t_mid, &h );
for( i = 1; i <= No_of_eqs; i++ ) ya[i] = y[i] + k2[i]/2;
f( k3, ya, &t_mid, &h );
for( i = 1; i <= No_of_eqs; i++ ) ya[i] = y[i] + k3[i];
f( k4, ya, &t_new, &h );
Plasma Application
for( i = 1; i <= No_of_eqs; i++ ) y[i] = y[i] + (k1[i] + k2[i]*2 + k3[i]*2 + k4[i])/6;
Modeling, POSTECH
}
Results of Program 9-3
CSL/C9-3
Fourth-Order Runge-Kutta Scheme
for a Set of Equations
Interval of t for printing ? 1
Interval of t for printing ? 1
Number of steps in one print interval ? 100
Number of steps in one print interval ? 10
Maximum t to stop calculations ? 5.0
Maximum t to stop calculations ? 5.0
h= 0.01
h= 0.1
t
y(1),
y(2), .....
t
y(1),
y(2), .....
0.0000
1.00000e+00
0.00000e+00
0.0000 1.00000e+00 0.00000e+00
1.0000 5.40302e-01 -8.41471e-01
1.0000 5.40303e-01 -8.41470e-01
2.0000 -4.16147e-01 -9.09298e-01
2.0000 -4.16145e-01 -9.09298e-01
3.0000 -9.89992e-01 -1.41120e-01
3.0000 -9.89992e-01 -1.41123e-01
4.0000 -6.53644e-01 7.56802e-01
4.0000 -6.53646e-01 7.56800e-01
5.0000 2.83662e-01 9.58924e-01
5.0000 2.83658e-01 9.58925e-01
Type 1 to continue, or 0 to stop.
6.0000 9.60168e-01 2.79420e-01
1
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.00000e+00 0.00000e+00
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.
Plasma Application
Modeling, POSTECH
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)
KM=K/M
BM=B/M
DO 20 N=1,20
DO 10 KOUNT=1,50
T=T+H
K1=H*Z
L1=-H*(BM*Z + KM*Y)
K2=H*(Z+L1)
L2=-H*(BM*(Z+L1) + KM*(Y+K1))
Y=Y+(K1+K2)/2
Z=Z+(L1+L2)/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
Y=1
PRINT *, 'H=', H
XP=0
HH=H/2
PRINT *
PRINT *,'-------------------PRINT *,'
T
Y'
PRINT *,'---------------------PRINT 82, XP,Y
82 FORMAT( 1X,F10.6, 7X,1PE15.6)
30 DO 40 J=1,I
XB=XP
XP=XP+H
YN=Y
XM=XB+HH
K1=H*FUN(YN,XB)
K2=H*FUN(YN+K1/2,XM)
K3=H*FUN(YN+K2/2,XM)
K4=H*FUN(YN+K3,XP)
Y=YN + (K1+K2*2+K3*2+K4)/6
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/(X*X+Y*Y )
RETURN
END
Executing Programs for the simulation class
An user ID will be made as team name at a PAM computer.
A password is the same as an user ID
You can execute your programs to connect your PC to pam 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.
Plasma Application
Modeling, POSTECH
Boundary Value Problems of
Ordinary Differential Equations &
Scharfetter-Gummel method
Sung Soo Yang and Jae Koo Lee
Plasma
Application
Modeling
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.
* Three types of boundary conditions
dy (a)
dt
1)   0   0 :
1)   0   0 :
1)   0   0 :
y (a)  

Dirichlet type boundary condition
Neumann type boundary condition
Mixed type boundary condition
Plasma Application
Modeling @ POSTECH
Program 10-1
Solve difference equation,
 2 y( x)  y( x)  exp(0.2 x)
With the boundary conditions,
y(0)  1, y(10)   y(10)



x=0
i=0
1
1
2
2
known
y(0)=1






9
9
10
10
y(10)   y(10)
y( x) 
yi 1  2 yi  yi 1
 yi 1  2 yi  yi 1
2
h
 2 y( x)  y( x)  exp(0.2 x)
Especially for i = 1, y(0)  1
2( yi 1  2 yi  yi 1 )  yi  exp(0.2i)
5 y1  2 y2  exp(0.2)  2
Plasma Application
Modeling @ POSTECH
Program 10-1
For i = 10, y(10)   y(10)
y 
y(10)  y(9.5) y(10)   y (10)  y (9)

0.5
0.5
 2 y( x)  y( x)  exp(0.2 x)
 2 y9  4.5 y10  0.5 exp(2)
Summarizing the difference equations obtained, we write
5 y1  2 y2
 exp(0.2)  2
 2 yi 1  5 yi  2 yi 1  exp(0.2xi )
 2 y9  4.5 y10  0.5 exp(2)
Tridiagonal matrix
 5 2
  y1  exp(0.2)  2
 2 5  2
  y   exp(0.2  2) 

 2  


  y3    exp(0.2  3) 
2 5 2

  





  


 2 4.5  y10   0.5 exp(2)  Plasma Application
Modeling @ POSTECH
Solution Algorithm for Tridiagonal Equations (1)
Based on Gauss elimination
 B1
A
 2
 0
 A2
A
 2
 0
C1
B2
A3
R2C1
B2
A3
 A2 R2C1
 R3 A
3
B2
0
B2

A3
 0
0  1   D1 
C2  2    D2 
B3  3   D3 
0  1   R2 D1 
C2  2    D2 
B3  3   D3 
0     R2 D1 
1
A3     A3 
C2  2    D2 
B2 
B2 

 
B3  3   D3 
R2
 A2
 B B1
 1
 A2
 0


 A2
0

 0
 A2
0

 0
A2

 A2 
C1 0  
D1 



1
B1
B
   1 
B2
C2  2    D2 
A3
B3  3   D3 






R2C1
B2
B2  R2C1
A3
R2C1
A3
0
0  1   R2 D1 
D2




C2  2    D2  R2 D1 

B3  3  
D3
 1   R2 D1 
R3C2  2    R3 D2 
B3  R3C2  3   D3  R3 D2 
B3
D3
0
Plasma Application
Modeling @ POSTECH
Solution Algorithm for Tridiagonal Equations (2)
 A2
0

 0
R2C1
A3
0
0  1   R2 D1 
R3C2  2    R3 D2 
B3  3   D3 
A32  R3C23  R3 D2  2 
 2 
D3
 3 
B3
R3
A
( D2  C23 ) , R3  3
A3
B2
1
( D2  C23 )
B2
R2
A2
A21  R2C12  R2 D1  1 
( D1  C12 ) , R2 
A2
B1
1
1
 1  ( D1  C12 )  ( D1  C12 )
B1
B1
Plasma Application
Modeling @ POSTECH
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 <= n; i++ )
{
Recurrently calculate the equations in
r = a[i]/b[i - 1];
increasing order of i until i=N is reached
b[i] = b[i] - r*c[i - 1];
d[i] = d[i] - r*d[i - 1];
}
Calculate the solution for the last
d[n] = d[n]/b[n];
unknown by   D B
N
for ( i = n - 1; i >= 1; i-- )
{
d[i] = (d[i] - c[i]*d[i + 1])/b[i];
}
return;
N
N
Calculate the following equation in
decreasing order of i
 i  ( Di  Cii 1 ) Bi
}
Plasma Application
Modeling @ POSTECH
Program 10-1
/* CSL/c10-1.c Linear Boundary Value Problem */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
/* a[i], b[i], c[i], d[i] : a(i), b(i), c(i), and d(i) n: number of grid points */
int main()
{
int i, n;
float a[20], b[20], c[20], d[20], 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 <= n; i++ )
{
x = i;
a[i] = -2; b[i] = 5; c[i] = -2;
d[i] = exp( -0.2*x );
}
 2 yi 1  5 yi  2 yi 1  exp(0.2xi )
Plasma Application
Modeling @ POSTECH
Program 10-1
5 y1  2 y2
 exp(0.2)  2
d[1] = d[1] + 2;
d[n] = d[n]*0.5;
 2 y9  4.5 y10  0.5 exp(2)
b[n] = 4.5;
trdg( a, b, c, d, n );
d[0] = 1; /* Setting d[0] for printing purpose */
printf( "\n Grid point no. Solution\n" );
for ( i = 0; i <= n; i++ )
{
CSL/C10-1 Linear Boundary Value Problem
Grid point no.
Solution
printf( " %3.1d %12.5e \n", i, d[i] );
0
1.00000e+00
}
1
8.46449e-01
exit(0);
}
2
3
4
5
6
7
8
9
10
7.06756e-01
5.85282e-01
4.82043e-01
3.95162e-01
3.21921e-01
2.59044e-01
2.02390e-01
1.45983e-01
7.99187e-02
Plasma Application
Modeling @ POSTECH
Program 10-3
An eigenvalue problem of ordinary differential equation

1 d 

m d
p
(
x
)
x
f
(
x
)

  q( x) f ( x)  v( x) f ( x) , 0  x  H
m
x dx 
dx

We assume
p( x)  1, q( x)  0, v( x)  1, m  0

i-1
a

i
b

i+1
d2
 2 f ( x )  f ( x )
dx
Ai fi (t1)  Bi fi (t )  Ci fi (t1)  (t 1)Gi fi (t 1)
f i 1  f i f i  f i 1

f  2 f i  f i 1
f (b)  f (a)
h
h
f  

 i 1
h
h
h2
 fi (t1)  2 fi (t )  fi (t1)  (t 1) fi (t 1)
Plasma Application
Modeling @ POSTECH
Program 10-3 (Inverse Power method)
Step1 : fi (0) for all i and ( 0 ) are set to an arbitrary initial guess.
(1)
(1)
(1)
( 0) ( 0 )
(1)
Step2 :  fi 1  2 fi  fi 1  fi  is solved by the tridiagonal solution for f i
Step3 : The next estimate for  is calculated by the equation:
(1)  ( 0)
( 0 ) (1)
G
f
 i i fi
i
G f
i i
(1)
f i (1)
i
( 2)
( 2)
( 2)
(1) (1)
Step4 :  fi 1  2 fi  fi 1   fi is solved by the tridiagonal solution for fi
( 2)
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.
(t 1) / ( t )  1  Criterion for convergence
Plasma Application
Modeling @ POSTECH
Program 10-3
main()
{
•••••
while (TRUE)
{
k = 0;
n = 10;
ei = 1.5;
Step 1
it = 30;
ep = 0.0001;
for (i=1; i <= n; i++)
{
as[i] = -1.0;
bs[i] = 2.0;
cs[i] = -1.0;
f[i] = 1.0;
ds[i] = 1.0;
}
printf("\n It. No. Eigenvalue\n");
while (TRUE) {
k = k+1;
for (i=1; i<=n; i++) fb[i] = f[i]*ei;
for (i=1; i<=n; i++) {
a[i] = as[i]; b[i] = bs[i]; c[i] = cs[i];
d[i] = ds[i]*fb[i];
}
Step 2,4
trdg(a, b, c, d, n);
sb = 0; s = 0;
for (i=1; i <= n; i++) {
f[i] = d[i];
s = s + f[i]*f[i];
Step 3
sb = sb + f[i]*fb[i];
}
eb = ei; ei = sb/s;
printf("%3.1d %12.6e \n", k, ei);
if (fabs(1.0-ei/eb) <= ep) break; Step 6
if (k >it) {
printf("Iteration limit exceeded.\n");
break;
}
Plasma Application
}
Modeling @ POSTECH
Program 10-3
z = 0;
for (i=1; i <= n; i++)
if (fabs(z) <= fabs(f[i])) z = f[i];
for (i=1; i <= n; i++) f[i] = f[i]/z;
eigen = ei;
printf("Eigenvalue = %g \n", eigen);
printf("\n Eigenfunction\n");
printf("i f(i)\n");
for (i=1; i<=n; i++)
printf("%3.1d %12.5e \n", i, f[i]);
printf("------------------------------------\n");
printf("Type 1 to continue, or0 to stop. \n");
scanf("%d", &kstop);
if (kstop != 1) exit (0);
}
}
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.00000e+00
6 1.00000e+00
7 9.19018e-01
8 7.63594e-01
9 5.46290e-01
10 2.84692e-01
-----------------------------------Type 1 to continue, or0 Plasma
to stop.Application
Modeling @ POSTECH
Scharfetter-Gummel method
n
  Γ  S
t
Γ   Dn  sgn(q)nE
Ai, j ni*1, j  Bi, j ni*, j  Ci, j ni*1, j  Di, j
Ai', j nik, j 11  Bi', j nik, j 1  Ci', j nik, j 11

Di', j
Tridiagonal matrix
• 2D discretized continuity eqn. integrated by the alternative direction implicit
(ADI) method
ni*, j  nik, j
t / 2
nik, j 1  ni*, j
t / 2


x*,i 0.5, j  x*,i 0.5, j
x
yk,i ,1j 0.5  yk,i,1j 0.5
y
 Sik, j 
 Sik, j 
yk,i, j 0.5  yk,i, j 0.5
y
x*,i 0.5, j  x*,i 0.5, j
x
x,i0.5,j   x,i,j ni, j   x,i, j ni1, j Scharfetter-Gummel method
Plasma Application
Modeling @ POSTECH
Scharfetter-Gummel method
i  1
i  1
2


ni 1
Γ
i
1
2
Flux is constant between half grid
points and calculated at the grid point
2
ni 1
ni
 ( sign) i ni Ei 
zi  1

ni Di
x

 ni Di
 ni Di
x
x
2
Analytic integration between i
and i+1 leads to
Flux i  12 is a linear combination of
ni and ni 1
zi  1

2
ni Di  ni Di
 Γi  1
2
x
x
1
zi  1 ni 1 Di 1  Γ i  1
zi  1
2
2
x
e 2
1
zi  1 ni Di  Γ i  1
2
2
x
Plasma Application
Modeling @ POSTECH
Scharfetter-Gummel method
where, zi  1  
(sign)i  1
Di  1
2
Γi  1 
2
i 
zi  1
x
zi  1
zi  1
e
2
x e
(Vi 1 Vi ) 
(sign)i  1
Di  1
2
ni Di e
2
2
zi  1
Di
zi  1
2
 ni 1Di 1
2
1
2
e
1
zi  1
2
2
Ei  1 xi
2
2
  i ni   i ni 1
i 
zi  1
2
Di 1
x e zi  12 1
The main advantage of SG scheme is that is provides numerically stable
estimates of the particle flux under all conditions
zi  1  1
2
zi  1  1
2
(Big potential difference) (small potential difference)
Plasma Application
Modeling @ POSTECH
Scharfetter-Gummel method
Γi  1
2
ni exp zi  1   ni 1
n expk   ni 1
2


i  1 Ei  1  i
i  1 Ei  1
2
2
2
2
expk   1
exp zi  1   1
2

In case of zi  1  1
2
lim Γi  1
k 
2
ni expk   ni 1
exp(k )
 lim
i  1 Ei  1
k 
2
2 exp( k )
expk   1
ni  ni 1 exp(k )
 lim
i  1 Ei  1  ni i  1 Ei  1
k 
2
2
2
2
1  exp(k )
lim Γi  1
k 
2
ni expk   ni 1
 lim
i  1 Ei  1  ni 1i  1 Ei  1
k 
2
2
2
2
expk   1
Drift flux
Plasma Application
Modeling @ POSTECH
Scharfetter-Gummel method
In case of zi  1  1
2
lim Γ i  1
k 0
2
ni expk   ni 1
ni expk   ni 1 Di  12
 lim
i  1 Ei  1  lim
k
k 0
k

0
2
2
expk   1
expk   1
xi
d
k ni expk   ni 1 
Di  1
2

lim dk
d
xi k 0
expk   1
dk

ni e(1  k ) expk   ni 1
xi k 0
exp(k )
Di  1
 Di  1
2
2
lim
ni 1  ni
xi
Diffusion flux
Plasma Application
Modeling @ POSTECH
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
y1 ' '  f ( x, y, y' ), y(a)=A and y(b)=B
* Three types of boundary conditions
dy (a)
dt
1)   0   0 :
1)   0   0 :
1)   0   0 :
y (a)  

Dirichlet type boundary condition
Neumann type boundary condition
Mixed type boundary condition
Plasma Application
Modeling, POSTECH
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 x=b 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.
Plasma Application
Modeling, POSTECH
Shooting Method for Boundary Value Problem ODEs
Rewrite the second-order ODE as two first-order ODEs:
y' (t )  z,
y (a)  A
z ' (t )  f ( x, y, z ),
z (a)  y ' (a)  ?
We should assume a initial value of z(a).
z(a) is determined for which y(b)=B by secant method.
B  (B) 
B
(n)
( A' )( n1)  ( A' )( n )
( A' )
( n 1)

 ( B)( n1)

 Slope
( A' )( n )  ( A' )( n1)
 ( A' )
(n)
B  ( B) 

( n)
(n)
Slope
Plasma Application
Modeling, POSTECH
Example of Shooting Method
Ex)
T ' ' 2T  0,
T(0)=0 and T(1.0)=100
Sol) Rewrite the second-order ODE as two first-order ODEs:
T '  S,
T(0)=0.0
S '   2T ,
S(0)=T’(0)
By modified Euler method
P
Ti 1
 Ti  tSi
1
Ti C1  Ti  t Si  SiP1 
2
SiP1  Si  t 2Ti
1
SiC1  Si  t 2 Ti  Ti P1 
2
Let x=0.25, S(0.0)(1)=50, S(0.0)(2)=100
Plasma Application
Modeling, POSTECH
Solution by the Shooting Method
T (1.0)( 2)  T (1.0)(1) 170.507813 85.253906
Slope

 1.705078
( 2)
(1)
S (0.0)  S (0.0)
100.0  50.0
S (0.0)
( 3)
 S (0.0)
100.0  T (1.0)   58.648339

( 2)
( 2)
Slope
Plasma Application
Modeling, POSTECH
Errors in the Solution by the Shooting Method
Plasma Application
Modeling, POSTECH
Equilibrium Method
y1 ' ' P( x) y'Q( x) y  F ( x),
y(a)=A and y(b)=B
 x 
 x 
2
1

P
y


2


x
Q
y

Pi  yi 1  x 2 Fi

1 
i  i 1
i
i
2 
2 



Ex) T ' ' 2T   2Ta ,

T(0)=0 and T(1.0)=100
Ti 1  (2   2x2 )Ti  Ti 1   2x2Ta
Let x=0.25, Ta=0, =2.0
Ti 1  2.25Ti  Ti 1  0
x=0.25 : T0  2.25T1  T2  0, T0  0
x=0.50 : T1  2.25T2  T3  0
x=0.75 : T2  2.25T3  T4  0, T4  100.0
-2.25 1.0 0
1.0 -2.25 1.0
0
1.0 -2.25
0
T1
0
T2 =
-100
T3
Plasma Application
Modeling, POSTECH