EMAT30008: Scientific Computing

Week 20: The method of lines and explicit Euler

Overview

This week we enter the realm of numerically solving partial differential equations (PDEs). We will focus on one of the most important PDEs of them all: the diffusion equation.

By using finite differences to discretise space, the diffusion equation can be converted into a system of ODEs. This provides an opportunity to link the finite-difference code you developed in Week 19 with the ODE time-steppers you developed in Week 14. For example, the euler_step and solve_to functions you created in Week 14 can be used to implement the explicit Euler method described in the Week 20 video.

Supplementary material

Use the links below to find additional notes on

Exercise

The goal of this week is to create a PDE solver that is capable of computing numerical solutions to diffusion equations of the form

ut=D2ux2+q(x,t,u;μ). \frac{\partial u}{\partial t} = D \frac{\partial^2 u}{\partial x^2} + q(x, t, u; \mu).

The domain of the problem is given by axba \leq x \leq b. In this equation, D>0D > 0 is a parameter. The term q(x,t,u;μ)q(x, t, u; \mu) represents a function that depends on the spatial coordinate xx, time tt, the solution u(x,t)u(x, t), and a parameter μ\mu.

The user should be able to specify which type of boundary conditions to impose (Dirichlet, Neumann, or Robin) as well as the method used to time step the equation (e.g. explicit Euler, RK, solve_ivp).

If qq is independent of time, then the steady-state solutions of (1) are given by

Dd2udx2+q(x,u;μ)=0, D \frac{\mathrm{d}^2 u}{\mathrm{d} x^2} + q(x, u; \mu) = 0,

which is exactly the equation you were solving in Week 19. Therefore, you can use your numerical solutions of (1) to determine the stability of the steady-state solutions computed from (2).

Steps

  1. (Essential) Use the explicit Euler method to solve the linear diffusion equation without a source term

    ut=D2ux2. \frac{\partial u}{\partial t} = D\frac{\partial^2 u}{\partial x^2}.

    Set the boundary conditions to u(a,t)=0u(a,t) = 0 and u(b,t)=0u(b,t) = 0 and the initial condition to

    u(x,0)=sin(π(xa)ba). u(x,0) = \sin\left(\frac{\pi (x - a)}{b - a}\right).

    The exact solution to the problem is given by

    u(x,t)=exp(Dπ2t(ba)2)sin(π(xa)ba). u(x, t) = \exp\left(-\frac{D \pi^2 t}{(b-a)^2}\right)\sin\left(\frac{\pi (x - a)}{b - a}\right).

    which you can use to test your code.

  2. (Essential) Use the 4th-order Runge-Kutta method from Week 14 to solve the diffusion equation above.

  3. Extend your code so that it can account for non-homogeneous Dirichlet boundary conditions and a source term in the diffusion equation. That is, your code should be able to solve

    ut=D2ux2+q(x,t,u;μ), \frac{\partial u}{\partial t} = D \frac{\partial^2 u}{\partial x^2} + q(x, t, u; \mu), \\

    with boundary conditions u(a,t)=αu(a, t) = \alpha, u(b,t)=βu(b,t) = \beta and initial condition u(x,0)=f(x)u(x,0) = f(x). Use the supplementary notes to find exact solutions that you can use to test your code with.

  4. Use your code to solve the dynamic Bratu problem

    ut=D2ux2+eμu,\begin{aligned} \frac{\partial u}{\partial t} = D \frac{\partial^2 u}{\partial x^2} + e^{\mu u}, \end{aligned}

    with u(0,t)=0u(0, t) = 0, u(1,t)=0u(1,t) = 0 and initial condition u(x,0)=0u(x,0) = 0. Take D=1.0D = 1.0. Compute numerical solutions when μ=2\mu = 2 and compare to the steady-state solution that you computed in the Week 19 exercises.

  5. Update your code so that it can also account for Dirichlet, Neumann, and Robin boundary conditions with time-dependent coefficients.