Numerical methods - The Elara Handbook

Topics in numerical methods¶

Numerical methods are focused on ways to solve problems that, in general, cannot be solved exactly.

We will describe both how to write these solvers from scratch, as well as how to write. For each, we will provide code examples using vanilla Python (with NumPy) and other libraries as needed.

ODEs and initial-value problems¶

In our discussion on differential equations, we mentioned that ordinary differential equations take the general form:

\dfrac{dy}{dt} = f(t, y)

(1)

(note that $y$ can be a vector, and reduction of order may be necessary to cast a particular ODE in this form). An initial-value problem is a description of a physical scenario in which an ODE is provided along with the initial conditions, that is, $y(0)$ , the value of $y$ at $t = 0$ . Often, we cannot solve initial-value problems exactly, and some initial-value problems are so complicated that even approximate solutions cannot be found. We therefore turn to numerical methods.

The general basis of solving an initial-value problem numerically is to discretize $y(x)$ into a set of line segments. Each segment connects the $n$ th-point $(t_n, y_n)$ to its next point $(t_{n + 1}, y_{n+1})$ . By starting at the known initial condition $(0, y(0))$ , we may then use the fundamental of calculus to calculate the next point:

y_{n+1} = \int_{t_n}^{t_{n+1}}f(t, y_n) d t

(2)

This process is known as numerical integration and is the principle behind solving ODEs. Different ways of computing the integral numerically result in a diverse class of specialized numerical integration schemes, each of which has its own advantages and disadvantages.

forwards Euler, backwards Euler, midpoint, leapfrog, runge-kutta 2 & 4, also some bit about symeptic integrators

PDEs and boundary-value problems¶

The most common methods of solving boundary-value problems numerically are as follows:

Finite difference methods, which approximate derivatives (which are limits of the difference quotient) with finite difference quotients, so a PDE can be reformulated as a system of equations and solved with some linear algebra. Subcategories of this include the method of lines. This category of numerical methods usually only work on a regular and discrete grid, so it is out of the question for our numerical work, which involves very nontrivial geometries for parabolic reflectors with openings and secondary mirrors
Finite element methods as well as the related spectral methods first turn a PDE into an integral equation. The solution is approximated as a linear combination of basis functions (like simple polynomials that can be easily integrated) and their coefficients. The domain is then discretized so the integral can then be evaluated, resulting (again) in a system of equations and solved with some linear algebra.
Finite volume methods involve cleverly rewriting the PDE in the form of a continuity equation, which can be converted into an equivalent integral equation that is then discretized and solved. This has the advantage of preserving conservation laws in the system, and in the future we will probably use these a lot more, but I haven’t tried them yet.

The finite difference method¶

More info here: https://labs.elaraproject.org/theoretical/Guide-to-partial-differential-equations.html#numerically-solving-pdes

The finite element method¶

a) Assume the solution takes the form $u(\mathbf{x}) = \sum_i \varphi_i \phi_i(\mathbf{x})$ where $\phi_i(x)$ is a set of functions that form a complete basis. A popular choice is the Chebyshev polynomials, you can also use Legendre or Lagrange polynomials, but any basis will do. The more terms in the series expansion/approximation, the more accurate your solution will be. b) By autodiff (or just hard-coded derivative rules), substitute $u(x) = \sum_i \varphi_i \phi_i(\mathbf{x})$ into the ordinary/partial differential equation in question where the ordinary/partial differential equation is written in the form $F(u, \nabla u, \nabla^2 u) = 0$ and $F$ is some combination of linear and nonlinear operators. This substitution will result in a system of (possibly nonlinear) equations in the form $\mathbf{G}(\mathbf{x}) = 0$ . c) Solve this system of homogenous equations either with a linear algebra solver or Newton’s method for nonlinear