Module NumericalMethods
Contents
9.3. Module NumericalMethods#
The following code cells are the content of file NumericalMethods.jl, used to define the module NumericalMethods
This file exists for two reasons:
It can be convenient to gather the cells defining various functions for the module in a notebook like this one, which allows documentation, and then convert to the the “.jl” file with the JupterLab command
File > Save and Export Notebook As ... > Executable ScriptThis gathers the contents of the code cells, ignorig any markdown cells.This description of the module’s definitions can be used as a section in a Jupyter Book.
Usage is:
include("NumericalMethods.jl")
then
using .NumericalMethods
or for a particular function, like
import .NumericalMethods: newtonmethod
# Module `NumericalMethods`
# Version of 2022-11-13
# The following code cells are the content of file `NumericalMethods.jl`, used to define the module `NumericalMethods`
# The notebook file version exists for two reasons:
#
# 1. It can be convenient to gather the cells defining various functions for the module in a notebook like this one,
# which allows documentation, and then convert to the the ".jl" file with the JupterLab command
# `File > Save and Export Notebook As ... > Executable Script`
#
# This gathers the contents of the code cells, ignorig any markdown cells.
# 2. This description of the module's definitions can be used as a section in a Jupyter Book.
module NumericalMethods
Root-finding#
Newton’s method#
function newtonmethod(f, Df, x0, errortolerance; maxiterations=20, demomode=false)
# Basic usage is:
# (rootApproximation, errorEstimate, iterations) = newton(f, Df, x0, errorTolerance)
# There is an optional input parameter "demomode" which controls whether to
# - println intermediate results (for "study" purposes), or to
# - work silently (for "production" use).
# The default is silence.
if demomode
println("Solving by Newton's Method.")
println("maxiterations = $maxiterations")
println("errortolerance = $errortolerance")
end
x = x0
global errorestimate # make it global to this function; without this it would be local to the "for" loop.
for iteration in 1:maxiterations
fx = f(x)
Dfx = Df(x)
# Note: a careful, robust code would check for the possibility of division by zero here,
# but for now I just want a simple presentation of the basic mathematical idea.
dx = fx/Dfx
x -= dx # Aside: this is shorthand for "x = x - dx"
errorestimate = abs(dx);
if demomode
println("At iteration $iteration, x = $x with estimated error $errorestimate and backward error $(abs(f(x)))")
end
if errorestimate <= errortolerance
if demomode
println("Done!")
end
return (x, errorestimate, iteration)
end
end
# Note: if we get to here (no "return" above), it completed maxIterations iterations without satisfying the accuracy target,
# but we still return the information that we have.
return (x, errorestimate, maxiterations)
end;
The secant method#
function secantmethod(f, a, b, errortolerance; maxiterations=20, demomode=false)
# Solve f(x)=0 in the interval [a, b] by the Secant Method.
if demomode
print("Solving by the Secant Method.")
end;
# Some more descriptive names
x_older = a
x_more_recent = b
f_x_older = f(x_older)
f_x_more_recent = f(x_more_recent)
for iteration in 1:maxiterations
global x_new, errorestimate
if demomode
println("\nIteration $(iteration):")
end;
x_new = (x_older * f_x_more_recent - f_x_older * x_more_recent)/(f_x_more_recent - f_x_older)
f_x_new = f(x_new)
(x_older, x_more_recent) = (x_more_recent, x_new)
(f_x_older, f_x_more_recent) = (f_x_more_recent, f_x_new)
errorestimate = abs(x_older - x_more_recent)
if demomode
println("The latest pair of approximations are $x_older and $x_more_recent,")
println("where the function's values are $f_x_older and $f_x_more_recent respectively.")
println("The new approximation is $x_new with estimated error $errorestimate and backward error $(abs(f_x_new))")
end;
if errorestimate < errortolerance
break
end;
end;
# Whether we got here due to accuracy of running out of iterations,
# return the information we have, including an error estimate:
return (x_new, errorestimate)
end;
Linear Algebra and Simultaneous Equations#
Row Reduction#
(with no pivoting)
function rowreduce(A, b)
# To avoid modifying the matrix and vector specified as input,
# they are copied to new arrays, with the function copy().
# Warning: it does not work to say "U = A" and "c = b";
# this makes these names synonyms, referring to the same stored data.
U = copy(A) # not "U=A", which makes U and A synonyms
c = copy(b)
n = length(b)
L = zeros(n, n)
for k in 1:n-1
for i in k+1:n
# compute all the L values for column k:
L[i,k] = U[i,k] / U[k,k] # Beware the case where U[k,k] is 0
for j in k+1:n
U[i,j] -= L[i,k] * U[k,j]
end
# Put in the zeros below the main diagonal in column k of U;
# this is not important for calculations, since those elements of U are not used in backward substitution,
# but it helps for displaying results and for checking the results via residuals.
U[i,k] = 0.
c[i] -= L[i,k] * c[k]
end
end
for i in 2:n
for j in 1:i-1
U[i,j] = 0.
end
end
return (U, c)
end;
Backward substitution#
function backwardsubstitution(U, c; demomode=false)
n = length(c)
x = zeros(n)
x[end] = c[end]/U[end,end]
if demomode
println("x_$n = $(x[n])")
end
for i in n-1:-1:1
if demomode
println("i=$i")
end
x[i] = ( c[i] - sum(U[i,i+1:end] .* x[i+1:end]) ) / U[i,i]
if demomode
print("x_$i = $(x[i])")
end
end
return x
end;
Solve a linear system (no pivoting)#
solvelinearsystem(A, b) = backwardsubstitution(rowreduce(A, b)...);
LU factorization#
function lu_factorize(A; demomode=false)
# Compute the Doolittle LU factorization of A.
# Sums like $\sum_{s=1}^{k-1} l_{k,s} u_{s,j}$ are done as matrix products;
# in the above case, row matrix L[k, 1:k-1] by column matrix U[1:k-1,j] gives the sum for a give j,
# and row matrix L[k, 1:k-1] by matrix U[1:k-1,k:n] gives the relevant row vector.
n = size(A)[1] # First component of the array's size; size(A) returns "(rows, columns)"
# Initialize U as a zero matrix;
# correct below the main diagonal, with the other entries to be computed and filled below.
U = zeros(n,n)
# Initialize L as a zero matrix;
# correct above the main diagonal, with the other entries to be computed and filled in below.
L = zeros(n,n)
# Column and row 1 are special:
U[1,:] = A[1,:]
L[1,1] = 1.
L[2:end,1] = A[2:end,1]/U[1,1]
if demomode
println("After step k=1")
println("U="); printmatrix(U)
println("L="); printmatrix(L)
end;
for k in 2:n-1
# Julia note: it is necessary to use indices "[k]" and so on to get a one-row matrix instead of a vector.
U[[k],k:end] = A[[k],k:end] - L[[k],1:k] * U[1:k,k:end]
L[k,k] = 1.
L[k+1:end,k] = ( A[k+1:end,k] - L[k+1:end,1:k] * U[1:k,k] )/U[k,k]
if demomode
println("After step k=$k")
println("U="); printmatrix(U)
println("L="); printmatrix(L)
end;
end;
# The last row is also special: nothing to do for L
L[end,end] = 1.
U[end,end] = A[end,end] - sum(L[[n],1:end-1] * U[1:end-1,end])
if demomode
println("After step k=$n")
println("U="); printmatrix(U)
end;
return [L, U]
end;
Forward substitution#
(without pivoting)
function forwardsubstitution(L, b)
# Solve L c = b for c.
n = length(b)
c = zeros(n)
c[1] = b[1]
for i in 2:n
c[i] = b[i] - sum(L[i:i,1:i] * c[1:i])
end;
return c
end;
PLU factorization#
function plu(A; demomode=false)
# Compute the Doolittle PA=LU factorization of A —
# but with the permutation recorded as permutation vector, not as the permutation matrix P.
# Sums like $\sum_{s=1}^{k-1} l_{k,s} u_{s,j}$ are done as matrix products;
# in the above case, row matrix L[k, 1:k-1] by column matrix U[1:k-1,j] gives the sum for a give j,
# and row matrix L[k, 1:k-1] by matrix U[1:k-1,k:n] gives the relevant row vector.
n = size(A)[1] # gives the number of rows in the 2D array.
π = zeros(Int, n)
# Julia can use Greek letters (and in fact, UNICODE):
# to insert character π, type \pi, hit tab, and select "π" from the menu.
# Or just call it "perm" or such.
π = collect(1:n)
# Julia language note: function "collect" converts the abstract entity "1:n" into an array of numbers.
# Initialize U as the zero matrix;
# correct below the main diagonal, with the other entries to be computed below.
U = zeros(n,n)
# Initialize L as zeros;
# correct above the main diagonal, with the other entries to be computed below,
# including the ones on the diagonal.
L = zeros(n,n)
for k in 1:n-1
if demomode; println("k=$k"); end
# Find the pivot element in column k:
pivotrow = k
abs_u_ik_max = abs(A[π[k],k])
for row in k+1:n
abs_u_ik = abs(A[π[row],k])
if abs_u_ik > abs_u_ik_max
pivotrow = row
abs_u_ik_max = abs_u_ik
end
end
if pivotrow > k # swap rows, virtually
if demomode; println("Swap row $k with row $pivotrow"); end
(π[k], π[pivotrow]) = (π[pivotrow], π[k])
else
if demomode; println("No row swap needed."); end
end
U[k,k:end] = A[[π[k]],k:end] - L[[π[k]],1:k] * U[1:k,k:end]
L[π[k],k] = 1.
for row in k+1:n
L[π[row],k] = ( A[π[row],k] - L[π[row],1:k] ⋅ U[1:k,k] ) / U[k,k]
# Julia note: To enter the centered dot notation for the dot product, type "\cdot" and then hit the tab key.
end
if demomode
println("permuted A is:")
for row in 1:n
println(A[π[row],:])
end
println("Intermediate L is"); printmatrix(L)
println("Intermediate U is"); printmatrix(U)
end
end
# The last row (index "end") is special: nothing to do for L except put in the 1 on the "permuted main diagonal"
L[π[end],end] = 1.
U[end,end] = A[π[end],end] - L[π[end],1:end-1] ⋅ U[1:end-1,end]
if demomode
println("After the final step, k=$(n-1)")
println("L is"); printmatrix(L)
println("U is"); printmatrix(U)
end
return (L, U, π)
end;
Forward substitution with pivoting#
function forwardsubstitution(L, b, π)
# Version 2: with permutation of rows
# Solve L c = b for c, with permutation of the rows of L and of b.
n = length(b)
c = zeros(n)
c[1] = b[π[1]]
for i in 2:n
c[i] = b[π[i]] - L[π[i], 1:i] ⋅ c[1:i]
end
return c
end;
Collocation and Data Fitting#
Polynomial collocation#
function polyfit(x, y)
# Version 1: exact collocation.
# Compute the coeffients c_i of the polynomial of lowest degree that collocates the points (x[i], y[i]).
# These are returned in an array c of the same length as x and y, even if the degree is less than the normal length(x)-1,
# in which case the array has some trailing zeroes.
# The polynomial is thus p(x) = c[1] + c[2]x + ... c[d+1] x^d where d =length(x)-1, the nominal degree.
n_nodes = length(x)
degree = n_nodes - 1
V = zeros(n_nodes,n_nodes)
for i in 0:degree
for j in 0:degree
V[i+1,j+1] = x[i+1]^j # Shift the array indices up by one, since Julia counts from 1, not 0.
end
end
c = solvelinearsystem(V, y)
return c
end;
Least squares polynomial approximation#
function polyfit(x, y, n)
# Version 2: least squares fitting.
# Compute the coeffients c_i of the polynomial of degree n that give the best least squares fit to data (x[i], y[i]).
N = length(x)
m = zeros(2n+1)
for k in 0:2n
m[k+1] = sum(x.^k) # Here and below, shift the indices up by one, since Julia counts from 1, not 0.
end
M = zeros(n+1,n+1)
for i in 0:n
for j in 0:n
M[i+1, j+1] = m[i+j+1]
end
end
p = zeros(n+1)
for k in 0:n
p[k+1] = sum(x.^k .* y)
end
c = solvelinearsystem(M, p)
return c
end;
Evaluate a polynomial#
function polyval(x; coeffs) # coeffs has to be a keyword argument in order that only x gets vectorized
# Evaluate the polynomial with coefficients in c (as given by polyfit, for example).
# If x is an array, the usage becomes y = polyval.(c, x)
# for each element of array x.
y = coeffs[1]
for i in 2:length(coeffs)
y += coeffs[i]*x^(i-1)
end
return y
end;
Derivatives and Definite Integrals#
Minimization#
Differential Equations#
Euler’s method#
function eulermethod(f, a, b, u_0; n=100)
# Solve du/dt = f(t, u) for t in [a, b], with initial value u(a) = u_0
h = (b-a)/n
t = range(a, b, n+1) # Note: "n" counts steps, so there are n+1 values for t.
u = zeros(n+1)
u[1] = u_0
for i in 1:n
u[i+1] = u[i] + f(t[i], u[i])*h
end
return (t, u)
end;
function eulermethod_errorcontrol(f, a, b, u_0; errortolerance=1e-3, h_min=1e-6, h_max=0.1, steps_max=1000, demomode=false)
steps = 0
t_i = a
U_i = u_0
t = [t_i]
U = [U_i]
h = h_max # Start optimistically!
while t_i < b && steps < steps_max
K_1 = h*f(t_i, U_i)
K_2 = h*f(t_i + h/2, U_i + K_1/2)
errorestimate = abs(K_1 - K_2)
s = 0.9 * sqrt(errortolerance/errorestimate)
if errorestimate <= errortolerance # Success!
t_i += h
U_i += K_1
append!(t, t_i)
append!(U, U_i)
# Adjust step size up, but not too big
h = min(s*h, h_max)
else # Innacurate; reduce step size and try again
h = max(s*h, h_min)
if demomode
println("t_i=$t_i: Decreasing step size to $(about(h)) and trying again.")
end
end
# A refinement not mentioned above; the next step should not overshoot t=b:
if t_i + h > b
h = b - t_i
end
steps += 1
end
return (t, U)
# Note: if the step count ran out, this does not reach t=b, but at least it is correct as far as it goes
end;
function eulermethod_system(f, a, b, u_0, n)
# TO DO: one could use multiple dispatch to keep the name "eulermethod".
# The conversion for the system version is mainly "U[i] -> U[i,:]"
h = (b-a)/n
t = range(a, b, n+1)
# The following three lines and the one in the for loop below change for the system version
n_unknowns = length(u_0)
U = zeros(n+1, n_unknowns)
U[1,:] = u_0 # Only for system version
for i in 1:n
U[i+1,:] = U[i,:] + f(t[i], U[i,:])*h # For the system version
end
return (t, U)
end;
The explicit trapezoid method#
function explicittrapezoid(f, a, b, u_0; n=100, demomode=false)
# Use the Explict Trapezoid Method (a.k.a Improved Euler) to solve
# du/dt = f(t, u)
# for t in [a, b], with initial value u(a) = u_0
h = (b-a)/n
t = range(a, b, n+1) # Note: "n" counts steps, so there are n+1 values for t.
u = zeros(n+1)
u[1] = u_0
for i in 1:n
K_1 = f(t[i], u[i])*h
K_2 = f(t[i]+h, u[i]+K_1)*h
u[i+1] = u[i] + (K_1 + K_2)/2.0
end;
return (t, u)
end;
function explicittrapezoid_system(f, a, b, u_0, n)
# Use the Explict Trapezoid Method (a.k.a Improved Euler) to solve the system
# du/dt = f(t, u) for t in [a, b], with initial value u(a) = u_0
# The conversion for the system version is mainly "u[i] -> u[i,:]"
h = (b-a)/n
t = range(a, b, n+1)
n_unknowns = length(u_0)
u = zeros(n+1, n_unknowns)
u[1,:] = u_0
for i in 1:n
K_1 = f(t[i], u[i,:])*h
K_2 = f(t[i]+h, u[i,:]+K_1)*h
u[i+1,:] = u[i,:] + (K_1 + K_2)/2.0
end
return (t, u)
end;
The explicit midpoint method#
function explicitmidpoint(f, a, b, u_0; n=100, demomode=false)
# Use the Explicit Midpoint Method (a.k.a Modified Euler) to solve
# du/dt = f(t, u) for t in [a, b], with initial value u(a) = u_0
h = (b-a)/n
t = range(a, b, n+1) # Note: "n" counts steps, so there are n+1 values for t.
u = zeros(length(t))
u[1] = u_0
for i in 1:n
K_1 = f(t[i], u[i])*h
K_2 = f(t[i]+h/2, u[i]+K_1/2)*h
u[i+1] = u[i] + K_2
end;
return (t, u)
end;
function explicitmidpoint_system(f, a, b, u_0, n)
# Use the Explict Midpoint Method (a.k.a Modified Euler) to solve the system
# du/dt = f(t, u) for t in [a, b], with initial value u(a) = u_0
# The conversion for the system version is mainly "u[i] -> u[i,:]"
h = (b-a)/n
t = range(a, b, n+1)
n_unknowns = length(u_0)
u = zeros(n+1, n_unknowns)
u[1,:] = u_0
for i in 1:n
K_1 = f(t[i], u[i,:])*h
K_2 = f(t[i]+h/2, u[i,:]+K_1/2)*h
u[i+1,:] = u[i,:] + K_2
end
return (t, u)
end;
The Runge-Kutta method#
function rungekutta(f, a, b, u_0; n=100, demomode=false)
# Use the (classical) Runge-Kutta Method to solve
# du/dt = f(t, u) for t in [a, b], with initial value u(a) = u_0
h = (b-a)/n
t = range(a, b, n+1) # Note: "n" counts steps, so there are n+1 values for t.
u = zeros(length(t))
u[1] = u_0
for i in 1:n
K_1 = f(t[i], u[i])*h
K_2 = f(t[i]+h/2, u[i]+K_1/2)*h
K_3 = f(t[i]+h/2, u[i]+K_2/2)*h
K_4 = f(t[i]+h, u[i]+K_3)*h
u[i+1] = u[i] + (K_1 + 2*K_2 + 2*K_3 + K_4)/6
end;
return (t, u)
end;
function rungekutta_system(f, a, b, u_0, n)
# Use the (classical) Runge-Kutta Method to solve
# du/dt = f(t, u) for t in [a, b], with initial value u(a) = u_0
# The conversion for the system version is mainly "u[i] -> u[i,:]"
h = (b-a)/n
t = range(a, b, n+1)
n_unknowns = length(u_0)
u = zeros(n+1, n_unknowns)
u[1,:] = u_0
for i in 1:n
K_1 = f(t[i], u[i,:])*h
K_2 = f(t[i]+h/2, u[i,:]+K_1/2)*h
K_3 = f(t[i]+h/2, u[i,:]+K_2/2)*h
K_4 = f(t[i]+h, u[i,:]+K_3)*h
u[i+1,:] = u[i,:] + (K_1 + 2*K_2 + 2*K_3 + K_4)/6
end
return (t, u)
end;
Some auxilliary functions#
For examples, presentation of results, etc.
Helper function printmatrix#
function printmatrix(A)
# A helper function to "pretty print" matrices
(rows, cols) = size(A)
print("[ ")
for row in 1:rows
if row > 1
print(" ")
end
for col in 1:cols
print(A[row,col], " ")
end
if row < rows;
println()
else
println("]")
end
end
end;
# A helper function for rounding some output to three significant digits
approx3(x) = round(x, sigdigits=3);
# A helper function for rounding some output to four significant digits
approx4(x) = round(x, sigdigits=4);
For some examples in Chapter Initial Value Problems for Ordinary Differential Equations#
f_mass_spring(t, u) = [ u[2], -(K/M)*u[1] - (D/M)*u[2] ];
function y_mass_spring(t; t_0, u_0, K, M, D)
(y_0, v_0) = u_0
discriminant = D^2 - 4K*M
if discriminant < 0 # underdamped
omega = sqrt(4K*M - D^2)/(2M)
A = y_0
B = (v_0 + y_0*D/(2M))/omega
return exp(-D/(2M)*(t-t_0)) * ( A*cos(omega*(t-t_0)) + B*sin(omega*(t-t_0)))
elseif discriminant > 0 # overdamped
Delta = sqrt(discriminant)
lambda_plus = (-D + Delta)/(2M)
lambda_minus = (-D - Delta)/(2M)
A = M*(v_0 - lambda_minus * y_0)/Delta
B = y_0 - A
return A*exp(lambda_plus*(t-t_0)) + B*exp(lambda_minus*(t-t_0))
else
lambda = -D/(2M)
A = y_0
B = v_0 - A * lambda
return (A + B*t)*exp(lambda*(t-t_0))
end
end;
function damping(K, M, D)
if D == 0
println("Undamped")
else
discriminant = D^2 - 4K*M
if discriminant < 0
println("Underdamped")
elseif discriminant > 0
println("Overdamped")
else
println("Critically damped")
end
end
end;
end;