Boundary value problems: not all explicit solutions are useful

Consider this linear differential equation: {y''(t) + 4y'(t) + 2t y(t) = 7} with boundary conditions {y(0)=1} and {y(1)=0}. Nothing looks particularly scary here. Just one nonconstant coefficient, and it’s a simple one. Entering this problem into Wolfram Alpha produces the following explicit solution:


I am not sure how anyone could use this formula for any purpose.

Let us see what simple linear algebra can do here. The differential equation can be discretized by placing, for example, {4} equally spaces interior grid points on the interval: {t_k = k/5}, {k=1, \dots, 4}. The yet-unknown values of {y} at these points are denoted {y_1,\dots, y_4}. Standard finite-difference formulas provide approximate values of {y'} and {y''}:

{\displaystyle y'(t) \approx \frac{y(t+h) - y(t-h)}{2h}}

{\displaystyle y''(t) \approx \frac{y(t+h) - 2y(t) + y(t-h)}{h^2}}

where {h} is the step size, {1/5} in our case. Stick all this into the equation: we get 4 linear equations, one for each interior point. Namely, at {t_1 = 1/5} it is

{\displaystyle \frac{y_2 - 2y_1 + 1}{(1/5)^2} + 4 \frac{y_2 - 1}{2/5} + 2\cdot \frac15 y_1 = 7 }

(notice how the condition {y(0)=1} is used above), at {t_2 = 2/5} it is

{\displaystyle \frac{y_3 - 2y_2 + y_1}{(1/5)^2} + 4 \frac{y_3 - y_1}{2/5} + 2 \cdot \frac25 y_2 = 7 }

and so on. Clean up this system and put it in matrix form:

{\displaystyle \begin{pmatrix} -49.6 & 35 & 0 & 0 \\ 15 & -49.2 & 35 & 0 \\ 0 & 15 & -48.8 & 35 \\ 0 & 0 & 15 & -48.4 \end{pmatrix} \vec y = \begin{pmatrix} -8 \\ 7 \\ 7 \\ 7 \end{pmatrix} }

This isn’t too hard to solve even with pencil and paper. The solution is

{\displaystyle y = \begin{pmatrix} -0.2859 \\ -0.6337\\ -0.5683\\ -0.3207 \end{pmatrix}}

It can be visualized by plotting 4 points {(t_k, y_k)}:

Discrete solution

Not particularly impressive is it? And why are all these negative y-values in a problem with boundary condition {y(0)=1}? They do not really look like they want to approach {1} at the left end of the interval. But let us go ahead and plot them together with the boundary conditions, using linear interpolation in between:

Linear algebra + linear interpolation

Or better, use cubic spline interpolation, which only adds another step of linear algebra (see Connecting dots naturally) to our computations.

Same points, cubic spline interpolation

This begins to look believable. For comparison, I used a heavier tool: BVP solver from SciPy. Its output is the red curve below.

Cubic spline and BVP solver

Those four points we got from a 4-by-4 system, solvable by hand, pretty much tell the whole story. At any rate, they tell a better story than the explicit solution does.

Graphics made with: SciPy and Matplotlib using Google Colab.

What difference a difference scheme makes

One of the simplest nonlinear PDE is

{\displaystyle u_t = -uu_x}

inviscid Burgers equation in one dimension. The function {u} represents the velocity of a compressible fluid. The equation can model a shock wave: for discontinuous initial data such as

{\displaystyle u(x,0) = \begin{cases} 1, \quad x\le 0 \\ 0, \quad x>0 \end{cases}}

the solution (understood in appropriate weak sense) has the jump from {1} to {0} at all times, but the position of the jump moves to the right. The equation can also demonstrate the creation of shocks from continuous initial data (nonlinear steepening of waves), but I am not going into that here.

Let’s try to solve this PDE numerically. Given the initial {u(x,0)} as a sequence of numbers {\dots 1 1 1 1 1 0 0 0 0 0 \dots} in the top row of a spreadsheet, I use the difference of adjacent entries as approximation to {u_x}, multiply it by {u}, and use the result to calculate the solution at the next moment of time. The time step is set at {\Delta t=0.1}. On the first try, I use the forward difference scheme for {u_x}, which in spreadsheet terms amounts to formulas such as

F2 = F1 - 0.1*F1*(G1-F1)

The result of running this scheme is nothing short of disastrous:

Forward (downwind) difference
Forward (downwind) difference

Let’s try the backward (upwind) difference (“wind direction” is positive, since we work with positive values of velocity). This means using formulas such as

F2 = F1 - 0.1*F1*(F1-E1)

The numerical solution is less ridiculous, but is still quite wrong: it is constant in time.

Backward (downwind) difference
Backward (downwind) difference

Going back to the PDE, rewrite it as

{\displaystyle u_t = -(u^2/2)_x}

This is the conservation form of the PDE, which is more suitable for interpretation of discontinuous solutions. Indeed, it says exactly that the vector field {V(x,t) = (u^2/2,u)} has zero divergence. The property of having zero divergence can be understood without taking derivatives: it means that the flux integral over every closed curve (e.g., over every rectangle) is zero.

Discretizing the conservation form of the PDE leads to finite differences of squares, such as

F2 = F1 - 0.1*(G1^2-F1^2)/2

(conservative downwind method):

Conservative downwind method
Conservative downwind method

Still not good. But the conservative upwind method, with formulas such as

F2 = F1 - 0.1*(F1^2-E1^2)/2

turns out to work:

Conservative upwind method
Conservative upwind method

Here is the plot of the data:

Numerical shock propagation with the conservative upwind method.
Numerical shock propagation with the conservative upwind method.

The shock front looks a bit diffuse, but it sustains itself: note how the value at “4” converges to {1}. Recalling that every row corresponds to time interval of {0.1}, we can estimate the propagation speed to be about {0.5}, which is what the analysis of PDE would predict.