Eigenfunctions - Executable Engineering: Interactive Numerical Methods in Python

Let’s reinterpret eigenvectors in terms of our usage.

In previous chapters we saw how to:

represent a solution as a vector of unknows.
represent a differential operator (e.g. a PDE) as:
- a matrix (for linear operators)
- a vector function (nonlinear operators) which could be solved by root finding, which involves linear operators.

Another way to think about this is that we represented continuous functions in a finite basis (weights / coordinates collected into our vector of unkowns). If we used infinite basis we could represent the full function, but then our vector would be infinite too.

Turn this around, and consider that eigenfunctions (the function approximated by the eigenvectors) are special functions which are not changed by the operator / differential equation.

You already know all about eigenfunctions!

Consider swinging on a swing - your frequency is independant of the amplitude
Consider resonance: If you pump your legs at the same frequency, you will go higher!

In our applications with PDEs, time stepping was an operator that transformed $y^t$ to $y^{t+1}$ . The eigenfunctions of that operator only change in amplitude, not shape!

Example: Coupled oscillators¶

\begin{align} m \ddot{x}_1 &=-k x_1 +k [x_2-x_1]\\ m \ddot{x}_2 &= -k[x_2-x_1] - k x_2 \end{align}

(1)

or in vector form with $\vec{x} = \begin{bmatrix} x_1 \\ x_2 \end{bmatrix}$

\ddot{\vec{x}} = -\frac{k}{m} \begin{bmatrix} 2 & -1 \\ -1 & 2 \end{bmatrix} \vec{x}

(2)

which is just an IVP which we know how to solve!

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
from ipywidgets import interact, interactive, fixed, interact_manual
import ipywidgets as widgets

# Define the system of differential equations

def coupled_oscillators(t, y, k, m):
    x1, x2, v1, v2 = y
    dx1dt = v1
    dx2dt = v2
    dv1dt = -(k/m) * (2*x1 - x2)
    dv2dt = -(k/m) * (-x1 + 2*x2)
    return [dx1dt, dx2dt, dv1dt, dv2dt]

# Define parameters
k = 1  # Spring constant
m = 1  # Mass

# Time span for integration
t_span = (0, 100)

# Function to solve and plot
def solve_and_plot(x1_0, x2_0, v1_0=0, v2_0=0):
    # Initial conditions
    y0 = [x1_0, x2_0, v1_0, v2_0]

    # Solve the IVP
    sol = solve_ivp(coupled_oscillators, t_span, y0, args=(k, m), max_step = .1,  dense_output=True)

    # Time points for plotting
    t = np.linspace(t_span[0], t_span[1], 500)

    # Get the solution at the specified time points
    x1, x2, v1, v2 = sol.sol(t)


    # Calculate FFT
    N = len(t)
    dt = t[1] - t[0]
    xf = np.fft.fftfreq(N, d=dt)
    x1_fft = np.fft.fft(x1)
    x2_fft = np.fft.fft(x2)


    plt.figure(figsize=(10,6))

    # Plot the results
    plt.subplot(2,1,1)
    plt.plot(t, x1, label='x1(t)')
    plt.plot(t, x2, label='x2(t)')
    plt.xlabel('Time')
    plt.ylabel('Displacement')
    plt.title('Coupled Oscillators')
    plt.legend()
    plt.grid(True)
    plt.show()

    #Plot FFT
    plt.subplot(2,1,2)
    plt.plot(xf, np.abs(x1_fft), label='FFT of x1(t)')
    plt.plot(xf, np.abs(x2_fft), label='FFT of x2(t)')
    plt.xlabel('Frequency')
    plt.ylabel('Magnitude')
    plt.title('FFT of Coupled Oscillators')
    plt.legend()
    plt.xlim(-5,5)
    plt.grid(True)
    plt.show()

# Create interactive plot with sliders for initial conditions
interact(solve_and_plot, x1_0=widgets.FloatSlider(min=-2, max=2, step=0.1, value=1),
         x2_0=widgets.FloatSlider(min=-2, max=2, step=0.1, value=0));

Loading...

Note the ineresting behaviour of $x_1 = x_2$ and $x_1 = -x_2$ !

So what is happenning?

Eigenvector form¶

