Basic Fix Point Iteration

This exercise is about the exponential function expressed in the Taylor Series \[ \exp(x) \approx f^{(N)}(x) = \sum_{n=0}^{N} \frac{1}{n!} x^{n} \] and to find the highest eigenvalue, the lowest eigenvalue and a specific eigenvalue closest to a constant \(\lambda \).

Theory

Taylor Series

The Taylor series of a real or complex-valued function f(x), that is infinitely differentiable at a real or complex number a, is the power series \[ f(a) + \frac{f'(a)}{1!} (x-a) + \frac{f''(a)}{2!} (x-a)^2 + \frac{f'''(a)}{3!} (x-a)^3 + \ldots = \sum_{n=0}^{\infty} \frac{f^{(n)}(a)}{n!} (x-a)^n \]

With \(a = 0\), the Maclaurin series takes the form \[ f(0) + \frac{f'(0)}{1!} (x) + \frac{f''(0)}{2!} (x)^2 + \frac{f'''(0)}{3!} (x)^3 + \ldots = \sum_{n=0}^{\infty} \frac{f^{(n)}(0)}{n!} (x)^n \]

It is relevant to look at the Macluarian series for the exponent function (so take \(a = 0 \)), because of the first term gives \(f(0) = 1 \) and the derivative of the exponent function is the exponent function again \(\frac{d}{dx} \left(e^{x} \right) = e^{x} \).

Deriving the Taylor (Maclaurin) series of the exponent function

First we choose to expand around the point \(a = 0 \). We then have to find \(f(0)\), \(f'(0)\), \(f''(0) \), etc. for the exponential function. \[f(0) = e^{0} = 1 \] \[f'(0) = e^{0} = 1 \] \[f''(0) = e^{0} = 1 \] etc. etc. ... so we now can express the Taylor (Macluarian) series as \[f(x) = 1 + \frac{1}{1!} (x) + \frac{1}{2!} (x)^2 + \frac{1}{3!} (x)^3 + \ldots \] \[f(x) = 1 + x + \frac{1}{2} x^2 + \frac{1}{6} x^3 + \ldots \]

This can be simplified as a series \[f(x) = \exp(x) = \sum_{n=0}^{\infty} \frac{1}{n!} x^{n} \]

Relative error and accuracy

Relative error is the ratio of the absolute error and the actual value (so a measure of accuracy) in this case. \[\mathbf{RE}_{\text{accuracy}} = \frac{\text{Absolute Error}}{\text{Actual/True Error}} \] We will expect that for each additional term, we will come closer to the actual value of the exponential function. Since the actual value is the sum with infinite many terms.

Implementation

This implementation involves the approximation of the exponential function from a series to a finite sum with \(N\) terms. We have to choose the following

The positive and negative values of \(x\)
Number of terms \(N \)
A trick to get a good accuracy for all points \(a\)

We have to define a function, e.g. \(\text{approx_exp(x, N)}\) to approximate through Taylor series. We can also define a function to calculate the relative error or just make a loop.

What I did:

Values of \(x\): We pick a few positive values and negative values. The lecturer took 1,3 and 10 and the same for the negatives, -1, -3 and -10.
Number of terms \(N\): It should drop down fast, but the lecturer took 100 terms. We might need more if we increase the values of \(x\).
Trick for good accuracy for all points \(a \): Since we know \(a = 0 \) is a good point for Taylor series. As we move farther from zero (either positive or negative), the terms in the series grow, and we need significantly more terms to maintain accuracy. For large values the truncated series cannot capture the rapid growth or decay and leading to higher errors.

Results

The results from the command line are shown in Figure 2 and corresponds with the given values of \(\lambda \)

Figure 1: Plot of basic fix point iterations.

Figure 2: Results of basic fix point iterations.

For the matrix \(\mathbf{B}\) \[ \mathbf{B} = \frac{1}{5} \left(\begin{matrix} 3 & 4 & 0 \\ -4 & 3 & 0 \\ 0 & 0 & 5 \end{matrix} \right) \]

The convergence plot (with 10 iterations) showing the three different iteration methods are shown in Figure 3. The results from the command line are shown in Figure 4

Code

Show/Hide Python Code

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import factorial  # Use scipy's factorial for array support

# Values of x
x = np.array([1, 3, 10], dtype=float)
# x = np.array([-1, -3, -10], dtype=float)
term_size = 60
err_list = []

