About

Numerical Computing with Julia

Lesson 1: Solving Systems of Linear Equations

Objective: Understand direct and iterative methods for solving systems of linear equations and implement them in Julia.

Outline: 1. Introduction (10 mins) - Brief overview of linear systems - Applications of solving linear systems in various fields

2. Direct Methods (20 mins)

   • Gaussian Elimination

   • LU Decomposition

   • Example and Julia implementation:

     using LinearAlgebra
     
     A = [3 2; 1 2]
     b = [4; 2]
     x = A \ b  # Solving Ax = b

3. Iterative Methods (20 mins)

   • Jacobi Method

   • Gauss-Seidel Method

   • Example and Julia implementation:

     function jacobi(A, b, tol=1e-10, max_iter=1000)
         n = size(A, 1)
         x = zeros(n)
         x_old = copy(x)
         for _ in 1:max_iter
             for i in 1:n
                 x[i] = (b[i] - A[i, 1:i-1] * x_old[1:i-1] - A[i, i+1:end] * x_old[i+1:end]) / A[i, i]
             end
             if norm(x - x_old, Inf) < tol
                 return x
             end
             x_old = copy(x)
         end
         return x
     end

Lesson 2: Matrix Factorizations

Objective: Learn different matrix factorization techniques (LU, QR, Cholesky) and their implementations in Julia.

Outline: 1. Introduction (5 mins) - Importance of matrix factorizations

2. LU Decomposition (15 mins)

   • Definition and applications

   • Example and Julia implementation:

     using LinearAlgebra
     
     A = [4 3; 6 3]
     L, U = lu(A)

3. QR Decomposition (15 mins)

   • Definition and applications

   • Example and Julia implementation:

     using LinearAlgebra
     
     A = [1 2; 3 4; 5 6]
     Q, R = qr(A)

4. Cholesky Decomposition (15 mins)

   • Definition and applications

   • Example and Julia implementation:

     using LinearAlgebra
     
     A = [4 2; 2 3]
     L = cholesky(A).L

Lesson 3: Eigenvalues and Eigenvectors

Objective: Understand the concepts of eigenvalues and eigenvectors, and learn how to compute them in Julia.

Outline: 1. Introduction (10 mins) - Definition of eigenvalues and eigenvectors - Applications in science and engineering

2. Computing Eigenvalues and Eigenvectors (20 mins)

   • Power Method

   • Example and Julia implementation:

     using LinearAlgebra
     
     A = [4 1; 2 3]
     λ, v = eigen(A)

3. Applications (20 mins)
   • Principal Component Analysis (PCA)
   • Vibrations and stability analysis
   • Example and discussion

Lesson 4: Implementing Linear Algebra Algorithms in Julia

Objective: Practice implementing various linear algebra algorithms in Julia and understand their practical applications.

Outline: 1. Introduction (5 mins) - Review of previous lessons

2. Practice Implementations (25 mins)
   • Gaussian Elimination
   • Jacobi Method
   • Eigenvalue computation
   • Example implementations
3. Applications (20 mins)

   • Solving real-world problems using linear algebra

   • Case studies and examples

   • Julia code implementation:

     using LinearAlgebra
     
     # Example: Solving a system of equations for a real-world problem
     A = [3 2 -1; 2 -2 4; -1 0.5 -1]
     b = [1; -2; 0]
     x = A \ b

Learning Objectives:

• Understand and solve Initial Value Problems (IVPs)
• Apply Euler’s method and improved Euler’s method
• Implement Runge-Kutta methods
• Implement ODE solvers in Julia

Lesson Plan

Lecture 1: Introduction to Initial Value Problems (IVPs)

1. Concepts Covered:
   • Definition and formulation of Initial Value Problems
   • Examples of real-world IVPs
   • Basic numerical integration concepts
2. Julia Functions:
   • Basic plotting functions from Plots.jl: plot, plot!
   • Functions to define differential equations: f(t, y) = ...
3. Activities:
   • Discuss and derive the mathematical formulation of IVPs
   • Illustrate with examples and graphical representations
   • Assign a simple IVP to solve analytically and numerically

Lecture 2: Euler’s Method and Improved Euler’s Method

1. Concepts Covered:
   • Euler’s method: Step-by-step numerical solution of ODEs
   • Improved Euler’s method (Heun’s method): Enhancement and accuracy
2. Julia Functions:
   • Loops: for, while
   • Function definition: function euler_method(f, t0, y0, h, n)
   • Basic arithmetic operations and plotting
3. Activities:
   • Derive Euler’s method formula
   • Implement Euler’s method in Julia
   • Compare analytical solutions with Euler’s method results
   • Introduce and implement Improved Euler’s method

Lecture 3: Runge-Kutta Methods

1. Concepts Covered:
   • Overview of Runge-Kutta methods
   • Derivation and implementation of the classical 4th-order Runge-Kutta method (RK4)
2. Julia Functions:
   • Advanced function definitions: function rk4(f, t0, y0, h, n)
   • Plotting solutions for comparison: plot(sol_analytical, sol_rk4)
