Author: Douglas Wilhelm Harder
Note that many algorithms are now implemented on repl.it, so this page should be considered somewhat archaic.
In the source directory are some of the algorithms implemented in the course, especially related to initial-value problem (IVP) solvers.
Associated with many of the algorithms is a header file (.h) and a template include file (.tpp). In addition to these are either example or test programs (.cpp) with the same name. The output of the example or test program is in the corresponding file with the extension .output.
All files are in the source directory.
As one example where the forcing function is discontinuous, suppose we have the differential equation $y'(t) + y(t) = {\rm u}(t - 1) - {\rm u}(t - 2)$ where ${\rm u}(t)$ is the Heaviside unit-step function. The solution is continuous but not differentiable at $t = 1$ and $t = 2$. Consequently, adaptive techniques will necessarily require smaller and smaller step sizes, as no step size will be sufficiently small to deal with the discontinuity in the derivative. It is at these points that the minimum step size is met.
Figure 1. An approximation to the solution evaluated at $t = 0.01k$ for $k = 0, ..., 300$.
Subtracting off the error, we see that where the approximation points are, the approximation is very close; however, where the interpolating polynomial is far from the points of approximation, the error increases. Note that near the discontinuity, the approximation points get significantly closer.
Figure 2. The error of each value of the approximation at the points $t = 0.01k$ for $k = 0, ..., 300$.
You can see the source code for this in discont.cpp.
As an example of a second-order initial-value problem converted to a system of two first-order initial-value problems, consider the differential equation $y''(t) + y(t) = 0$, the solution to which is sinusoid. Specifically, if the initial conditions are that $y(0) = y'(0) = 1$, the solution is $y(t) = \sin(t) + \cos(t)$. A plot of approximating this solution at points separated by a distance of $h = 0.01$ on the interval $[0, \pi]$ is shown in Figure 3.
Figure 3. An approximation to the solution evaluated at $t = 0.01k$ for $k = 0, ..., 628$.
We can also plot the error, and you will see that at the approximation points, the error is minimized, but between these points, the error is cubic.
Figure 4. The errors of the approximations of the solution to the sinusoidal initial-value problem.
You can see the source code for this in sinusoid.cpp.