When solving numerical problems like partial differential equations, knowing how accurate your solution is is always important. But how can you know if you don’t have the exact solution for the problem? After all, this is why you resorted to a numerical solution scheme in the first place, right? Well, there is a powerful way how you can get around this problem: the method of manufactured solutions. Here is how it works.

The Method of Manufactured Solutions (MMS) involves creating an artificial, or “manufactured,” solution that serves as a known answer to the equation. By working backward from this predetermined solution, we can accurately verify the numerical solution process’s correctness. The first step is to propose a solution that meets the boundary and initial conditions, and also satisfies all the physical requirements of the problem. This arbitrary solution is then substituted into the equation, yielding a new source term. The equation, coupled with the source term, represents a well-posed problem that can be solved numerically. This method presents the advantage of complete control over the equation’s outcome, thereby offering a comprehensive examination of numerical schemes’ accuracy and stability. The primary purpose of using MMS is not to mirror real-world situations, but to ensure the mathematical and computational soundness of the methods used to solve numerical problems with are actually too complex to solve analytically.

Let’s see how this works by a simple example.

Consider the 2D steady-state heat conduction equation:

$$ \frac{\partial^2 T}{\partial x^2}+\frac{\partial^2 T}{\partial y^2}=\rho $$

where T is temperature, and x and y are spatial dimensions, and ρ is some heat source or sink. And suppose that have implemented a numerical solution scheme for that. We will manufacture a solution that will work as a benchmark for our implementation, how accurate it is.

We first manufacture an exact solution that satisfies the problem’s boundary conditions. For simplicity, let’s choose a polynomial:

$$ T_{\text {exact }}=x^2+y^2 $$

Substituting this solution into the heat conduction equation, we get a very simple source term:

$$ \frac{\partial^2\left(x^2+y^2\right)}{\partial x^2}+\frac{\partial^2\left(x^2+y^2\right)}{\partial y^2}=4 $$

Therefore, the problem we need to solve numerically becomes:

$$ \frac{\partial^2 T}{\partial x^2}+\frac{\partial^2 T}{\partial y^2}=4 $$

The manufactured solution T_exact can then be used to check the accuracy of our numerical scheme.

Finally, let's implement a simple solver and check with our manufactured solution, if it works.

import numpy as np

# Set grid dimensions
nx, ny = 20, 20
x = y = np.linspace(-1, 1, nx)
dx = dy = x[1] - x[0]

# Set source term
source = 4.0

# Set up the exact solution for comparison
X, Y = np.meshgrid(x, y, indexing='ij')
T_exact = X**2 + Y**2

# Initialize solution grids
T = np.zeros((nx, ny))

# Set boundary conditions for manufactured solution
bdry = np.ones_like(T, dtype=bool)
bdry[1:-1, 1:-1] = False
T[bdry] = X[bdry]**2 + Y[bdry]**2

# Iterate until convergence, comparing with exact solution
max_diff = 1e-9
max_iter = 10000
it = 0
diff = 1.0
while True:
    T_new = T.copy()
    T_new[1:-1, 1:-1] = 0.25 * (T[:-2, 1:-1] +
                                T[2:, 1:-1] +
                                T[1:-1, :-2] +
                                T[1:-1, 2:] - dx**2 * source)
    diff = np.max(np.abs(T_new - T))
    if diff < max_diff:
        break
    if it == max_iter:
        raise Exception('No convergence')
    T = T_new
    it += 1

err = np.max(np.abs(T_exact - T))
print(it, err)

1160 7.258984531348378e-08

The problem is defined on a 20x20 grid with x and y coordinates ranging from -1 to 1. The grid spacing dx and dy are computed from these equidistant coordinate arrays. The source term is set to 4.0, according to the manufactured solution for the Poisson equation.

An exact solution T_exact is created using the manufactured solution T_exact, for later comparison with the numerical solution.

The boundary conditions are set based on the manufactured solution. To achieve this, a boolean array bdry is created to conveniently select the boundary. The boundary values of T are then filled using the formula for the exact solution.

The solver enters an iterative loop where the Jacobi method is applied until the maximum difference between the new and old solutions at each grid point is below a predefined threshold max_diff. Alternatively, if the solver hasn't converged after max_iter iterations, it raises an exception.

In each iteration, the new temperature T_new is computed using the Jacobi formula. This involves averaging the values from the four neighboring points and subtracting the scaled source term.

The difference between the old and new solution is then computed and compared to the threshold. If it’s below the threshold, the loop breaks, and the solver has converged. If not, the old solution T is replaced with the new solution, and the solver proceeds to the next iteration.

After convergence, the maximum absolute difference between the exact and numerical solutions (the error) is computed and printed, along with the number of iterations taken to converge.

Conclusion

In conclusion, the Method of Manufactured Solutions is an essential tool for verifying the accuracy and stability of numerical methods applied to complex problems. By formulating a theoretical solution that fits within the constraints of a problem and using it to manufacture a source term, we can establish a definitive benchmark for evaluating our computational methods.

Our practical application of the MMS to the 2D steady-state heat conduction equation underpins its effectiveness. We generated an artificial solution, adapted it to our problem, and successfully used it to test the accuracy of our numerical solver.