Demo: The 1-D heat equation

This code illustrates the numerical solution of the 1-D heat equation using a finite difference scheme.

  • For conciseness, we skip the description of the heading of the main program. This has been presented previously. We assign values of the length of the interval. We also declare a file where the solution is stored

       const double L = 1.;
   ofstream pf("output.dat");

  • We now declare discretization data. nx is the number of grid sub-intervals and h is the grid size. The time step is obtained as a program argument and the final time is set to 1. These two parameters are assigned to global variables in .

       size_t nx = atoi(argv[1]);
   double h = L/double(nx);
   theTimeStep = atof(argv[2]);
   theFinalTime = 1.;

  • We declare the problem matrix, which is tridiagonal (instance of class TrMatrix<double>) and then declare x, f, u and v which are respectively the vector of grid point coordinates, soure term, solution at previous time step and solution at actual one. Then we define grid as a uniform one, and define the initial condition as the function sin(πx):

       Vect<double> x, f(nx+1), u(nx+1), v(nx+1);
   x.setUniform(0.,L,nx+1);
   u.set("sin(pi*x)",x);

  • We next start the time stepping procedure. This is simplified in by using the macro TimeLoop. Note that this macro uses the global variables theTimeStep, theFinalTime and theTime that have obvious meaning. 

       TimeLoop {

  • For each time step, we construct the linear system of equations. First we consider the source vector, then we build the matrix and the right-hand side vector. Note that it is mandatory to assign the actual time value to the vector f in order to properly intepret the algebraic expression assigned to it. This expression depends indeed on the time variable.

          f.setTime(theTime);
      f.set("exp(t)*sin(pi*x)*(1+pi^2)",x);
      for (int i=2; i<=nx; i++) {
         A(i,i  ) =  1./theTimeStep + 2./(h*h);
         A(i,i+1) = -1./(h*h);
         A(i,i-1) = -1./(h*h);
         v(i) = u(i)/theTimeStep + f(i);
      }

  • We impose Dirichlet boundary conditions:

          A(1,1) = 1./theTimeStep;
      A(1,2) = 0.;
      v(1) = 0.;
      A(nx+1,nx+1) = 1./theTimeStep;
      A(nx+1,nx) = 0.;
      v(nx+1) = 0.;

  • We solve the linear system and transfer vector v to u: 

          A.solve(v);
      u = v;

  • We can now store the solution at current time step. 

          for (size_t i=1; i<=nx+1; i++)
         pf << x(i) << "   " << u(i) << endl;
         pf << endl;
      }