Transcript 슬라이드 1

Boundary Value Problems of
Ordinary Differential Equations &
Scharfetter-Gummel method
(Computational E&M: 490D)
Sung Soo Yang, and Jae Koo Lee
March 23, 2004
Plasma
Application
Modeling
Homepage : http://jkl.postech.ac.kr
Boundary value problems
Type of Problems
Advantages
Disadvantages
Shooting method
An existing program for
Trial-and-error basis. Application
initial value problems may is limited to a narrow class of
be used
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

 A2 
C1 0  
D1 



1
B1
B
   1 
B2
C2  2    D2 
A3
B3  3   D3 






 A2
0

 0
0  1   R2 D1 

B2    
D2
B2  R2C1 C2  2    D2  R2 D1 

A3
B3  3  
D3
 A2
0

 0
R2C1
R2C1
A3
0
 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