3. Activities:
   • Derive the RK4 method step by step
   • Implement RK4 in Julia
   • Compare results from Euler’s method, Improved Euler’s method, and RK4

Lecture 4: Implementing ODE Solvers in Julia

1. Concepts Covered:
   • Utilizing Julia packages for solving ODEs
   • Introduction to DifferentialEquations.jl
2. Julia Functions:
   • Using DifferentialEquations.jl: @ode_def, solve, ODEProblem, plot from Plots.jl
3. Activities:
   • Install and import DifferentialEquations.jl package
   • Define an ODE problem using DifferentialEquations.jl
   • Solve and plot the solution using built-in solvers

Useful Julia Functions

• Differential Equations: DifferentialEquations.jl, @ode_def
• Plotting: Plots.jl, plot, plot!
• Basic Julia: function, for, while, basic arithmetic operators

1. Classification of PDEs

Partial Differential Equations (PDEs) can be classified based on their linearity, order, and the nature of the coefficients. Here are some common classifications:

• Linearity: Linear vs. Nonlinear
• Order: First-order, Second-order, etc.
• Nature of Coefficients: Constant coefficients vs. Variable coefficients

In this tutorial, we’ll focus on linear second-order PDEs with constant coefficients.

2. Finite Difference Methods for PDEs

Finite difference methods are numerical techniques for solving PDEs by approximating derivatives with difference equations. Here’s a simple example using Julia:

# Define grid size and time step
nx, nt = 100, 100
dx, dt = 1.0/nx, 0.01

# Initialize solution matrix
u = zeros(nx, nt)

# Initial condition
u[:, 1] = sin.(2 * π * (0:dx:1))

# Finite difference loop
for n in 1:nt-1
    for i in 2:nx-1
        u[i, n+1] = u[i, n] + dt * (u[i+1, n] - 2*u[i, n] + u[i-1, n]) / dx^2
    end
end

3. Solving Heat and Wave Equations

Let’s tackle the heat equation and the wave equation using finite difference methods.

Heat Equation:

The heat equation is a PDE that describes the distribution of heat over time. The equation is:

\[\frac{\partial u}{\partial t} = \alpha \frac{\partial^2 u}{\partial x^2}\]

Here’s how to solve it using Julia:

# Parameters
α = 0.01
dx, dt = 0.1, 0.001
nx, nt = 100, 1000

# Initialize solution matrix
u = zeros(nx, nt)

# Initial condition
u[:, 1] = sin.(2 * π * (0:dx:1))

# Finite difference loop
for n in 1:nt-1
    for i in 2:nx-1
        u[i, n+1] = u[i, n] + α * dt * (u[i+1, n] - 2*u[i, n] + u[i-1, n]) / dx^2
    end
end

Wave Equation:

The wave equation describes the propagation of waves. The equation is:

\[\frac{\partial^2 u}{\partial t^2} = c^2 \frac{\partial^2 u}{\partial x^2}\]

Here’s how to solve it using Julia:

# Parameters
c = 1.0
dx, dt = 0.1, 0.01
nx, nt = 100, 100

# Initialize solution matrices
u = zeros(nx, nt)
u_new = copy(u)
u_old = copy(u)

# Initial condition
u[:, 1] = sin.(2 * π * (0:dx:1))

# Finite difference loop
for n in 2:nt-1
    for i in 2:nx-1
        u_new[i, n+1] = 2*u[i, n] - u_old[i, n] + (c^2 * dt^2 * (u[i+1, n] - 2*u[i, n] + u[i-1, n]) / dx^2)
    end
    u_old = copy(u)
    u = copy(u_new)
end

4. Implementing PDE Solvers in Julia

To make the PDE solvers more reusable and modular, we can wrap them into functions.

function heat_equation_solver(α, dx, dt, nx, nt, initial_condition)
    u = zeros(nx, nt)
    u[:, 1] = initial_condition

    for n in 1:nt-1
        for i in 2:nx-1
            u[i, n+1] = u[i, n] + α * dt * (u[i+1, n] - 2*u[i, n] + u[i-1, n]) / dx^2
        end
    end
    return u
end

function wave_equation_solver(c, dx, dt, nx, nt, initial_condition)
    u = zeros(nx, nt)
    u_new = copy(u)
    u_old = copy(u)
    u[:, 1] = initial_condition

    for n in 2:nt-1
        for i in 2:nx-1
            u_new[i, n+1] = 2*u[i, n] - u_old[i, n] + (c^2 * dt^2 * (u[i+1, n] - 2*u[i, n] + u[i-1, n]) / dx^2)
        end
        u_old = copy(u)
        u = copy(u_new)
    end
    return u
end

# Example usage:
initial_condition = sin.(2 * π * (0:dx:1))
u_heat = heat_equation_solver(0.01, 0.1, 0.001, 100, 1000, initial_condition)
u_wave = wave_equation_solver(1.0, 0.1, 0.01, 100, 100, initial_condition)

That’s a brief overview of solving PDEs using finite difference methods in Julia.