From Documentation
Jump to: navigation, search
Note: This page hasn't been updated for a long time, and is most likely outdated.
If you need it to be updated, please submit a ticket to help@sharcnet.ca with the link to the page.

RUNGE-KUTTA for ORDINARY DIFFERENTIAL EQUATIONS

Nick Chepurniy, Ph.D.
e-mail: nickc@sharcnet.ca
HPC Consultant - SHARCNET http://www.sharcnet.ca
Compute Canada http://www.computecanada.org


Overview

This tutorial illustrates the Runge-Kutta method for solving systems of Ordinary Differential Equations (ODE) and Partial Differential Equations (PDE). Examples of FIRST-ORDER, SECOND-ORDER and SIXTH-ORDER ODEs are given and solved using a c-program. Results are discussed. PDE in two and three independent variables can be discretized in all but one direction and then solved as a system of ODEs.





A FIRST-ORDER ordinary differential equation (ODE) has the first derivative as its highest derivative and it usually can be cast in the following form:


         y'(x) = F(x,y(x))                                           (1)

where x is the independent variable, y(x) is the dependent variable and y'(x) is the first derivative.

The dependent variable at x+h, i.e. y(x+h), when using fourth-order Runge-Kutta of the first-order differential equation (1), is given by:

       y(x+h) = y(x) + (k1 +2*k2 + 2*k3 + k4)/6                      (3)

where

       k1 = h*F(x,y(x))

       k2 = h*F(x+h/2,y(x)+k1/2)

       k3 = h*F(x+h/2,y(x)+k2/2)

       k4 = h*F(x+h,y(x)+k3)

Therefore, to start the solution we need to have the Initial Conditions:

       y(x) = y           at x = x
               0                  0

This will then yield the dependent variable at

       x = x  + h
            0

Once we know

       y(x  + h)
          0

we proceed to calculate

       y(x  + 2*h)        
          0

and so on.

FIRST ORDER ODE EXAMPLE

In this section we will present a c program which solves ODE's of any order. However, we will start using it with an ODE of FIRST ORDER. Furthermore, to ensure that the program works we will build the ODE using a known dependent variable. This will allow us to verify our results against the exact solution.

Assume that

      y(x) = exp(-2*x)      (EXACT SOLUTION)

is the dependent variable. Then the first derivative of y(x) with respect to x is given by

      y'(x) = -2*exp(-2*x)

But the expression exp(-2*x) - 3*y(x) is also equal to -2*exp(-2*x). Therefore, we arrive at the following ODE:

  y'(x) = exp(-2*x) - 3*y(x)

and for the Initial Conditions we evaluate y(x) at x = 0 which can be stated as

  y(0) = 1
The RUNGE-KUTTA c-program
/*
 
    Assume y(t) = exp(-2*t) then y(0) = 1.0 and dy/dt = -2*y(t) so we can use following
    FIRST-ORDER ODE for testing:
 
           dy/dt = exp(-2*t) - 3*y(t)
            y(0) = 1.0
 
*/
 
 
#include <stdio.h>
#include <math.h>
 
#define N     1            /* number of first order equations */
#define dist  0.005        /* stepsize in t*/
#define MAX   1.0          /* max for t */
 
FILE *output;              /* internal filename */
 
void runge4(double x, double y[], double step);   /* Runge-Kutta function */
double f(double x, double y[], int i);            /* function for derivatives */
 
main()
{
  double t, y[N], exactf,perc_diff;
  int j;
 
  output=fopen("output1.dat005", "w");            /* external filename */
 
  t   = 0.0;
  y[0]= 1.0;                                      /* initial position */
  exactf = exp (-2.0*t);
 
  fprintf(output, "%f\t%f\t%f\n", t, y[0],exactf);
 
  for (j=1; j*dist<=MAX ;j++)                     /* time loop */
     {
       t=j*dist;
       runge4(t, y, dist);
       exactf = exp (-2.0*t);
       perc_diff = (exactf-y[0])*100.0/exactf;
       fprintf(output, "%f\t%f\t%lf\t%lf\n", t, y[0],exactf,perc_diff);
     }
 
  fclose(output);
}
void runge4(double x, double y[], double step)
{
  double h=step/2.0,            /* the midpoint */
  t1[N], t2[N], t3[N],          /* temporary storage arrays */
  k1[N], k2[N], k3[N],k4[N];    /* for Runge-Kutta */
  int i;
 
  for (i=0;i<N;i++) t1[i]=y[i]+0.5*(k1[i]=step*f(x, y, i));
  for (i=0;i<N;i++) t2[i]=y[i]+0.5*(k2[i]=step*f(x+h, t1, i));
  for (i=0;i<N;i++) t3[i]=y[i]+ (k3[i]=step*f(x+h, t2, i));
  for (i=0;i<N;i++) k4[i]= step*f(x+step, t3, i);
 
  for (i=0;i<N;i++) y[i]+=(k1[i]+2*k2[i]+2*k3[i]+k4[i])/6.0;
}
 
double f(double x, double y[], int i)
{
  if (i==0) return(exp(-2*x) - 3*y[0]);     /* derivative of first equation */
  if (i!=0) exit(2);
}

Above c-program was compiled using the command:

pathcc  RK_EX_1.c  -lm

where RK_EX_1.c is the filename containing the c-program. It produced following output written to file output1.dat005:

