Julia Functions & Processes for Numerical Methods

Overview

This page provides a collection of Julia code examples using only built-in functions. These examples can help you solve problems such as LU Decomposition, root finding using Newton's Method, computing matrix norms, approximating functions with Taylor series, and computing Jacobians for systems of equations.

LU Decomposition (2×2 Example)

The following code shows a manual LU Decomposition for a 2×2 matrix and uses it to solve a linear system \( Ax = b \):

# LU Decomposition for a 2×2 matrix
function lu2x2(A)
    a11, a12 = A[1,1], A[1,2]
    a21, a22 = A[2,1], A[2,2]
    L = [1.0 0.0; a21/a11 1.0]
    U = [a11 a12; 0.0 a22 - (a21/a11)*a12]
    return L, U
end

# Solve Ax = b using the computed L and U
function solve_lu2x2(A, b)
    L, U = lu2x2(A)
    # Forward substitution: L*y = b
    y1 = b[1]
    y2 = b[2] - L[2,1] * y1
    # Backward substitution: U*x = y
    x2 = y2 / U[2,2]
    x1 = (y1 - U[1,2] * x2) / U[1,1]
    return [x1, x2], L, U
end

# Example usage:
A = [4.0 2.0; 6.0 3.0]
b = [10.0, 15.0]
x, L, U = solve_lu2x2(A, b)
println("Solution x: ", x)
      

Newton's Method for Root Finding

This code implements Newton's Method to find a root of a function \( f(x) = 0 \) using a provided derivative \( f'(x) \):

# Newton's Method implementation
function newtons_method(f, df, x0, tol=1e-6, maxiter=100)
    x = x0
    for i in 1:maxiter
        fx = f(x)
        if abs(fx) < tol
            return x, i
        end
        x -= fx / df(x)
    end
    error("Newton's method did not converge within the maximum iterations")
end

# Example: finding sqrt(2)
f(x) = x^2 - 2
df(x) = 2x
root, iterations = newtons_method(f, df, 1.5)
println("Root: ", root, " found in ", iterations, " iterations")
      

Computing Matrix Norms

A simple function to compute the Frobenius norm of a matrix:

# Frobenius norm calculation
function frobenius_norm(A)
    return sqrt(sum(abs2, A))
end

# Example usage:
A = [1.0 2.0; 3.0 4.0]
println("Frobenius norm of A: ", frobenius_norm(A))
      

Taylor Series Expansion (Example)

This example computes a Taylor polynomial approximation for \( f(x) = e^x \) around a point \( a \). (For a general function, you would compute derivatives manually.)

# Taylor polynomial for f(x) = e^x around point a
function taylor_exp_e(a, n, x)
    sum = 0.0
    for k in 0:n
        sum += exp(a) / factorial(k) * (x - a)^k
    end
    return sum
end

# Example usage:
a = 0.0
n = 5
x = 1.0
println("Taylor approximation of e^x at x=", x, ": ", taylor_exp_e(a, n, x))
      

Numerical Jacobian for Systems

This function computes the Jacobian matrix of a vector-valued function \( F(x) \) using finite differences:

# Numerical Jacobian using finite differences
function numerical_jacobian(F, x, h=1e-8)
    n = length(x)
    J = zeros(n, n)
    F0 = F(x)
    for j in 1:n
        xh = copy(x)
        xh[j] += h
        Fh = F(xh)
        J[:, j] = (Fh - F0) / h
    end
    return J
end

# Example: Jacobian for a system of two functions
F(x) = [x[1]^2 + x[2]; sin(x[1]) + cos(x[2])]
x0 = [1.0, 1.0]
println("Jacobian at x0:")
println(numerical_jacobian(F, x0))
      

Notes on Plotting in Julia

While Julia does not include a built-in plotting library, if external packages are allowed you can use the Plots package. For example:

# Uncomment the following lines if you have the Plots package installed.
# using Plots
# x = 0:0.1:2π
# y = sin.(x)
# plot(x, y, label="sin(x)", title="Sine Function", xlabel="x", ylabel="sin(x)")
      

If you are restricted to built-in functions only, consider exporting computed data and visualizing it with another tool.