2022-09-21 Transient Intro#

Last time#

  • Structured by-hand differentiation

  • Concept of PDE-based inference (inverse problems)

  • The Blasius problem (activity)

Today#

  • Blasius working session

  • Intro to transient problems

  • Stability diagrams

  • Energy dissipation

  • \(A\) and \(L\) stability

using Plots
default(linewidth=3)
using LinearAlgebra
using SparseArrays
using Zygote

function vander(x, k=nothing)
    if k === nothing
        k = length(x)
    end
    V = ones(length(x), k)
    for j = 2:k
        V[:, j] = V[:, j-1] .* x
    end
    V
end
vander (generic function with 2 methods)

Compressible Blasius boundary layer#

  • Activity will solve this 1D nonlinear PDE

Ordinary Differential Equations#

Given initial condition \(u_0 = u(t=0)\), find \(u(t)\) for \(t > 0\) that satisfies

\[ \dot u \equiv \frac{\partial u}{\partial t} = f(t, u) \]

Application

\(u\)

\(f\)

Orbital dynamics

position, momentum

conservation of momentum

Chemical reactions

concentration

conservation of atoms

Epidemiology

infected/recovered population

transmission and recovery

Heat transfer

temperature

conservation of energy

Seismology

displacement, momentum

conservative of momentum

Solving differential equations#

Linear equations#

\[ \dot u = A(t) u + \text{source}(t)\]
  • Autonomous if \(A(t) = A\) and source independent of \(t\)

  • Suppose \(u\) and \(a = A\) are scalars: \(u(t) = e^{at} u_0\)

Can do the same for systems#

\[ y(t) = e^{A t} y_0 \]

What does it mean to exponentiate a matrix?#

Taylor series!

\[ e^A = 1 + A + \frac{A^2}{2} + \frac{A^3}{3!} + \dotsb \]

and there are many practical ways to compute it.

Question#

Suppose that the diagonalization \(A = X \Lambda X^{-1}\) exists and derive a finite expression for the matrix exponential using the scalar exp function.

Forward Euler method#

function ode_euler(f, u0; tfinal=10., h=0.1)
    u = copy(u0)
    t = 0.
    thist = [t]
    uhist = [u0]
    while t < tfinal
        tnext = min(t+h, tfinal)
        h = tnext - t
        u += h * f(t, u)
        t = tnext
        push!(thist, t)
        push!(uhist, u)
    end
    thist, hcat(uhist...)
end
ode_euler (generic function with 1 method)
f1(t, u; k=11) = -k * (u .- cos(t))

thist, uhist = ode_euler(f1, [.5], tfinal=10, h=.2)
scatter(thist, uhist[1,:])
plot!(cos)
../_images/2022-09-21-transient_12_0.svg