Plasma Application Modeling Group POSTECH

1 / 18

Plasma Application Modeling Group POSTECH - PowerPoint PPT Presentation

Plasma Application Modeling Group POSTECH. Programs of Initial Value Problem &amp; Shooting Method for BVP ODEs (Computational E&amp;M: 490D). Sung Jin Kim, and Jae Koo Lee March 9. 2004. Plasma Application Modeling, POSTECH. Initial Value Problems of Ordinary Differential Equations.

I am the owner, or an agent authorized to act on behalf of the owner, of the copyrighted work described.

PowerPoint Slideshow about 'Plasma Application Modeling Group POSTECH' - kamran

Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author.While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.

- - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - -
Presentation Transcript

Plasma Application

Modeling Group

POSTECH

Programs of Initial Value Problem

& Shooting Method for BVP ODEs

(Computational E&M: 490D)

Sung Jin Kim, and Jae Koo Lee

March 9. 2004

Plasma Application

Modeling, POSTECH

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)

Plasma Application

Modeling, POSTECH

Second-Order Runge-Kutta Method

By 2nd order RK Method

Plasma Application

Modeling, POSTECH

Program 9-1

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;

}

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

2nd order RK

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

h=0.1

Exact Solution:

h=0.01

h=0.001

Plasma Application

Modeling Group

POSTECH

Error Estimation

Error estimation

Fourth order Runge-Kutta

Plasma Application

Modeling, POSTECH

Program 9-2 Fourth-order Runge-Kutta Scheme

A first order Ordinary differential equation

y(0)=1

Fourth-order RK Method

Plasma Application

Modeling, POSTECH

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

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

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

Interval of t for printing ?

1

Number of steps in one printing interval?

100

Maximum t?

5

h=0.01

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

t y

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

0.00000 0.000000e+00

1.00000 1.410686e+00

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(0)=1,

y1(0)=1, y1’(0)=0

* The 4th order RK method for the set of two equations(Nakamura’s book p332)

do{

for( n = 1; n <= 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 <= No_of_eqs; i++ ) ya[i] = y[i];

f( k1, ya, &t_old, &h );

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 );

for( i = 1; i <= No_of_eqs; i++ ) y[i] = y[i] + (k1[i] + k2[i]*2 + k3[i]*2 + k4[i])/6;

}

void f(float k[], float y[], float *t, float *h)

{

k[1] = y[2]**h;

k[2] = -y[1]**h;

/* More equations come here if the number of equations

are greater.*/

return;

}

Plasma Application

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

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.00000e+00 0.00000e+00

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

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.00000e+00 0.00000e+00

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

Plasma Application

Modeling, POSTECH

* Three types of boundary conditions

1) : Dirichlet type boundary condition

1) : Neumann type boundary condition

1) : Mixed type boundary condition

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

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:

We should assume a initial value of z(a).

z(a) is determined for which y(b)=B by secant method.

Plasma Application

Modeling, POSTECH

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 x=0.25, S(0.0)(1)=50, S(0.0)(2)=100

Plasma Application

Modeling, POSTECH

Solution by the Shooting Method

Plasma Application

Modeling, POSTECH

Errors in the Solution by the Shooting Method

Plasma Application

Modeling, POSTECH

Equilibrium Method

y(a)=A and y(b)=B

Ex)

T(0)=0 and T(1.0)=100

Let x=0.25, Ta=0, =2.0

x=0.25 :

0

0

-100

-2.25 1.0 0

1.0 -2.25 1.0

0 1.0 -2.25

T1

T2

T3

=

x=0.50 :

x=0.75 :

Plasma Application

Modeling, POSTECH

Executing Programs for 490D class

A user ID will be made as team name at pam10 computer.

A password is the same as a 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
• ./a

where, a is a changeable execution sentence.

TA’s E-mail: hsus77@postech.ac.kr (고한석), medic@postech.ac.kr (김성진)