Eigenvalue problems - Executable Engineering: Interactive Numerical Methods in Python

Numerical methods (e.g.: IVPs) will quickly give you a specific solution to a specific problem, and for complex systems sometimes that is all you can hope for.

Sometimes they can give you information about all the solutions! Nondimensionalization is one way, eigenvalue analysis is another.

The Eigenvalue problem is:

A \vec{x} = \lambda \vec{x}

(1)

with $A$ the operator matrix, $\lambda$ the eigenvalue with corresponding vector $\vec{x}$ . This means that the operator $A$ mearly stretches vector $\vec{x}$ by $\lambda$ in such a way that you may consider the set of eigenvectors as a natural basis for the operator $A$ . There is an excellent video from 3 Blue, 1 Brown on the subject.

Example of general eigenvectors¶

# prompt: make a function that plots 8 points on the unit circle marked in different colours and joined with a solid thin line, transforms them all by a given matrix A and then plots the result overlapping in the same colours with a broken line. . Make an ipython widget that takes the values of the matrix as text boxes and then calls the function on it

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

def plot_transformation(a, b, c, d):
    # Generate n points on the unit circle
    num_points = 17
    angles = np.linspace(0, 2*np.pi, num_points, endpoint=True)
    points = np.exp(1j*angles)
    x = np.real(points)
    y = np.imag(points)

    # Transformation matrix
    A = np.array([[a, b], [c, d]])

    # Apply the transformation
    transformed_points = np.dot(A, np.array([x, y]))
    transformed_x = transformed_points[0, :]
    transformed_y = transformed_points[1, :]

    # Plotting
    plt.figure(figsize=(6, 6))
    colors = plt.cm.jet(np.linspace(0, 1, num_points))

    plt.plot(x, y, marker='o', linestyle='-', color='k', linewidth=0.5, markerfacecolor='black', markeredgecolor='black')
    plt.plot(transformed_x, transformed_y, marker='o', linestyle='--', color='k', linewidth=0.5, markerfacecolor='black', markeredgecolor='black')
    for i in range(num_points):
        plt.plot(x[i], y[i], marker='o', markersize=8, color=colors[i])
        plt.plot(transformed_x[i], transformed_y[i], marker='o', markersize=8, color=colors[i])

    w, v = np.linalg.eig(A)
    print(v)
    # Plot eigenvectors
    plt.quiver(0, 0, v[0,0], v[1,0], angles='xy', scale_units='xy', scale=1/w[0], color='r', label='Eigenvector 1')
    plt.quiver(0, 0, v[0,1], v[1,1], angles='xy', scale_units='xy', scale=1/w[1], color='g', label='Eigenvector 2')

    plt.xlim(-2, 2)
    plt.ylim(-2, 2)
    plt.gca().set_aspect('equal')
    plt.title('Transformation of Points on the Unit Circle')
    plt.xlabel('X')
    plt.ylabel('Y')
    plt.grid(True)
    plt.show()

# Create interactive widgets for matrix elements
a_widget = widgets.FloatText(value=1, description='a:')
b_widget = widgets.FloatText(value=0, description='b:')
c_widget = widgets.FloatText(value=0, description='c:')
d_widget = widgets.FloatText(value=1, description='d:')

interact(plot_transformation, a=a_widget, b=b_widget, c=c_widget, d=d_widget);

Example: Symmetric matricies¶

# prompt: make a function that plots 17 points on the unit circle marked in different colours and joined with a solid thin line, transforms them all by a given symmetric matrix A and then plots the result overlapping in the same colours with a broken line. put that plot beside a surface plot of the function with hessian A, and draw lines along the surface in the eigenvector directions. Make an ipython widget that takes the values of the matrix as text boxes and then calls the function on it

import numpy as np
import matplotlib.pyplot as plt
from ipywidgets import interact, interactive, fixed, interact_manual
import ipywidgets as widgets
from mpl_toolkits.mplot3d import Axes3D

def plot_transformation(a, b, d):
    # Generate n points on the unit circle
    c = b
    num_points = 17
    angles = np.linspace(0, 2*np.pi, num_points, endpoint=True)
    points = np.exp(1j*angles)
    x = np.real(points)
    y = np.imag(points)

    # Transformation matrix
    A = np.array([[a, b], [c, d]])

    # Apply the transformation
    transformed_points = np.dot(A, np.array([x, y]))
    transformed_x = transformed_points[0, :]
    transformed_y = transformed_points[1, :]

    # Plotting
    fig = plt.figure(figsize=(12, 6))

    # Left plot: point transformations
    ax1 = fig.add_subplot(121)
    colors = plt.cm.jet(np.linspace(0, 1, num_points))
    ax1.plot(x, y, marker='o', linestyle='-', color='k', linewidth=0.5, markerfacecolor='black', markeredgecolor='black')
    ax1.plot(transformed_x, transformed_y, marker='o', linestyle='--', color='k', linewidth=0.5, markerfacecolor='black', markeredgecolor='black')
    for i in range(num_points):
        ax1.plot(x[i], y[i], marker='o', markersize=8, color=colors[i])
        ax1.plot(transformed_x[i], transformed_y[i], marker='o', markersize=8, color=colors[i])

    w, v = np.linalg.eig(A)
    ax1.quiver(0, 0, v[0,0], v[1,0], angles='xy', scale_units='xy', scale=1/w[0], color='r', label='Eigenvector 1')
    ax1.quiver(0, 0, v[0,1], v[1,1], angles='xy', scale_units='xy', scale=1/w[1], color='g', label='Eigenvector 2')
    ax1.set_xlim(-2, 2)
    ax1.set_ylim(-2, 2)
    ax1.set_aspect('equal')
    ax1.set_title('Transformation of Points')
    ax1.grid(True)

    # Right plot: surface plot
    ax2 = fig.add_subplot(122, projection='3d')
    X = np.arange(-2, 2, 0.1)
    Y = np.arange(-2, 2, 0.1)
    X, Y = np.meshgrid(X, Y)
    Z = 0.5 * (a*X**2 + 2*b*X*Y + d*Y**2) # Function with Hessian A


    # Line along eigenvector 1
    t = np.linspace(-2,2,50)
    x_line1 = v[0,0] * t
    y_line1 = v[1,0] * t
    z_line1 = 0.5 * (a*x_line1**2 + 2*b*x_line1*y_line1 + d*y_line1**2)
    ax2.plot(x_line1, y_line1, z_line1*1.1, color='r',  linewidth=5, alpha=1,label='Eigenvector 1')


    # Line along eigenvector 2
    x_line2 = v[0,1] * t
    y_line2 = v[1,1] * t
    z_line2 = 0.5 * (a*x_line2**2 + 2*b*x_line2*y_line2 + d*y_line2**2)
    ax2.plot(x_line2, y_line2, z_line2*1.1, color='g', linewidth=5,alpha=1, label='Eigenvector 2')

    ax2.plot_surface(X, Y, Z)
    plt.show()


# Create interactive widgets for matrix elements
a_widget = widgets.FloatText(value=1, description='a:')
b_widget = widgets.FloatText(value=0, step=.1, description='b,c:')
d_widget = widgets.FloatText(value=1, description='d:')

interact(plot_transformation, a=a_widget, b=b_widget, d=d_widget);

In general:

Eigenvalues are easier to calculate and don’t require solving the problem
Eigenvectors require solving the system which can be signkificantly more computationally expensive.