Differential Equations
Mathwrist Code Example Series
Abstract
This document gives some code examples on how to use differential equation solvers in Mathwrist’s C++ Numerical Programming Library (NPL).
1 HIRES Benchmark Problem
The HIRES benchmark test problem refers to “High Irradiance RESponse”. It originates from plant physiology and describes how light is involved in morphogenesis. It is a stiff nonlinear ODE system of 8 equations,
The solution at
Since the HIRES problem is nonlinear, at line 2 we create a struct Hires derived from NPL’s general ODE interface class ODE. Lines 4 to 6 implement base class virtual functions to specify the ODE’s size, start and end time. We supply an implementation of the y_prime virtual function from line 9 to 20 to compute
At line 45, we create an IDeC solver to use
2 An Example of Linear ODE
This is an example of linear ODE to demonstrate NPL’s implicit IDeC method. The problem is defined as the following:
It has an analytic solution
At line 2, we now derive from NPL’s linear ODE interface ODE_Linear and define the behavior of this ODE via a set of virtual functions. For large linear ODE systems,
From line 15 to 22, we supply an implementation to compute
From line 37 to 49, we solve an linear system
Lastly, line 52 to 60 computes
3 PDE - Heat Equations
In this example we have 2 unknown functions
They differ from each other on the following initial values
, and have the analytic solutions
We will create different solvers to solve
At line 4, we create a derived struct Heat under the PDEParabolic1D interface. From line 9 to 18, we implement the coefficient pure virtual function by filling the
Struct F0 from line 24 to 44 implements the initial value function interface and fill

Between line 98 and 102, we implement the Dirichlet boundary condition. Line 113 and 114 create a PDE problem specification object, where we specify that there are 2 problems to be solved together. Code block 117 to 132 uses Crank-Nicolson method. The numerical solutions are stored in 2 separate rows of the vx matrix parameter on return. We then check the solution error from line 127 to 131.
Code blocks from line 135 to 163 use the method of lines with the 2nd order finite difference scheme and the 4th order finite difference scheme. In each case, we use 3 iterations of error correction in implicit IDeC method. The 2nd order method obtains similar result as the Crank-Nicolson method. The 4th order method on the other hand produces 1.e4 times more accurate result. Figure 2 is a plot of the final numerical solutions.

Finally code block between line 165 and 206 illustrates the spectral collocation method, where we choose Chebyshev approximation polynomial of degree 11 and 3 error correction iterations. Note that the output matrix C at line 170 consists of two rows of Chebyshev basis coefficients. We need construct 2 Chebyshev approximation polynomials with these coefficients at line 175 and 180 in order to compute the actual numerical solutions. Code block 182 to 197 uses these two approximation polynomials to evaluate the maximum numerical error over 100 discrete points. As we can see from Figure 1, the initial values are very oscillating. Chebyshev collocation method is not a good choice for this situation. Not surprisingly, Chebyshev collocation method produces large error for not very smooth functions.