01 — Fall 2024

Numerical Analysis

These are notes & functions from Numerical Analysis, 10th Edition by Richard L. Burden, Douglas J. Faires, and Annette M. Burden. Some functions are slightly altered for functionality (and use numpy imported as np).

Lagrange Interpolation

P_n = L₀ f(x₀) + L₁ f(x₁) + ... + L_n f(x_n)

The Langrange polynomial interpolates the given nodes, because the structure makes it such that when calculating the polynomial at any node x_i, all L become 0 except L_i , which becomes 1.

f(x) = P_n + ¹/_(n+1)! f ⁽ⁿ⁺¹⁾(ζ) Π(x-x_i)

Newton's Divided Differences

P_n = f(x₀) + Δ^kf(x₁) + ... + Δⁿf(x_n)

This is another method of calculating interpolating polynomials, using a system of divided differences.

def divdiff(x, f):
    '''
    returns the coefficients Δ^0, Δ^1, ... Δ^n
    '''
    n = len(f)
    delta = np.zeros([n, n])
    delta[:,0] = f

    for j in range(1,n):
    for i in range(n-j):
        delta[i][j] = (delta[i+1][j-1] - delta[i][j-1]) / (x[i+j]-x[i])

    return delta

Hermite Interpolation

H_2n+1 = H_i,i + H_1,1(x-x₀) + ... +

The Hermite polynomial osculates the given points, meaning it matches the function value and the first derivative value at each node. Therefore it has smaller error than the previous interpolating polynomials:

f(x) = H_2n+1 + ¹/_(2n+1)! f ⁽²ⁿ⁺¹⁾(ζ) Π(x-x_i)²

These essentially act as the coefficients in the Divided Difference method, only accounting for 2n+1th degree of accuracy considering the use of known derivatives.

def hermite(x,f,d):
    '''
    prints the coefficients H_(i,i)
    '''
    n = len(f)
    q = np.zeros([2*n + 1, 2*n + 1])
    z = np.zeros(2*n + 1)
    
    for i in range(0, n):
        z[2*i] = x[i]
        z[2*i + 1] = x[i]
        q[2*i][0] = f[i]
        q[2*i + 1][0] = f[i]
        q[2*i + 1][1] = d[i]
        if i!=0:
        q[2*i][1] = (q[2*i][0] - q[2*i - 1][0]) / (z[2*i] - z[2*i - 1])
    
    for i in range(2, 2*n + 1):
        for j in range(2,i+1):
        q[i][j] = (q[i][j - 1] - q[i - 1][j -1]) / (z[i] - z[i - j])
    
    for i in range(0, 2*n):
        print(q[i][i])

Cubic Splines

S_i(x) = a_i + b_ix + c_ix² + d_ix³

Cubic splines are a piecewise solution: one third-degree polynomial for each pair of neighbor nodes, where the function value, first derivative value, and second derivative value all match at each node. This is used more in computer graphics rather than to approximate actual functions.

def natspline(x,a):
    n = len(x)
    i = np.zeros(2*n)
    u = np.zeros(2*n)
    z = np.zeros(2*n)
    h = np.zeros(2*n)
    alph = np.zeros(2*n)
    
    b = np.zeros(n+1)
    c = np.zeros(n+1)
    d = np.zeros(n+1)
    
    for j in range(0,n-1):
        h[j] = x[j+1] - x[j]
    
    for j in range(1,n-1):
        alph[j] = (3/h[j]) * (a[j+1] - a[j]) - (3/h[j-1]) * (a[j] - a[j-1])
    
    i[0] = 1
    u[0] = 0
    z[0] = 0
    
    for j in range(1,n-1):
        i[j] = 2*(x[j+1] - x[j-1]) - h[j-1] * u [j-1]
        u[j] = h[j] / i[j]
        z[j] = (alph[j] - h[j-1]*z[j-1])/i[j]
    
    i[n] = 1
    z[n] = 0
    c[n] = 0
    
    for j in reversed(range(0,n-1)):
        c[j] = z[j] - u[j]*c[j+1]
        b[j] = (a[j+1] - a[j])/h[j] - h[j]*(c[j+1] + 2*c[j])/3
        # d[j] = (a[j+1] - a[j])/h[j] - h[j]*(c[j+1] + 2*c[j])/3
        d[j] = (c[j+1]-c[j])/(3*h[j])
    
    for j in range(0,n-1):
        print(a[j],b[j],c[j],d[j])

Bezier Curves

x(t) = a₀ + a₁t + a₂t² + a₃t³

y(t) = b₀ + b₁t + b₂t² + b₃t³

Bezier curves parameterize lines in order to bypass the one-to-one limits of functions (allowing for shapes like loops and circles).

def bezier(n,x,y,left,right):
    a = np.zeros([n,4])
    b = np.zeros([n,4])
    
    i=0
    while i < n:
        a[i][0] = x[i]
        a[i][1] = 3*(left[i][0] - x[i])
        a[i][2] = 3*(x[i] + right[i][0] - 2*left[i][0])
        a[i][3] = x[i+1] - x[i] + 3*left[i][0] - 3*right[i][0]
    
        b[i][0] = y[i]
        b[i][1] = 3*(left[i][1] - y[i])
        b[i][2] = 3*(y[i] + right[i][1] - 2*left[i][1])
        b[i][3] = y[i+1] - y[i] + 3*left[i][1] - 3*right[i][1]
    
        print(i, a[i], b[i])
        i += 1