Plasma Application Modeling Group POSTECH

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 (김성진)

Plasma Application Modeling Group POSTECH