Say the solution is of the form $x = e^{-i\omega t}$ :

\begin{align} \vec{x} &= e^{-i \omega t} \\ \ddot{\vec{x}} &= -\omega^2 e^{-i \omega t} = -\omega^2 \vec{x} \end{align}

(3)

and our equation becomes:

-\omega^2 \rm{I} \ \vec{x} = -\frac{k}{m} \begin{bmatrix} 2 & -1 \\ -1 & 2 \end{bmatrix} \vec{x} \\ \bigg[\begin{bmatrix} 2 & -1 \\ -1 & 2 \end{bmatrix} -\frac{m \omega^2}{k} \rm{I}\bigg] \vec{x} = 0

(4)

The solution to this equation must be:

Trivial, $\vec{x} = \vec{0}$
Non-zero, but then the matrix must be singular!

For a matrix $A - \lambda \rm{I}$ to be singular, $det(A-\lambda \rm{I}) = 0$ . The scalar values $\lambda$ are the eigenvalues of $A$ , and the corresponding $\vec{x}$ is the eigenvector.

In our example, the eigenvalues of $\begin{bmatrix} 2 & -1 \\ -1 & 2 \end{bmatrix}$ are $\frac{m \omega^2}{k}$ . Note that we have grouped the unknowns so we really find a general solution!

The determinant is convenient because it is simply a polynomial, and we know how to find roots of a polynomial!

# prompt: Write the above equation in sympy and find en expression for hte determinant.

import sympy as sp
import numpy as np

# Define symbolic variables
lam = sp.symbols('lambda')

# Define the matrix
matrix = sp.Matrix([[2, -1], [-1, 2]])

# Define the identity matrix
identity_matrix = sp.eye(2)

# Form the matrix A - λI
matrix_minus_lambda_I = matrix - lam * identity_matrix

# Calculate the determinant
determinant = matrix_minus_lambda_I.det()

# Display the determinant
print("Determinant:", determinant)

# Convert the symbolic determinant to polynomial coefficients
poly = sp.Poly(determinant, lam)
coeffs = poly.all_coeffs()
print("Coefficients of the polynomial:", coeffs)

# Solve for lambda using numpy.roots
roots = np.roots(coeffs)
print("Roots (eigenvalues) of the characteristic polynomial:", roots)

Determinant: lambda**2 - 4*lambda + 3
Coefficients of the polynomial: [1, -4, 3]
Roots (eigenvalues) of the characteristic polynomial: [3. 1.]

Let’s check the Eigenvalues, then use them to find the corresponding vectors!

from scipy.optimize import root
A = np.array([[2, -1], [-1, 2]])
for l in [1,3]:
  print(np.linalg.det((A-l*np.eye(2))))
  #print(np.linalg.solve(A-1*np.eye(2), np.zeros(2)))
  print(root(lambda x: x@(A-l*np.eye(2)), np.array([[1,1]])).x)

0.0
[1. 1.]
0.0
[ 1. -1.]

The eigenvalues plugged back into $x = e^{-i\omega t}$ implies the two harmonics are:

\begin{align} x &= e^{-i\sqrt{\frac{k}{m}} t} \\ x &= e^{-i\sqrt{3 \frac{k}{m}} t} \end{align}

(5)

For the first harmonic: $\frac{m \omega^2}{k} = 1$ :

\begin{align} \bigg[\begin{bmatrix} 2 & -1 \\ -1 & 2 \end{bmatrix} - \rm{I}\bigg] \vec{x} &= 0 \\ \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix} \vec{x} &= 0 \\ x_1 &= x_2 \end{align}

(6)

Similarly the second harmonic, $\frac{m \omega^2}{k} = 3$ :

\begin{align} \begin{bmatrix} -1 & -1 \\ -1 & -1 \end{bmatrix} \vec{x} &= 0 \\ x_1 &= -x_2 \end{align}

(7)

which is what we observed above.

This is all fine and good, but there must be a better way! For small matricies ( $n \le 4$ ) there are analytical solutions for the roots which we can use.

You may recal that the numerical roots function finds the roots as the eigenvalues of a companion matrix so it obviously isn’t the correct approach but does reveal the euqivilance between polynomial roots and eigenvalues!