0.000000        1.000000        1.000000
0.005000        0.990001        0.990050        0.004963
0.010000        0.980102        0.980199        0.009901
0.015000        0.970302        0.970446        0.014814
0.020000        0.960600        0.960789        0.019703
0.025000        0.950996        0.951229        0.024567
0.030000        0.941488        0.941765        0.029407
0.035000        0.932075        0.932394        0.034223
0.040000        0.922756        0.923116        0.039015
0.045000        0.913531        0.913931        0.043783
0.050000        0.904398        0.904837        0.048528
...
0.985000        0.138587        0.139457        0.623438
0.990000        0.137206        0.138069        0.625292
0.995000        0.135838        0.136695        0.627136
1.000000        0.134484        0.135335        0.628970

The first column is the independent variable x ranging from 0.0 to 1.0 in steps of 0.005. The second column is the dependent variable, y(x), computed by the c-program while column 3 is the EXACT SOLUTION. Finally, column 4 is the percentage error between the computed and EXACT solution.

SECOND ORDER ODE EXAMPLE

Let us continue with the dependent variable y(x) = exp(-2*x) and calculate the second order derivative y"(x) = 4*exp(-2*x). Since the expression y'(x) + 5*y(x) + exp(-2*x) is also equal to 4*exp(-2*x) we arrive at the following SECOND-ORDER ODE for our example:

         y"(x) = y'(x) + 5*y(x) + exp(-2*x)

To use the c-program we need to convert this SECOND-ORDER ODE to an equivalent system of two equations of FIRST-ORDER ODE's:

         y0'(x) = y1(x)
         y1'(x) = y1(x) + 5*y0(x) + exp(-2*x)

Here, we have set y(x) = y0(x) and introduced a new dependent variable y1(x)=y'(x). The INITIAL CONDITIONS would be:

         y0(0) =  1
         y1(0) = -2

If you make a copy of the earlier c-program and call it RK_EX_2.c then you must make the following changes for it to solve the above INITIAL VALUE PROBLEM:

#define N     2            /* number of first order equations */
...
  y[0]=  1.0;              /* initial conditions */
  y[1]= -2.0;
...
double f(double x, double y[], int i)
{
  if (i==0) return(y[1]);                          /* derivative of first  equation */
  if (i==1) return(y[1]+5.0*y[0]+exp(-2.0*x));     /* derivative of second equation */
  if (i > 1) exit(2);
}

Compile the RK_EX_2.c program using the command:

       pathcc  RK_EX_2.c  -lm

Execute the program by typing:

       ./a.out

You should get the following output written to file output2.dat005:

0.000000        1.000000        1.000000
0.005000        0.990050        0.990050        0.000013
0.010000        0.980198        0.980199        0.000051
0.015000        0.970444        0.970446        0.000115
0.020000        0.960787        0.960789        0.000206
0.025000        0.951226        0.951229        0.000324
0.030000        0.941760        0.941765        0.000471
0.035000        0.932388        0.932394        0.000647
0.040000        0.923108        0.923116        0.000852
0.045000        0.913921        0.913931        0.001087
0.050000        0.904825        0.904837        0.001354
...
0.985000        0.132766        0.139457        4.797496
0.990000        0.131277        0.138069        4.919336
0.995000        0.129800        0.136695        5.044162
1.000000        0.128336        0.135335        5.172045

Note that the percentage error apearing in column 4 for the SECOND-ORDER ODE example is larger than that from the FIRST-ORDER ODE example. To reduce this error we could specify a smaller step size by making following changes to the RK_EX_2.c program:

#define dist  0.001        /* stepsize in t*/
...
  output=fopen("output2.dat001", "w");                     /* external filename */

Recompiling RK_EX_2.c and executing ./a.out will produce the output file output2.dat001:

0.000000        1.000000        1.000000
0.001000        0.998002        0.998002        0.000000
0.002000        0.996008        0.996008        0.000000
0.003000        0.994018        0.994018        0.000001
0.004000        0.992032        0.992032        0.000002
0.005000        0.990050        0.990050        0.000003
...
0.985000        0.138113        0.139457        0.963341
0.986000        0.137831        0.139178        0.968187
0.987000        0.137549        0.138900        0.973056
0.988000        0.137267        0.138623        0.977949
0.989000        0.136986        0.138346        0.982866
0.990000        0.136705        0.138069        0.987807
0.991000        0.136425        0.137793        0.992771
0.992000        0.136146        0.137518        0.997760
0.993000        0.135867        0.137243        1.002773
0.994000        0.135589        0.136969        1.007810
0.995000        0.135311        0.136695        1.012872
0.996000        0.135034        0.136422        1.017958
0.997000        0.134757        0.136150        1.023069
0.998000        0.134481        0.135878        1.028205
0.999000        0.134205        0.135606        1.033365
1.000000        0.133930        0.135335        1.038551

SIXTH ORDER ODE EXAMPLE

