FNC Sections 11.4 to 11.5

Author

Andrew Pua

Published

February 7, 2024

Loading required package: JuliaCall

Julia version 1.9.2 at location /home/apua/julia-1.9.2/bin will be used.

Loading setup script for JuliaCall...

Finish loading setup script for JuliaCall.

Learning objectives

Describe stiffness is and what to pay attention to when solving systems of differential equations numerically
Move on from periodic end conditions to more general boundary conditions

Stiffness

“Problems for which the time step is dictated by stability rather than accuracy are referred to as stiff.”
“Stiffness is not a binary condition but a spectrum. It can arise in nonlinear problems and in problems having nothing to do with diffusion. Except for the mildest instances of stiffness, an implicit time stepping method is the best choice.”

More general PDEs

Consider more generally \[\mathbf{u}'=\mathbf{f}(t,\mathbf{u})\]
- To use Section 11.3 and talk about absolute stability, we need to linearize around an exact solution to produce something similar to a first-order system of linear differential equations.
- But the Jacobian depends on time, so the resulting system does not have constant coefficients. Solution is to “freeze time”. Then pretty much analyze absolute stability in the way Section 11.3 does.
Also have to be mindful of magnitudes of the eigenvalues. Strikingly different magnitudes leads us to consider multiple time scales.
“Hence it is desirable in stiff problems generally, and diffusion problems in particular, to have a stability region that is unbounded in at least the negative real direction.”

Exercises

Problem 2

This exercise is about the IVP \(u'=\cos(t) - 200(u-\sin(t))\), \(u(0)=0\).

Show that \(u(t) = \sin(t)\) is the exact solution, and find the linearization about this solution.
- The linearization is given by \(v^\prime (t)=-200 v(t)\).
Find the lone eigenvalue of the Jacobian. What other time scale is also relevant in the solution?
- The lone eigenvalue is \(-200\). The relevant time scale is \(1/200\).
Use FNC.ab4 to solve the IVP over \(t\in[0,\pi/2]\) with \(n=800,850,900,\ldots,1200\) steps. By comparing to the exact solution, show that the method gets either no accurate digits or at least 11 accurate digits.

ivp = ODEProblem((u,p,t)-> cos(t)-200*(u-sin(t)),0.,(0.,pi/2));
n = 800:50:1200;
err = [];
for n in n
    t,u = FNC.ab4(ivp,n);
    push!( err, norm(sin.(t)-u,Inf) );
end
pretty_table([n err],["n","inf-norm error"])

Problem 4

The equation \(u'=u^2-u^3\) is a simple model for combustion of a flame ball in microgravity. After “ignition,” the exact solution rapidly approaches 1.

Solve the problem with initial condition \(u(0) =0.001\) for \(0\le t \le 2000\), using FNC.rk4 with \(n=2000\) steps. Plot the solution.
Find the \(1\times 1\) Jacobian of this system, and use it with the absolute stability figure for RK4 to derive an upper bound on the time step of RK4 when \(u=1\).
- The Jacobian is equal to \(-1\) because \(2u-3u^2\) is evaluated at \(u=1\).
- A rough upper bound for the time step is \(\tau\leq 3\).
Repeat the computations with \(n=200,300,\ldots,1000\), and make a table of the error at the final time, assuming the exact solution is 1. Discuss.

# (a)
ivp = ODEProblem((u,p,t)-> u^2-u^3,0.001,(0.,2000.));
t, u = FNC.rk4(ivp,2000);
plot(u)
# (c)
n = 200:100:1000
err = []
for n in n
    t,u = FNC.rk4(ivp,n)
    push!( err, u[n]-1 )
end
pretty_table([n err],["n","error at final time"])

Problem 5

The van der Pol equation is a famous nonlinear oscillator given by \[y'' - \mu(1-y^2)y' + y = 0,\] where \(\mu\ge 1\) is a constant.

Write the equation as a first-order system and find its Jacobian.
- Let \(u=y\) and \(v=y^\prime\). We can rewrite the equation as a system of first-order differential equations: \[\begin{eqnarray}u^\prime=v \\ v^\prime=\mu(1-u^2)*v-u\end{eqnarray}\]
- Matrix of first-order partial derivatives of \(f\) with respect to \(u\) and \(v\): \[\begin{pmatrix}0 & 1 \\ -2u\mu v-1 & \mu(1-u^2)\end{pmatrix}\]
Find the eigenvalues of the Jacobian when \(y=\pm 2\) and \(y'=0\).
- Here \(u=\pm 2\) and \(v=0\). So, \[\begin{pmatrix}0 & 1 \\ -1 & -3\mu\end{pmatrix}\]
Solve the problem solve with \(\mu=500\), \(y(0) = y'(0) = 1\), for \(0\le t \le 2000\). Plot \(y(t)\) as a function of \(t\).
Define \(M(t)\) as the minimum (i.e., most negative) real part of the eigenvalues of the Jacobian using the computed solution at time \(t\). Evaluate \(M\) for each \(t=0,2,4,6,\ldots,2000\), and plot \(M(t)\). Explain how your plot relates to the prior two findings.

function ode(u,p,t)
    f = [ 
        u[2],
        500*(1-u[1]^2*u[2])-u[1]
        ]
    return f
end
ivp = ODEProblem(ode,[1., 1.],(0.,2000.));
sol=solve(ivp);
plot(sol)
t,u = sol.t[1:2:end],sol.u[1:2:end];
λ = fill(0.0im,length(t),2);
for (k,u) in enumerate(u)
    J = [0 1; -2*u[1]*500*u[2]-1 500*(1-u[1]^2)];
    λ[k,:] .= eigvals(J);
end
plot(t, minimum(real(λ), dims=2))

Including boundary conditions

Instead of periodic end conditions from Section 11.2, which effectively “controlled”
Setup:
- Equation: \[u_t = \phi(t,x,u,u_x,u_{xx}), \quad a \le x \le b\]
- Boundary condition at endpoint \(a\): \[g_1\left( u(a,t), \frac{\partial u}{\partial x}(a,t) \right) = 0\]
- Boundary condition at endpoint \(b\): \[g_2\left( u(b,t), \frac{\partial u}{\partial x}(b,t) \right) = 0\]
Numerical setup:
- \(m\) is number of nodes in space
- Discretize in space, but this time include \(u(a,t)=u_0\) and \(u(b,t)=u_m\).
- Let \(\mathbf{v}\) be solution values over the interior of the interval \((u_1,\ldots,u_{m-1})\). Since \(u_0\) and \(u_m\) are unknown, we need to recover them from the system \[\begin{eqnarray} g_1\left( u(a,t), \frac{\partial u}{\partial x}(a,t) \right) = 0 \Rightarrow g_1\left(u_0, u_0^\prime\right)=0 \\ g_2\left( u(b,t), \frac{\partial u}{\partial x}(b,t) \right) = 0 \Rightarrow g_2\left(u_m, u_m^\prime\right)=0 \end{eqnarray}\]
- Form \[\mathbf{u} = \begin{bmatrix} u_0 \\ \mathbf{v} \\ u_m \end{bmatrix}\]
- \(\mathbf{v}\) proceeds as usual like \(\mathbf{u}\) without the boundary conditions in Section 11.2.

Implementation

Refer to Algorithm 11.5.3 and Function 11.5.4 for details.
Chebyshev discretization is used for the matrix of derivatives \(\mathbf{D}_x\) and \(\mathbf{D}_{xx}\).
Recovering \(u_0\) and \(u_m\) may require solving a system of nonlinear equations. But there are special cases.
Section 11.5 Problem 7 require modifying the code to accommodate homogeneous Neumann boundary conditions. Only solving a linear system is required here.
- The idea is to include \(u^\prime(a)=0\) and \(u^\prime(b)=0\).
- System of equations for \(u_0\) and \(u_m\): \[\begin{eqnarray}\frac{-\frac{3}{2}u_0+2u_1-\frac{1}{2}u_2}{h}=0 \\ \frac{\frac{1}{2}u_{m-2}-2u_{m-1}+\frac{3}{2}u_m}{h}=0 \end{eqnarray}\]

function parabolicHN(ϕ,xspan,m,tspan,init)
    x,Dₓ,Dₓₓ = FNC.diffcheb(m,xspan) # include FNC prefix here
    int = 2:m    # indexes of interior nodes
    function extend(v)
        # linear system of equations
        ubc = -[-3/2 0 ; 0 3/2]\([2 -1/2 zeros(Float64, 1, m-3); zeros(Float64, 1, m-3) 1/2 -2]*v) 
        return [ubc[1];v;ubc[2]]
    end
    function ode!(f,v,p,t)
        u = extend(v)
        uₓ,uₓₓ = Dₓ*u,Dₓₓ*u
        @. f = ϕ(t,x[int],u[int],uₓ[int],uₓₓ[int])
    end
    ivp = ODEProblem(ode!,init.(x[int]),float.(tspan))
    u = solve(ivp)
    return x,t->extend(u(t))
end

Section 11.5 Problem 6 require modifying the code to accommodate Dirichlet boundary conditions.
- The idea is to include \(u(a,t)=\alpha\) and \(u(b,t)=\beta\).
- Solving any system of equations is not needed here.

function parabolicD(ϕ,xspan,m,alpha,beta,tspan,init)
    x,Dₓ,Dₓₓ = FNC.diffcheb(m,xspan) # include FNC prefix
    int = 2:m    # indexes of interior nodes
    function extend(v)
        return [alpha;v;beta]
    end
    function ode!(f,v,p,t)
        u = extend(v)
        uₓ,uₓₓ = Dₓ*u,Dₓₓ*u
        @. f = ϕ(t,x[int],u[int],uₓ[int],uₓₓ[int])
    end
    ivp = ODEProblem(ode!,init.(x[int]),float.(tspan))
    u = solve(ivp)
    return x,t->extend(u(t))
end

Examples

If boundary conditions involve first derivatives, replace with forward second-order finite differences if at the left endpoint. Replace with backward second-order finite differences if at the right endpoint.
Refer to Example 11.5.1 for an example with a backward finite differences.
Section 11.5 Problem 1 require both forward and backward finite differences.
Useful exercise to write down mathematical expressions for the PDEs with initial and boundary conditions for Demos 11.5.5 to 11.5.7.

Exercises

Problem 2

Use FNC.parabolic to solve the heat equation for \(0\le x \le 5\) with initial condition \(u(x,0)=x(5-x)\) and subject to the boundary conditions \(u(0,t)=0\), \(u(5,t)-u_x(5,t)=5\). Plot the solution at \(t=1\) and find the value of \(u(2.5,1)\).

ϕ = (t,x,u,uₓ,uₓₓ) -> uₓₓ;
init = x -> x*(5-x);
g₁ = (u,uₓ) -> u;
g₂ = (u,uₓ) -> u-uₓ-5;
x,u = FNC.parabolic(ϕ,(0,5),60,g₁,g₂,(0,1),init)
plot(x,u(1), xaxis=(L"x"), yaxis=(L"u(x,1)",(0,10)))

Problem 3

Use separation of variables to solve the IVP \(u_t=u^2\), \(u(0) = A>0\). What happens as \(t\to 1/A\) from below?
- Refer to this exposition about the idea behind separation of variables.
- The idea is to assume that a solution will have the form \(u(x,t)=F(x)G(t)\).
- \(u_t=u^2\) implies \[\begin{eqnarray} F(x)G^\prime(t)=\left(F(x)G(t)\right)^2 \\ \Rightarrow \frac{G^\prime(t)}{\left(G(t)\right)^2}= F(x)\end{eqnarray}\]
- The only way for this equation to hold is when both sides are equal to some constant which does not depend on \(t\) or \(x\). Call this constant \(\theta\neq 0\).
- Let \(C\) be some arbitrary constant. Solving the ODE \[\frac{G^\prime(t)}{\left(G(t)\right)^2}=\theta \Rightarrow G(t)=-\frac{1}{\theta t +C}\]
- So, we have \[u(x,t)=-\frac{\theta}{\theta t +C}\]
- The initial condition \(u(0)=A>0\) means that \[A=u(x,0)=-\frac{\theta}{C}\Rightarrow C=-\frac{\theta}{A}\]
- As a result, \[u(x,t)=-\frac{1}{t-1/A}\]
- Therefore, as \(t\to 1/A\) from below, \(u(x,t)\to\infty\).
Try to continue the solution in Demo 11.5.6 to \(t=1\). What happens?
- Note that here \(u_t=u^2+u_{xx}\).

# Code directly copied from Demo 11.5.6
ϕ = (t,x,u,uₓ,uₓₓ) -> u^2 + uₓₓ
g₁ = (u,uₓ) -> u
g₂ = (u,uₓ) -> uₓ
init = x -> 400x^4*(1-x)^2
# Modification next line
x,u = FNC.parabolic(ϕ,(0,1),60,g₁,g₂,(0,1),init);
anim = @animate for t in range(0,1,length=101) 
    plot(x,u(t), xaxis=(L"x"), yaxis=(L"u(x,t)",(0,10)),dpi=100, 
        title="Heat equation with source",leg=:topleft)
end
# change FPS
mp4(anim,"boundaries-source.mp4",fps=10)

Let the initial condition be \(u(x,0) = C x^4(1-x)^2\); the demo uses \(C=400\). To the nearest 10, find a critical value \(C_0\) such that the solution approaches zero asymptotically if \(C < C_0\), but not otherwise.
- Crude experimentation to be demonstrated.

Problem 6

Repeat Demo 11.5.5 with modified code allowing for homogeneous Dirichlet boundary conditions.

ϕ = (t,x,u,uₓ,uₓₓ) -> uₓₓ
init = x -> 1 + sin(π*x/2) + 3*(1-x^2)*exp(-4x^2);
x,u = parabolicD(ϕ,(-1,1),60,0,2,(0,0.75),init)
plt = plot(xlabel=L"x",ylabel=L"u(x,t)",legend=:topleft,
          title="Solution of the heat equation")
for t in 0:0.1:0.4
    plot!(x,u(t),label=@sprintf("t=%.2f",t))
end
plt

Problem 7

Solve heat equation \(u_t=u_{xx}\) on \(x \in [0,4]\), \(t\in [0,1]\) with initial condition \(u(x,0)=x^2(4-x)^4\) and homogeneous Neumann boundary conditions.

ϕ = (t,x,u,uₓ,uₓₓ) -> uₓₓ
init = x -> x^2*(4-x)^4;
x,u = parabolicHN(ϕ,(0,4),60,(0,1),init)
plt = plot(xlabel=L"x",ylabel=L"u(x,t)",legend=:topleft,
          title="Solution of the heat equation")
for t in 0:0.05:1
    plot!(x,u(t),label=@sprintf("t=%.2f",t))
end
plt