ApproxFun and Burgers' equation

I am trying to use ApproxFun to solve Burgers’ equation, u_t + uu_x = \mu u_{xx}, on the interval -L \leq x \leq L with initial condition u(x, 0) = f(x) = 1/(1+x^2) and boundary conditions u(\pm L, t) = 0. I was trying to modify some of the code here Burgers Pseudospectral Methods Work-Precision Diagrams for my case, but it doesn’t explain how to incorporate boundary conditions / what to do in the non-periodic case where I’d assume ChebyshevInterval would be used instead.

Does anyone know of an example for my case, or have any guidance for incorporating such boundary conditions / setting up the differential operators required? Thank you. I’m trying to replicate the following Chebfun code from MATLAB:

M = 250;
N = 100;
mu = 0.05;
L = 40;
T = 8/sqrt(27);
t = linspace(0, T, N);
a = -L;
b = L;
xvals = linspace(a, b, M);
d = [a b];
pdefun = @(t, x, u) mu * diff(u, 2) - diff(u.^2/2);
bc.left = @(u) u;                     
bc.right = @(u) u;                     
x = chebfun(@(x) x, d);
f = 1/(1+x^2);
opts = pdeset('Eps', 1e-4);
[t, u_num] = pde15s(pdefun, t, f, bc, opts);

into ApproxFun.

Boundary conditions make time-stepping more complicated as it becomes an DAE; ie a mass matrix but the mass matrix singular. With some simple eliminations you can incorporate the BCs to get a non-singular mass matrix which can then be inverted to reduce to an ODE (though one with a dense matrix)

There used to be a package to help:


Thanks! In that package, would the “better ways” in the deprecation notice,

This repository has been deprecated. There are better ways to do this.

be the “simple eliminations” you mention? I’m not so familiar with this currently so I’m not sure what eliminations you are referring to / what I’d be eliminating. I’ve also edited in my original question the actual Chebfun code I’m trying to replicate in ApproxFun.

Basically, it would be better to do this by canonicalizing the input and using a canonical form rather than the type at all points of the internals.

For the boundary conditions, choose a space that implicitly supports the BC, or make it a DAE where the algebraic equation is what’s required from the basis set for the BC(s) to be satisfied.

There’s also the Chebyshev tau method, where you replace the last couple (highest-order) discretized differential equations with equations for the boundary conditions, resulting in a nonsingular LHS matrix. With the right transformation of the Chebyshev modes, the u_{xx} term plus BCs produces a banded tridiagonal diagonal system that can be solved in O(n) time.

This is discussed in detail in Canuto, Hussaini, Quarteroni, and Zang’s “Spectral Methods in Fluid Dynamics” text.


The “simple eliminations” are essentially, to use as the basis set linear combinations of Chebyshev polynomials that satisfy the boundary conditions individually. That’s what Chris is saying: “choose a space that implicitly supports the BC”.

For example, the even Chebyshevs satisfy T_{2n}(\pm1) = 1 and the odd satisfy T_{2n+1}(\pm 1) = \pm 1. So define a new set of basis functions \phi_n(x) = T_{n+2}(x) - T_n(x), so that \phi_n(\pm 1) = 0, and expansions over \phi_n automatically satisfy the BCs.

1 Like

This is actually easy to do in ClassicalOrthogonalPolynomials.jl (Your suggestion is equivalent to the basis (1-x^2) * U_n(x) which corresponds to JacobiWeight(1,1) .* ChebyshevU()).

I’ll try to cook up an example at some points