Let us build a sixth order ODE using the same technique we used earlier for the SECOND order ODE. Let the dependent variable y(x) be denoted by y0 = exp(-2*t) + sin(2*t). The details can be found in the listing of the c program RK_EX_6.c (i.e. the INITIAL CONDITIONS and the system of ODE's corresponding to our sixth order ODE. A MAPLE program was used to set up the system of ODE's since the calculations are rather lengthy. The program was compiled using:

pathcc  RK_EX_6.c  -lm

and when it executed it produced following output:

./a.out
 
Starting 6-th order ...
 
   100000000
   200000000
   300000000
   400000000
   500000000
   600000000
   700000000
   800000000
   900000000
  1000000000
 
Normal End of Job
 
 - See output in file:
 
       output6.dat002
 
Number of seconds (elapsed time): 781

The detailed output is in file output6.dat002:

0.000000        1.000000        1.000000
0.100000        1.017400        1.017400       -0.000000
0.200000        1.059738        1.059738        0.000000
0.300000        1.113454        1.113454        0.000000
0.400000        1.166685        1.166685        0.000000
0.500000        1.209350        1.209350        0.000000
0.600000        1.233233        1.233233        0.000000
0.700000        1.232047        1.232047        0.000000
0.800000        1.201470        1.201470        0.000000
0.900000        1.139147        1.139147        0.000000
1.000000        1.044633        1.044633        0.000000
Number of seconds (elapsed time): 781

Finally, we list the whole c program in file "RK_EX_6.c":

/*
 
    Assume
 
            y0(t) = exp(-2*t) + sin(2*t)
 
    then using a MAPLE program we got:
 
> y0:=exp(-2*t)+sin(2*t);
                      y0 := exp(-2 t) + sin(2 t)
 
> y1:=diff(y0,t);
 
                   y1 := -2 exp(-2 t) + 2 cos(2 t)
 
> y2:=diff(y1,t);
 
                    y2 := 4 exp(-2 t) - 4 sin(2 t)
 
> y3:=diff(y2,t);
 
                   y3 := -8 exp(-2 t) - 8 cos(2 t)
 
> y4:=diff(y3,t);
 
                   y4 := 16 exp(-2 t) + 16 sin(2 t)
 
> y5:=diff(y4,t);
 
                  y5 := -32 exp(-2 t) + 32 cos(2 t)
 
 
> # INITIAL CONDITIONS:
 
> y0_t0:=value(Eval(y0,t=0));
 
                              y0_t0 := 1
 
> y1_t0:=value(Eval(y1,t=0));
 
                              y1_t0 := 0
 
> y2_t0:=value(Eval(y2,t=0));
 
                              y2_t0 := 4
 
> y3_t0:=value(Eval(y3,t=0));
 
                             y4_t0 := 16
 
> y5_t0:=value(Eval(y5,t=0));
 
                              y5_t0 := 0
> # SYSTEM of ODE:
 
> Dy0:=y1;
 
                   Dy0 := -2 exp(-2 t) + 2 cos(2 t)
 
> Dy1:=y2;
 
                   Dy1 := 4 exp(-2 t) - 4 sin(2 t)
 
> Dy2:=y3;
 
                   Dy2 := -8 exp(-2 t) - 8 cos(2 t)
 
> Dy3:=y4;
 
                  Dy3 := 16 exp(-2 t) + 16 sin(2 t)
 
> Dy4:=y5;
 
                  Dy4 := -32 exp(-2 t) + 32 cos(2 t)
 
> # Must have RHS of Dy5 such that Dy5 == DDy5 !
 
> Dy5:=-y5+y4-y3+y2-y1+y0+exp(-2*t)+26*cos(2*t)-77*sin(2*t);
 
                  Dy5 := 64 exp(-2 t) - 64 sin(2 t)
 
> DDy5:=diff(y5,t);
 
                  DDy5 := 64 exp(-2 t) - 64 sin(2 t)
 
> ZERO:=DDy5-Dy5;
 
                              ZERO := 0
 
*/
 
#include <stdio.h>
#include <math.h>
#include <time.h>
 
 
#define N               6              /* number of first order equations */
#define dist            0.000000001    /* stepsize in t*/
#define MAX             1.0            /* max for t */
#define PFREQ   100000000              /* Print frequency */
 
FILE *output;              /* internal filename */
 
void runge4(double x, double y[], double step);   /* Runge-Kutta function */
double f(double x, double y[], int i);            /* function for derivatives */
 
main()
{
  time_t sec_start, sec_finish, mytime;
  double t, y[N], exactf, perc_diff;
  int j;
 
  sec_start =  time(NULL);
  printf("\nStarting 6-th order ...\n\n");
  fflush( stdout );
 
  output=fopen("output6.dat002", "w");                     /* external filename */
 
  t   =   0.0;
 
  y[0]=   1.0;              /* initial conditions */
  y[1]=   0.0;
  y[2]=   4.0;
  y[3]= -16.0;
  y[4]=  16.0;
  y[5]=   0.0;
 
  exactf = exp (-2.0*t) + sin(2*t);
 
  fprintf(output, "%f\t%f\t%f\n", t, y[0],exactf);
 
 
  for (j=1; j*dist<=MAX ;j++) /* time loop */
     {
       t=j*dist;
       runge4(t, y, dist);
       if ( (j % PFREQ) == 0 )
         {
           printf("%12d \n",j);
           exactf = exp (-2.0*t) + sin(2*t);
           perc_diff = (exactf-y[0])*100.0/exactf;
           fprintf(output, "%f\t%f\t%lf\t%lf\n", t, y[0],exactf,perc_diff);
         }
     }
 
 
  printf("\nNormal End of Job\n\n - See output in file:\n\n       output6.dat002\n\n");
 
  sec_finish =  time(NULL);
 
  mytime = sec_finish - sec_start;
 
  printf ("Number of seconds (elapsed time): %ld \n", mytime);
  fprintf (output,"Number of seconds (elapsed time): %ld \n", mytime);
 
  fclose(output);
}
 
void runge4(double x, double y[], double step)
{
  double h=step/2.0,         /* the midpoint */
  t1[N], t2[N], t3[N],       /* temporary storage arrays */
  k1[N], k2[N], k3[N],k4[N]; /* for Runge-Kutta */
  int i;
 
  for (i=0;i<N;i++) t1[i]=y[i]+0.5*(k1[i]=step*f(x, y, i));
  for (i=0;i<N;i++) t2[i]=y[i]+0.5*(k2[i]=step*f(x+h, t1, i));
  for (i=0;i<N;i++) t3[i]=y[i]+ (k3[i]=step*f(x+h, t2, i));
  for (i=0;i<N;i++) k4[i]= step*f(x+step, t3, i);
 
  for (i=0;i<N;i++) y[i]+=(k1[i]+2*k2[i]+2*k3[i]+k4[i])/6.0;
}
 
 
double f(double x, double y[], int i)
  {
    if (i==0) return(y[1]);                              /* derivative of first  equation */
    if (i==1) return(y[2]);                              /* derivative of second equation */
    if (i==2) return(y[3]);                              /* derivative of third  equation */
    if (i==3) return(y[4]);                              /* derivative of fourth equation */
    if (i==4) return(y[5]);                              /* derivative of fifth  equation */
    if (i==5) return(-y[5]+y[4]-y[3]+y[2]-y[1]+y[0]+exp(-2.0*x)+26*cos(2*x)-77*sin(2*x) );
       /* derivative of sixth  equation */
  }

RUNGE-KUTTA for PARTIAL DIFFERENTIAL EQUATIONS

Formulation of Initial Value Problem (IVP)

For Partial Differential Equations (PDE) the dependent variable is a function of two or more independent variables. We will consider examples where time t and the coordinate variables x are the independent variables. If the dependent variable is a function of two such independent variables then we call them TWO dimensional, etc ...

For example:

      u(t,x) is 2-D
      v(t,x,y) is 3-D
      w(t,x,y,z) is 4-D


The technique we will use for solving partial differential equations with such dependent variable requires that we first discretize the dependent variables in all but one dimension, and then integrate the semi-discrete problem as a system of ODEs.

To illustrate the technique we will start by solving a 2-D PDE with a know exact solution and specify the Initial Conditions (IC) and Boundary Conditions (BC).

This will give us the advantage of comparing our numerical solution with the EXACT closed form solution.

Here again we will use MAPLE to carry out the calculations.

Consider the 2-D Heat Conduction Equation:


    U  = ThD*U                                                    (1)
     t        xx

where ThD is a constant (Thermal Diffusivity). The IC is:

    U(x,0) = sin(m*pi*x/L)  for m=positive integer                (2)

and the BC are:

    U(0,t) = 0                                                    (3)
    U(L,t) = 0                                                    (4)


Following MAPLE program verifies that the EXACT closed form solution:

                                2   2
                               m  pi  ThD t      m pi x
              u(x, t) := exp(- ------------) sin(------)
                                     2             L
                                    L

satisfies the PDE (1), IC (2) and BC (3) and (4) for any positive integer m. This Initial Value Problem (IVP) models the following physical problem: We have a circular bar of length L which at time t=0 has the initial temperature distribution given by eq. (2). The lateral side of the rod are insulated, therefore heat flows only in the x-direction. The ends of the circular bar (i.e. at x=0 and x=L) are kept at the temperatures specified by the B.C. (3) and (4).


> with(Student[Calculus1]):
 
> with(Student[Calculus1]):
> m=1;
 
                                m = 1
 
>
>
> u(x,t) := exp(-m^2*pi^2*ThD*t/L^2)*sin(m*pi*x/L);
 
                                2   2
                               m  pi  ThD t      m pi x
              u(x, t) := exp(- ------------) sin(------)
                                     2             L
                                    L
 
> u_t:=diff(u(x,t),t);
 
                                    2   2
                   2   2           m  pi  ThD t      m pi x
                  m  pi  ThD exp(- ------------) sin(------)
                                         2             L
                                        L
         u_t := - ------------------------------------------
                                       2
                                      L
 
> u_x:=diff(u(x,t),x);
 
                           2   2
                          m  pi  ThD t      m pi x
                    exp(- ------------) cos(------) m pi
                                2             L
                               L
             u_x := ------------------------------------
                                     L
 
> u_xx:=diff(u_x,x);
 
 
                            2   2
                           m  pi  ThD t      m pi x   2   2
                     exp(- ------------) sin(------) m  pi
                                 2             L
                                L
           u_xx := - --------------------------------------
                                        2
                                       L
 
> LHS:=u_t:
> RHS:=ThD*u_xx:
> ZERO:=LHS-RHS;
 
                              ZERO := 0
 
> u_t0:=sin(n*pi*x/L);
 
                                     n pi x
                         u_t0 := sin(------)
                                       L
 
> u_x0:=0;
 
                              u_x0 := 0
 
> u_xL:=0;
 
                              u_xL := 0
 
>


We will compare this EXACT closed form solution with the Numerical results obtained using the Runge-Kutta method. Here is how the Numerical solution is derived:

Eq. (1) will be discretized with respect to the variable x using the second-order finite difference approximation:

                    U(x+h,t) -2*U(x,t) + U(x-h,t)
       U_xx(x,t) = -------------------------------                  (5)
                               h**2

and first-order finite difference approximation (if needed):

                    U(x+h,t) - U(x-h,t)
        U_x(x,t) = ---------------------                            (6)
                            2h


Assume L=1. Choose a uniform grid xi, 0<=i<=N with spacing h=1/N such that xi=i*h. Let Ui[t] be the value of U(xi,t). Say, N=10, then


    u0[t]  u1[t]  u2[t]  u3[t]  u4[t]  u5[t]  u6[t]  u7[t]  u8[t]  u9[t]  u10[t]
 
     .------.------.------.------.------.------.------.------.------.------.
 
     x=0                                                                  x=1
    x0     x1     x2     x3     x4     x5     x6     x7     x8     x9     x10

The BC (2) and (3) imply that u0(t)=0 and uN(t)=0 for all times t. For the interior points (i.e. for 0 < i < N) we use:

       u[i]' = ThD*(u[i+1]-2*u[i]+u[i-1])/hx2

where ' denotes time derivative and hx2=h**2.

PROGRAM LISTING

/*
  2-D PDE  - Heat Conduction Equation
*/
 
#include <stdio.h>
#include <math.h>
#include <time.h>
 
#define N                  10             /* number of interval in x-direction */
#define N1                 N+1            /* number of first order equations */
#define delta_t            0.000001       /* stepsize in t*/
#define MAX                5.0            /* max for t */
#define PFREQ              500000         /* Print frequency */
 
FILE *output;              /* internal filename */
 
void runge4(double t, double hx, double ThC, double u[], double del_t); /* Runge-Kutta function */
double f(double t, double hx, double ThC, double u[], int i);       /* function for derivatives */
 
main()
{
  time_t sec_start, sec_finish, mytime;
  double t, x, hx, u[N1], del_t, exactf, perc_diff;
  double ThC, pi, L, L2, steps;
  int i, j, m, m2;
 
  output=fopen("output6.dat010", "w");                     /* external filename */
 
  sec_start =  time(NULL);
  printf("\nStarting 2-D PDE Runge-Kutta ...\n\n");
  fprintf(output,"\nStarting 2-D PDE Runge-Kutta ...\n\n");
  fflush( stdout );
 
  t = 0.0;
  hx = 1.0/(double) N;
  del_t = delta_t;
 
  steps = ((double) MAX)/delta_t;
  m = 1;
  m2 = m*m;
  pi = 2.0* acos(0.0);
  ThC = 1.0;                               /* Thermal Conductivity */
  L = N*hx;
  L2 = L*L;
 
  printf("N = %d  L = %8.3f\n",N,L);
  fprintf(output,"N = %d  L = %8.3f\n",N,L);
  printf("hx = %6.4f  steps = %12.0f\n\n",hx,steps);
  fprintf(output,"hx = %6.4f  steps = %12.0f\n\n",hx,steps);
 
/* Initial Conditions for 2-D PDE */
 
  printf("Initial Conditions for 2-D PDE - Heat Conduction Equation:\n\n");
  fprintf(output,"Initial Conditions for 2-D PDE - Heat Conduction Equation:\n\n");
 
  for (i=0; i<=N ; i++)
    {
          x = i*hx;
          u[i] = exp(-m2*pi*pi*ThC*t/L2)*sin(m*pi*x/L);
          printf("i = %3d  x = %6.4f  u[i] = %12.10f\n",i,x,u[i]);
          fprintf(output,"i = %3d  x = %6.4f  u[i] = %12.10f\n",i,x,u[i]);
    }
 
  for (j=1; j*del_t<=MAX ;j++) /* time loop */
     {
       t=j*del_t;
       runge4(t, hx, ThC, u, del_t);
 
/*     INTERMEDIATE print  */
 
       if ( (j % PFREQ) == 0 )
         {
           printf("\nstep = %12d  time = %9.6f\n",j,t);
           fprintf(output,"\nstep = %12d  time = %9.6f\n",j,t);
           for (i=0; i<=N ; i++)
             {
                   x = i*hx;
                   exactf = exp(-m2*pi*pi*ThC*t/L2)*sin(m*pi*x/L);
 
                   if ( i == 0 || i == N )
                     {
                       printf("i = %3d  x = %6.4f  u[i] = %12.10f   ES = %12.10f\n",i,x,u[i],exactf);
                       fprintf(output,"i = %3d  x = %6.4f  u[i] = %12.10f   ES = %12.10f\n",i,x,u[i],exactf);
                     }
                   else
                     {
                       perc_diff = (exactf-u[i])*100.0/exactf;
                      printf("i = %3d  x = %6.4f  u[i] = %12.10f   ES = %12.10f  %12.10f\n",i,x,u[i],exactf,perc_diff);
                      fprintf(output,"i = %3d  x = %6.4f  u[i] = %12.10f   ES = %12.10f  %12.10f\n",i,x,u[i],exactf,perc_diff);
                     }
             }
           }
     }
 
/*    FINAL print  */
 
           printf("\nFINAL STEP:\n");
           fprintf(output,"\nFINAL STEP:\n");
           printf("\nstep = %12.0f  time = %9.6f\n\n",steps,t);
           fprintf(output,"\nstep = %12.0f  time = %9.6f\n\n",steps,t);
           for (i=0; i<=N ; i++)
             {
               x = i*hx;
               exactf = exp(-m2*pi*pi*ThC*t/L2)*sin(m*pi*x/L);
               if ( i == 0 || i == N )
                 {
                   printf("i = %3d  x = %6.4f  u[i] = %12.10f   ES = %12.10f\n",i,x,u[i],exactf);
                   fprintf(output,"i = %3d  x = %6.4f  u[i] = %12.10f   ES = %12.10f\n",i,x,u[i],exactf);
                 }
               else
                 {
                   perc_diff = (exactf-u[i])*100.0/exactf;
                   printf("i = %3d  x = %6.4f  u[i] = %12.10f   ES = %12.10f  %12.10f\n",i,x,u[i],exactf,perc_diff);
                   fprintf(output,"i = %3d  x = %6.4f  u[i] = %12.10f   ES = %12.10f  %12.10f\n",i,x,u[i],exactf,perc_diff);
                 }
 
             }
 
  printf("\nNormal End of Job\n\n - See output in file:\n\n       output6.dat002\n\n");
 
  sec_finish =  time(NULL);
 
  mytime = sec_finish - sec_start;
 
  printf ("Elapsed time: %ld \n", mytime);
  fprintf (output,"Elapsed time: %ld \n", mytime);
 
  fclose(output);
}
 
void runge4(double t, double hx, double ThC, double u[], double del_t)
{
  double ht=del_t/2.0,             /* the midpoint */
  t1[N1], t2[N1], t3[N1],          /* temporary storage arrays */
  k1[N1], k2[N1], k3[N1], k4[N1];  /* for Runge-Kutta */
  int i;
 
  u[0]=0.0;
  u[N]=0.0;
 
  for (i=1;i<N;i++) t1[i]=u[i]+0.5*(k1[i]=del_t*f(t, hx, ThC, u, i));
  for (i=1;i<N;i++) t2[i]=u[i]+0.5*(k2[i]=del_t*f(t+ht, hx, ThC, t1, i));
  for (i=1;i<N;i++) t3[i]=u[i]+ (k3[i]=del_t*f(t+ht, hx, ThC, t2, i));
  for (i=1;i<N;i++) k4[i]= del_t*f(t+del_t, hx, ThC, t3, i);
 
  for (i=1;i<N;i++) u[i]+=(k1[i]+2*k2[i]+2*k3[i]+k4[i])/6.0;
 
}
 
double f(double t, double hx, double ThC, double u[], int i)
  {
    double hx2, x;
 
    hx2 = hx*hx;
    x = i*hx;
    return(ThC*(u[i+1]-2*u[i]+u[i-1])/hx2 );
  }
Results for N = 10
Starting 2-D PDE Runge-Kutta ...
 
N = 10  L =    1.000
hx = 0.1000  steps =      5000000
 
Initial Conditions for 2-D PDE - Heat Conduction Equation:
 
i =   0  x = 0.0000  u[i] = 0.0000000000
i =   1  x = 0.1000  u[i] = 0.3090169944
i =   2  x = 0.2000  u[i] = 0.5877852523
i =   3  x = 0.3000  u[i] = 0.8090169944
i =   4  x = 0.4000  u[i] = 0.9510565163
i =   5  x = 0.5000  u[i] = 1.0000000000
i =   6  x = 0.6000  u[i] = 0.9510565163
i =   7  x = 0.7000  u[i] = 0.8090169944
i =   8  x = 0.8000  u[i] = 0.5877852523
i =   9  x = 0.9000  u[i] = 0.3090169944
i =  10  x = 1.0000  u[i] = 0.0000000000
 
step =       500000  time =  0.500000
i =   0  x = 0.0000  u[i] = 0.0000000000   ES = 0.0000000000
i =   1  x = 0.1000  u[i] = 0.0023141626   ES = 0.0022224142  -4.1283232608
i =   2  x = 0.2000  u[i] = 0.0044017989   ES = 0.0042272830  -4.1283232608
i =   3  x = 0.3000  u[i] = 0.0060585564   ES = 0.0058183559  -4.1283232608
i =   4  x = 0.4000  u[i] = 0.0071222602   ES = 0.0068398875  -4.1283232608
i =   5  x = 0.5000  u[i] = 0.0074887875   ES = 0.0071918834  -4.1283232608
i =   6  x = 0.6000  u[i] = 0.0071222602   ES = 0.0068398875  -4.1283232608
i =   7  x = 0.7000  u[i] = 0.0060585564   ES = 0.0058183559  -4.1283232608
i =   8  x = 0.8000  u[i] = 0.0044017989   ES = 0.0042272830  -4.1283232608
i =   9  x = 0.9000  u[i] = 0.0023141626   ES = 0.0022224142  -4.1283232608
i =  10  x = 1.0000  u[i] = 0.0000000000   ES = 0.0000000000
 
step =      1000000  time =  1.000000
i =   0  x = 0.0000  u[i] = 0.0000000000   ES = 0.0000000000
i =   1  x = 0.1000  u[i] = 0.0000173303   ES = 0.0000159833  -8.4270770511
i =   2  x = 0.2000  u[i] = 0.0000329641   ES = 0.0000304021  -8.4270770511
i =   3  x = 0.3000  u[i] = 0.0000453712   ES = 0.0000418449  -8.4270770511
i =   4  x = 0.4000  u[i] = 0.0000533371   ES = 0.0000491917  -8.4270770511
i =   5  x = 0.5000  u[i] = 0.0000560819   ES = 0.0000517232  -8.4270770511
i =   6  x = 0.6000  u[i] = 0.0000533371   ES = 0.0000491917  -8.4270770511
i =   7  x = 0.7000  u[i] = 0.0000453712   ES = 0.0000418449  -8.4270770511
i =   8  x = 0.8000  u[i] = 0.0000329641   ES = 0.0000304021  -8.4270770511
i =   9  x = 0.9000  u[i] = 0.0000173303   ES = 0.0000159833  -8.4270770511
i =  10  x = 1.0000  u[i] = 0.0000000000   ES = 0.0000000000
 
step =      1500000  time =  1.500000
i =   0  x = 0.0000  u[i] = 0.0000000000   ES = 0.0000000000
i =   1  x = 0.1000  u[i] = 0.0000001298   ES = 0.0000001150  -12.9032972941
i =   2  x = 0.2000  u[i] = 0.0000002469   ES = 0.0000002186  -12.9032972941
i =   3  x = 0.3000  u[i] = 0.0000003398   ES = 0.0000003009  -12.9032972941
i =   4  x = 0.4000  u[i] = 0.0000003994   ES = 0.0000003538  -12.9032972941
i =   5  x = 0.5000  u[i] = 0.0000004200   ES = 0.0000003720  -12.9032972941
i =   6  x = 0.6000  u[i] = 0.0000003994   ES = 0.0000003538  -12.9032972941
i =   7  x = 0.7000  u[i] = 0.0000003398   ES = 0.0000003009  -12.9032972941
i =   8  x = 0.8000  u[i] = 0.0000002469   ES = 0.0000002186  -12.9032972941
i =   9  x = 0.9000  u[i] = 0.0000001298   ES = 0.0000001150  -12.9032972941
i =  10  x = 1.0000  u[i] = 0.0000000000   ES = 0.0000000000
 
step =      2000000  time =  2.000000
i =   0  x = 0.0000  u[i] = 0.0000000000   ES = 0.0000000000
i =   1  x = 0.1000  u[i] = 0.0000000010   ES = 0.0000000008  -17.5643103785
i =   2  x = 0.2000  u[i] = 0.0000000018   ES = 0.0000000016  -17.5643103785
i =   3  x = 0.3000  u[i] = 0.0000000025   ES = 0.0000000022  -17.5643103785
i =   4  x = 0.4000  u[i] = 0.0000000030   ES = 0.0000000025  -17.5643103785
i =   5  x = 0.5000  u[i] = 0.0000000031   ES = 0.0000000027  -17.5643103785
i =   6  x = 0.6000  u[i] = 0.0000000030   ES = 0.0000000025  -17.5643103785
i =   7  x = 0.7000  u[i] = 0.0000000025   ES = 0.0000000022  -17.5643103785
i =   8  x = 0.8000  u[i] = 0.0000000018   ES = 0.0000000016  -17.5643103785
i =   9  x = 0.9000  u[i] = 0.0000000010   ES = 0.0000000008  -17.5643103785
i =  10  x = 1.0000  u[i] = 0.0000000000   ES = 0.0000000000
 
...
Results for N = 100
Starting 2-D PDE Runge-Kutta ...
 
N = 100  L =    1.000
hx = 0.0100  steps =      5000000
 
Initial Conditions for 2-D PDE - Heat Conduction Equation:
 
i =   0  x = 0.0000  u[i] = 0.0000000000
i =   1  x = 0.0100  u[i] = 0.0314107591
i =   2  x = 0.0200  u[i] = 0.0627905195
i =   3  x = 0.0300  u[i] = 0.0941083133
i =   4  x = 0.0400  u[i] = 0.1253332336
i =   5  x = 0.0500  u[i] = 0.1564344650
i =   6  x = 0.0600  u[i] = 0.1873813146
i =   7  x = 0.0700  u[i] = 0.2181432414
i =   8  x = 0.0800  u[i] = 0.2486898872
i =   9  x = 0.0900  u[i] = 0.2789911060
i =  10  x = 0.1000  u[i] = 0.3090169944
i =  11  x = 0.1100  u[i] = 0.3387379202
i =  12  x = 0.1200  u[i] = 0.3681245527
i =  13  x = 0.1300  u[i] = 0.3971478906
...
i =  98  x = 0.9800  u[i] = 0.0627905195
i =  99  x = 0.9900  u[i] = 0.0314107591
i = 100  x = 1.0000  u[i] = 0.0000000000
 
step =       500000  time =  0.500000
i =   0  x = 0.0000  u[i] = 0.0000000000   ES = 0.0000000000
i =   1  x = 0.0100  u[i] = 0.0002259942   ES = 0.0002259025  -0.0405940232
i =   2  x = 0.0200  u[i] = 0.0004517654   ES = 0.0004515821  -0.0405940232
i =   3  x = 0.0300  u[i] = 0.0006770908   ES = 0.0006768160  -0.0405940232
i =   4  x = 0.0400  u[i] = 0.0009017479   ES = 0.0009013820  -0.0405940232
i =   5  x = 0.0500  u[i] = 0.0011255151   ES = 0.0011250584  -0.0405940232
i =   6  x = 0.0600  u[i] = 0.0013481716   ES = 0.0013476246  -0.0405940232
i =   7  x = 0.0700  u[i] = 0.0015694976   ES = 0.0015688607  -0.0405940232
i =   8  x = 0.0800  u[i] = 0.0017892747   ES = 0.0017885487  -0.0405940232
i =   9  x = 0.0900  u[i] = 0.0020072860   ES = 0.0020064715  -0.0405940232
i =  10  x = 0.1000  u[i] = 0.0022233163   ES = 0.0022224142  -0.0405940232
i =  11  x = 0.1100  u[i] = 0.0024371525   ES = 0.0024361636  -0.0405940232
i =  12  x = 0.1200  u[i] = 0.0026485836   ES = 0.0026475088  -0.0405940232
i =  13  x = 0.1300  u[i] = 0.0028574008   ES = 0.0028562413  -0.0405940232
...
i =  98  x = 0.9800  u[i] = 0.0004517654   ES = 0.0004515821  -0.0405940232
i =  99  x = 0.9900  u[i] = 0.0002259942   ES = 0.0002259025  -0.0405940232
i = 100  x = 1.0000  u[i] = 0.0000000000   ES = 0.0000000000
 
step =      1000000  time =  1.000000
i =   0  x = 0.0000  u[i] = 0.0000000000   ES = 0.0000000000
i =   1  x = 0.0100  u[i] = 0.0000016260   ES = 0.0000016247  -0.0812045251
i =   2  x = 0.0200  u[i] = 0.0000032504   ES = 0.0000032477  -0.0812045251
i =   3  x = 0.0300  u[i] = 0.0000048715   ES = 0.0000048676  -0.0812045251
i =   4  x = 0.0400  u[i] = 0.0000064879   ES = 0.0000064826  -0.0812045251
i =   5  x = 0.0500  u[i] = 0.0000080979   ES = 0.0000080913  -0.0812045251
i =   6  x = 0.0600  u[i] = 0.0000096998   ES = 0.0000096920  -0.0812045251
i =   7  x = 0.0700  u[i] = 0.0000112922   ES = 0.0000112831  -0.0812045251
i =   8  x = 0.0800  u[i] = 0.0000128735   ES = 0.0000128630  -0.0812045251
i =   9  x = 0.0900  u[i] = 0.0000144420   ES = 0.0000144303  -0.0812045251
i =  10  x = 0.1000  u[i] = 0.0000159963   ES = 0.0000159833  -0.0812045251
i =  11  x = 0.1100  u[i] = 0.0000175348   ES = 0.0000175206  -0.0812045251
i =  12  x = 0.1200  u[i] = 0.0000190560   ES = 0.0000190406  -0.0812045251
i =  13  x = 0.1300  u[i] = 0.0000205584   ES = 0.0000205418  -0.0812045251
...
i =  98  x = 0.9800  u[i] = 0.0000000234   ES = 0.0000000234  -0.1218315124
i =  99  x = 0.9900  u[i] = 0.0000000117   ES = 0.0000000117  -0.1218315124
i = 100  x = 1.0000  u[i] = 0.0000000000   ES = 0.0000000000
 
step =      2000000  time =  2.000000
i =   0  x = 0.0000  u[i] = 0.0000000000   ES = 0.0000000000
i =   1  x = 0.0100  u[i] = 0.0000000001   ES = 0.0000000001  -0.1624749919
i =   2  x = 0.0200  u[i] = 0.0000000002   ES = 0.0000000002  -0.1624749919
i =   3  x = 0.0300  u[i] = 0.0000000003   ES = 0.0000000003  -0.1624749919
i =   4  x = 0.0400  u[i] = 0.0000000003   ES = 0.0000000003  -0.1624749919
i =   5  x = 0.0500  u[i] = 0.0000000004   ES = 0.0000000004  -0.1624749919
i =   6  x = 0.0600  u[i] = 0.0000000005   ES = 0.0000000005  -0.1624749919
i =   7  x = 0.0700  u[i] = 0.0000000006   ES = 0.0000000006  -0.1624749919
i =   8  x = 0.0800  u[i] = 0.0000000007   ES = 0.0000000007  -0.1624749919
i =   9  x = 0.0900  u[i] = 0.0000000007   ES = 0.0000000007  -0.1624749919
i =  10  x = 0.1000  u[i] = 0.0000000008   ES = 0.0000000008  -0.1624749919
i =  11  x = 0.1100  u[i] = 0.0000000009   ES = 0.0000000009  -0.1624749919
i =  12  x = 0.1200  u[i] = 0.0000000010   ES = 0.0000000010  -0.1624749919
i =  13  x = 0.1300  u[i] = 0.0000000011   ES = 0.0000000011  -0.1624749919
...

Since hx = 1.0/N and the truncation error of the second-order derivative in the x-direction is of the order hx it follows that the results for N=10 (hx=0.1) have a much higher truncation error as the results for N=100 (hx=0.01). This can be seen by observing the last column of the above results.