Solution of the heat equation

Directly to the input form

• The heat equation is the most prominent example of a parabolic PDE. Here we are considering the case of one space dimension only. The equation (with constant coeffcients) has the form
  (d/dt)u = a (d/dx)**2 u + f(x,t)
  where x in ]0,1[ and t > 0.
• Parabolic partial differential equations describe diffusion processes.
• In order to get a well defined problem one needs initial values for t=0 and boundary conditions for x=0,1 and t > 0. A classical smooth solution requires additionally compatibility conditions between the initial and the boundary conditions. Here we use Dirichlet data only, that is we prescribe the values of the solution at the named boundaries. Observe that it is impossible to prescribe data for an ''end time''.
• You may think of u as a temperature of an insulated slab (supply or loss of energy takes place at the two endpoints x=0,1 of the slab) with an initial (given) temperature distribution which is heated or cooled at its two ends. The function f describes an internal heat source.
• The initial values u(x,0), x in [0,1] define the initial temperature and the boundary values at x = 0 resp. x = 1 the temperature provided ''outside'' for any given time t.
• We solve the problem by a two stage process, the so called vertical method of lines. In a first step we restrict x to a finite grid
  0=x₀ < x₁ < ... < x_N+1=1,
  considering only v_i(t) = u(x_i,t) and replacing u_xx(x_i,t) by the second order symmetric difference quotient
  (v_i+1(t)-2 v_i(t)+v_i-1(t))/dx² .
  Of course (d/dt)u(x_i,t)= (d/dt)v_i(t). Hence, using the known boundary values we get a coupled system of N ordinary first order differential equations for the v_i(t), i=1,..,N.
  We write this as
  (d/dt)v = F(t,v).
  v and F are now vectors of length N.
• This system of ordinary differential equations is a so called ''stiff'' system: the eigenvalues of its Jacobian matrix range between [-4a/dx²,-a*pi²+O(1/N²)]. Hence the use of an explicit time integrator is not advisable (although possible, and indeed, in space dimension 2 or 3 this can be advantageous using special integrators). Here we use an implicit linear multistep integrator, the so called ''BDF'' (backward differentiation) or ''Gear'' method.
• The BDF method of step number k (and convergence order k in the time step) obtains if one interpolates the data (t_j+1-k,v_j+1-k),... ,(t_j+1,v_j+1) by a polynomial of degree <=k (in the variable t) , differentiates this one and evaluates the derivative at t=t_j+1 (this is the new time where the solution is sought) and equates this with with F(t_j+1,v_j+1). This method is asymptotically stable and A(α)-stable only for 1<=k<=5. k=1,2 give methods which are A-stable and L-stable (α=π/2) where k=1 is ''Euler implicit'' whereas the higher order methods have a more and more shrinking α.
• The code which we are using here (vode from netlib/ode) works adaptively in timestep and order.

Input

• There are two predefined cases with known exact solution and you also have the possibility to specify a problem of your own.
• In case of a problem of your own you must specify
  1. The function f(x,t). This input must follow the rules of FORTRAN.
  2. The boundary functions φ₀(t) and φ₁(t) and the initial function u₀(x) (also in FORTRAN).
• In all three cases you then specify the thermal conduction coefficient a. Important : a > 0 (otherwise this would correspond to a process ''backwards in time'' which is illposed)!
• the number N of interior grid points in x-direction. Important : 1<=N<=200 !
• The endtime of the process tend. Important: tend > 0 !
• The desired relative precision tol for the time integrator. Attention: Since N is chosen beforehand, a very small tol doesn't mean that you get a precise solution of u but only of the vector v!
• You can limit the order of the time integrator ord, hence of the BDF-method. Of course : ord from {1,2,3,4,5} ! ord =1 means you are using Euler implicit only.
• The number grid of time lines you want to see on the plot.
  0 <= grid <= 200 ! grid = 0 means that you get all internally computed time lines (which can be quite large) whereas grid > 0 results in that number of t-equidistant grid lines until your specified endtime. The number of internally computed grid lines may be much larger.
• Since you get a 3D projected plot, you must also specify your viewpoint via two rotation angles for the x and the u axes.

Output

• You get a protocol of your chosen input parameters and a short statistics of the internal effort:
• This includes:
  1. The last time step used.
  2. The last used order of the BDF method.
  3. The total number of time steps.
  4. The number of Newton steps for the solution of the systems for each time step.
     This is the number of all tried time steps, since in our case the implicit equation is linear.
  5. The number of evaluations of the right hand side of the ordinary differential equation.
  6. The number of evaluations of the Jacobians of the ordinary differential equation which here reads I-δ_t*a/(δ_x)²*tridiag(1,-2,1)
• In the two predefined cases you get a plot of the discretization error and in case of a problem of your own of the computed solution u.

Questions ?!

• How behaves the discretization error in dependency of n and tol ?
• What takes place if one limits the maximal order of the time integrator?
• What must be observed in specifying the initial and boundary functions in order to get a smooth solution?

To the input form