Routine Name: Thomas Algorithm, Elliptic ODE, and Error Vector Norm
Author: Kyle Hovey
Language: C++
Description/Purpose:
This algorithm uses a finite difference approximation of a second order derivative to solve the equation:
This equation has solution:
For the interval \( (0, 1) \), initial conditions \( u(0) = 2.5, u(1) = 5 \), and driving function \( f(x) = \sin(\pi x) \), the exact solution becomes:
Input:
The interval \( (a, b) \), \( u(a), u(b) \), the driving function \( f \), and the size of the mesh.
Output:
A column vector representing the approximation of the solution \( u \) at the mesh points.
Usage/Example:
int main() {
// Mesh size
const uint meshSize = 10;
// Limits
const double a = 0;
const double b = 1;
// Function u at limits
const double ua = 2.5;
const double ub = 5.0;
// Driving function
std::function<double(double)> f = [](double x) {
return std::sin(M_PI * x);
};
// Generate exact solution vector
std::function<double(double)> exact = [](double x) {
return 2.5 * x + 2.5 - std::sin(M_PI * x) / (M_PI * M_PI);
};
const double h = (b - a) / (double) meshSize;
Mtx uExact(meshSize - 1, 1, [&](const uint& i, const uint& j) {
(void) j;
return exact(a + (i + 1) * h);
});
const auto soln = solveElliptic<double>(a, b, ua, ub, f, meshSize);
std::cout << "Solved solution\n";
std::cout << soln << std::endl;
std::cout << "Exact solution\n";
std::cout << uExact << std::endl;
std::cout << "Error vector\n";
const auto E = soln - uExact;
std::cout << E << std::endl;
std::cout << "1-norm of error vector\n";
std::cout << Matrix::Matrix<double>::vNorm(E, 1) << std::endl;
std::cout << "2-norm of error vector\n";
std::cout << Matrix::Matrix<double>::vNorm(E, 2) << std::endl;
std::cout << "infinity-norm of error vector\n";
std::cout << Matrix::Matrix<double>::vNorm(E, std::numeric_limits<uint>::max()) << std::endl;
return EXIT_SUCCESS;
}
Output:
Solved solution
2.71843
2.93995
3.16735
3.40284
3.64784
3.90284
4.16735
4.43995
4.71843
Exact solution
2.71869
2.94044
3.16803
3.40364
3.64868
3.90364
4.16803
4.44044
4.71869
Error vector
-0.00025879
-0.000492248
-0.000677521
-0.000796474
-0.000837462
-0.000796474
-0.000677521
-0.000492248
-0.00025879
1-norm of error vector
0.00528753
2-norm of error vector
0.00187262
infinity-norm of error vector
-0.00025879
Implementation/Code:
From the first problem on HW 2, I created a method that solves for finite difference coefficients of arbitrary order. To solve this code, I created a function that puts those coefficients into a matrix such that the coefficients are staggered across the diagonal of the matrix. For instance, we need a second order derivative here which yields a tri-diagonal matrix with \(1, -2, 1\) as the coefficients across the diagonal.
To handle boundary conditions, I subtract them off at the boundaries (as seen in the textbook). Then, to see the error of the method, I compute the exact solution over the mesh chosen and take the 2-norm of the difference between the exact solution and the approximated solution. I do this using the Thomas Algorithm, which implements a linear scan across the diagonal. Since the matrix is of dimension \(n \times n\), then the diagonal is of size \(n\) and the Thomas Algorithm solves this system in:
Following is the code required to complete these goals:
/**
* Solve u'' = f given boundary conditions
* @param a Left boundary
* @param b Right boundary
* @param ua u(a)
* @param ub u(b)
* @param f Driving function
* @param n Size of mesh
* @return Column vector of solution
*/
template <typename T>
Matrix::Matrix<T> solveElliptic(
const T& a,
const T& b,
const T& ua,
const T& ub,
const std::function<T(T)>& f,
const uint& n
) {
(void) ua;
(void) ub;
// Generate F vector
const T h = (b - a) / (T) n;
const Matrix::Matrix<T> F(n - 1, 1, [&](const uint& i, const uint& j) -> T {
(void) j;
return std::pow(h, 2) * f(a + (i + 1) * h);
});
// Use initial conditions
Matrix::Matrix<T> init(n - 1, 1);
init.setVal(0, 0, ua);
init.setVal(n - 2, 0, ub);
// Generate differential operator for second derivative
const auto D = Matrix::Matrix<T>::genFDMatrix(n - 1, 2);
// Solve for solution vector
const auto soln = Matrix::Matrix<T>::solve(
D,
F - init,
Matrix::Solve::Thompson
);
return soln;
}
Where the thomas algorithm is implemented as follows (this is just the Thomas part of the general solve
function):
template <typename T>
Matrix<T> Matrix<T>::solve(
const Matrix<T>& A,
const Matrix<T>& b
) {
// Only work for square matrices
if (!A.isSquare()) {
throw std::domain_error("Input matrix must be square.");
}
// Get size of matrix
const auto m = std::get<0>(A.getSize());
if (!A.isNDiagonal(3)) {
throw std::domain_error("Thomas method needs tri-diagonal matrix.");
}
// Throw diagonals into vectors (store in a vector of diags)
std::vector<std::vector<T>> diags(3, std::vector<T>());
diags[0].push_back(0);
for (uint row = 0; row < m; ++row) {
for (uint nDiag = 0; nDiag < 3; ++nDiag) {
const auto pos = row - 1 + nDiag;
if (pos < m) {
diags[nDiag].push_back(A.getVal(row, (uint) pos));
}
}
}
diags[2].push_back(0);
// Syntactic nicety
auto _a = diags[0];
auto _b = diags[1];
auto _c = diags[2];
// Copy solution vector
auto _d = std::vector<T>();
for (uint i = 0; i < m; ++i) {
_d.push_back(b.getVal(i, 0));
}
// Do elimination via algebra
_c[0] = _c[0] / _b[0];
_d[0] = _d[0] / _b[0];
for (uint i = 1; i < m; ++i) {
_c[i] = _c[i] / (_b[i] - _a[i] * _c[i - 1]);
_d[i] = (_d[i] - _a[i] * _d[i - 1]) / (_b[i] - _a[i] * _c[i - 1]);
}
// Create the solution
Matrix<T> x(m, 1);
x.setVal(m - 1, 0, _d[m - 1]);
for (int i = m - 2; i >= 0; --i) {
const uint _i = i;
x.setVal(_i, 0, _d[_i] - _c[_i] * x.getVal(_i + 1, 0));
}
return x;
}
The vector norm is calculated using a general function that will calculate any \( p \)-norm for a vector of arbitrary size. If \( p \) is the maximum value for its type, then the infinity norm is given instead.
template <typename T>
T Matrix<T>::vNorm(const Matrix<T>& v, const uint& n) {
const auto [ M, N ] = v.getSize();
if (M != 1 && N != 1 && M != N) {
throw std::domain_error("Need a row or column vector for vector norm.");
}
const auto isRow = (M == 1);
const auto size = isRow ? N : M;
T sum = 0;
T max = std::numeric_limits<T>::min();
for (uint i = 0; i < size; ++i) {
const auto val = isRow ? v.getVal(0, i) : v.getVal(i, 0);
max = (val > max) ? val : max;
sum += pow(std::abs(val), n);
}
if (n == std::numeric_limits<uint>::max()) {
// Infinity norm
return max;
} else {
return pow(sum, 1.0 / n);
}
}