# Function to approximate exp(x) using Taylor series expansion
def approx_exp(x, N):
    y = np.zeros_like(x, dtype=float)  # Initialize y as float type
    for n in range(N):
        y += (x ** n) / factorial(n, exact=False)  # Use scipy factorial for array operations
    return y

# Relative error
for N in range(1, term_size + 1):
    err = np.abs(approx_exp(x, N) - np.exp(x)) / np.exp(x)
    err_list.append(err)

# Convert list of errors to a numpy array for easier plotting
err_list = np.array(err_list)

# Plotting the results
plt.semilogy(range(1, term_size + 1), err_list)
plt.xlabel('Number of terms N')
plt.ylabel('Relative error')
plt.title('Convergence of Taylor Series Approximation of exp(x)')
plt.legend(['x=1', 'x=3', 'x=10'])
plt.grid(True)
plt.show()

Show/Hide MATLAB Code

clear all;
close all;

% Values of x
x = [1, 3, 10];
%x = [-1, -3, -10];
term_size = 60;
err_list = [];

% Relative error
for N = 1:term_size
    err = abs(approx_exp(x, N) - exp(x)) ./ exp(x);
    %err = abs(1./approx_exp(-x, N) - exp(x)) ./ exp(x);

    err_list = [err_list; err];
end

semilogy(1:term_size, err_list);  % Added 1:100 to the x-axis for clarity
xlabel('Number of terms N');
ylabel('Relative error');
title('Convergence of Taylor Series Approximation of exp(x)');
legend('x=1', 'x=3', 'x=10');

function y = approx_exp(x, N)
    % This function approximates the exponential function exp(x)
    % using the first N terms of the Taylor series expansion.

    y = zeros(size(x));  % Initialize the output to be the same size as x
    
    for n = 0:N-1
        y = y + x.^n / factorial(n);
    end
end

Show/Hide C++ Code

#include < iostream >
#include < Eigen/Dense >
#include < vector >
#include "matplotlibcpp.h"

namespace plt = matplotlibcpp;
using namespace Eigen;
using namespace std;

int main() {
    // Define the matrix A
    MatrixXd A(4, 4);
    A << 1, 2, 3, 4,
         1, 5, 6, 7,
         1, 1, 8, 9,
         1, 1, 1, 0;

    int iter = 100;
    int N = A.rows();

    // Fix point iteration
    VectorXd y = VectorXd::Random(N);
    vector lA;
    vector change;

    for (int i = 0; i < iter; ++i) {
        y.normalize();
        VectorXd x = A * y;
        double lambda = (y.transpose() * A * y) / (y.transpose() * y);
        lA.push_back(lambda);
        change.push_back((y - x / x.norm()).norm());
        y = x;
    }
    cout << "Fix Point Iteration Eigenvalue: " << (y.transpose() * A * y) / (y.transpose() * y) << endl;

    // Inverse iteration
    MatrixXd B = A.inverse();
    y = VectorXd::Random(N);
    vector lB;

    for (int i = 0; i < iter; ++i) {
        y.normalize();
        VectorXd x = B * y;
        double lambda = (y.transpose() * A * y) / (y.transpose() * y);
        lB.push_back(lambda);
        y = x;
    }
    cout << "Inverse Iteration Eigenvalue: " << (y.transpose() * A * y) / (y.transpose() * y) << endl;

    // Shift-and-invert
    double mu = -1.4;
    MatrixXd C = (A - mu * MatrixXd::Identity(N, N)).inverse();
    y = VectorXd::Random(N);
    vector lC;

    for (int i = 0; i < iter; ++i) {
        y.normalize();
        VectorXd x = C * y;
        double lambda = (y.transpose() * A * y) / (y.transpose() * y);
        lC.push_back(lambda);
        y = x;
    }
    cout << "Shift-and-Invert Iteration Eigenvalue: " << (y.transpose() * A * y) / (y.transpose() * y) << endl;

    // Plotting
    plt::figure();
    plt::plot(lA, "ko-", {{"label", "Fix Point Iteration"}});
    plt::plot(lB, "bo-", {{"label", "Inverse Iteration"}});
    plt::plot(lC, "ro-", {{"label", "Shift-and-Invert Iteration"}});
    plt::legend();
    plt::show();

    return 0;
}

Evaluation of the exponential function via the Taylor Series

Intro

Taylor Series

Deriving the Taylor (Maclaurin) series of the exponent function

Relative error and accuracy

Code