Parabolic PDEs in julia

Question

I am trying to solve a parabolic partial differential equation numerically using Julia, but I cannot find any accessible documentation that can help.

Here is an example: t, x are 1 dimensional real. I want to solve for u(t,x)=[u1(t,x) u2(t,x)]; u satisfies the PDE

du1/dt = d^2u1/dx^2 + a11(x,u) du1/dx + a12(x,u) du2/dx + c1(x,u)

du2/dt = d^2u2/dx^2 + a21(x,u) du1/dx + a22(x,u) du2/dx + c2(x,u)

Is it possible to do so in Julia? This type of problem is solvable using pdepe in Matlab.

Answer 1

Right now, we don't have "full-stop" PDE solvers, i.e. solvers where you put in the PDE and go. However, PDEs are solved by discretizing to ODEs, so the way a full-stop PDE solver would be written for this is as follows. Most of this is discussed in more depth in this blog post BTW.

Take your PDE. Now discretize the operators by dx. The second order finite difference discretization of the LaPlacian is the stencil u[i-1] - 2 u[i] + u[i+1] applied to our state. Then of course you have to take into account your boundary conditions when you hit the end. Usually a nice way to write this is as a matrix, so:

const Mx = Tridiagonal([1.0 for i in 1:N-1],[-2.0 for i in 1:N],[1.0 for i in 1:N-1])
# Do the reflections, different for x and y operators
Mx[2,1] = 2.0
Mx[end-1,end] = 2.0

has Mx*u/dx^2 perform the discretized LaPlacian.

The first order derivative terms are done similarly, but in this case it's common to use the upwinding scheme. You can take your du1/dx term and replace it by the kernel

a[i]*(u[i]-u[i-1])/dx

when a is positive, or

a[i]*(u[i]-u[i+1])/dx

when a is negative. Then incorporate boundary conditions of course. Then you just write your reactions as c1(x[i],u[i]). Together this will look like (in the non-matrix form:

function f(t,u,du)
    u1 = @view u[:,1]
    u2 = @view u[:,1]
    du1 = @view du[:,1]
    du2 = @view du[:,2]
    for i in 2:length(u)-1
        du1[i] = (u1[i-1] - 2u1[i] + u1[i+1])/dx^2 +
                a11(x[i],u1[i])*(u1[i]-u1[i-1])/dx +
                a12(x[i],u1[i])*(u1[i]-u1[i-1])/dx +
                c1(x1[i],u1[i])

        du2[i] = (u2[i-1] - 2u2[i] + u2[i+1])/dx^2 +
                a11(x[i],u2[i])*(u2[i]-u2[i-1])/dx +
                a12(x[i],u2[i])*(u2[i]-u2[i-1])/dx +
                c1(x1[i],u2[i])
    end

end

Notice that I did not do the ends because I don't know what boundary conditions you want. If it's Dirichlet with a zero constant, then you just write it at the end points but drop off the values that's off of space. And here x[i] = x0 + dx*i.

Now you have a set of ODEs where u[i,j] = u_j(x_i). Thus you discretize your initial condition into u0[i,j] and setup an ODE problem:

using DifferentialEquations
prob = ODEProblem(f,u0,tspan)

For this, see the DiffEq documentation, specifically the ODE tutorial. Now you just solve the discretized ODE representation of the PDE. For these kinds of equations, as mentioned in the blog post, the Sundials.jl CVODE_BDF method with the GMRES Krylov linear solver is a good choice, so we do:

sol = solve(prob,CVODE_BDF(linear_solver=:GMRES))

That gives back a continuous solution where sol(t)[i,j] is the numerical approximation to u_j(t,x_i). Of course, lower dx is more accurate, and you should adjust the tolerances of the ODE solver as you see fit.

We will have things to do this automatically for PDEs in the near future (any derivatives to any order), but it's currently a work-in-progress so for now manual discretization will have to do (and this is taught in any numerical methods course, so it's not too bad!). Hope this helps. If you need any more help, please check out our chat channel as most people in there will have experience doing this kind of discretization.