PENG9570 Course Summary

PENG9570: Course summary

Questions and suggested answers

Mathematical modelling

What is a mathematical model?

• a description of a system using mathematical concepts and language.

Why do we need mathematical models?

• To develop scientific understanding
• To predict outcomes of experiments
• To test the effect of changes in a system
• To estimate parameters in a system

Describe the stages of the modelling cycle

ODEs and systems of ODEs

What is a differential equation and what is an initial value problem?

• Differential equation:

$$F\left(t,y,y^{\prime },\dots ,y^{(n)}\right)=0$$

• Initial value problem:
• Differential equation plus initial values for $y$ and, possibly, some of its derivatives
• Example:

$$y^{\prime }=-y,y(0)=1$$

What is the steady state (or critical point or equilibrium) of a first order differential equation?

• In one-variable models $u^{\prime }=f(u)$ :

The value $u^{\ast }$ such that the derivative is zero: $u^{\prime }=f\left(u^{\ast }\right)=0$

• In two-variable models $x^{\prime }=P(x,y),y^{\prime }=Q(x,y)$ :

The point $\left(x^{\ast },y^{\ast }\right)$ such that the derivatives are zero: $x^{\prime }=P\left(x^{\ast },y^{\ast }\right)=0$ and $y^{\prime }=Q\left(x^{\ast },y^{\ast }\right)=0$

What does it mean that a steady state is (i) stable, (ii) asymptotically stable and (iii) unstable?

• Assume $u^{\ast }$ is a steady state
• Stable: small perturbation from $u^{\ast }$ has little effect (solution stays near $u^{\ast }$ )
• Unstable: small perturbation from $u^{\ast }$ has large effect (solution escapes from $u^{\ast }$ )
• Asymptotic stability: a solution that starts near enough to $u^{\ast }$ will converge to $u^{\ast }$ as $t\to \infty$

How are phase line plots used to provide insights into the dynamics of the solution?

$$u^{\prime }=f(u)$$

$$u^{\prime }=u(1-u)-h$$

Two-dimensional dynamical systems

$$\begin{array}{rl}& x^{\prime }=P(x,y)\\ & y^{\prime }=Q(x,y)\end{array}$$

Solution: $x(t),y(t)$ - functions of time

Solution Curres:

What is the linearization of an ODE system?

In what way does the linearization inform us about the nonlinear system?

Assumption: critical point in $(0,0)$.

Linearisation: (due to $P(x,y)\approx P_x(0,0)x+P_y(0,0)y)$

$$\left[\begin{array}{l}{x}^{\prime }\\ y^{\prime }\end{array}\right]=\left[\begin{array}{ll}{P}_x(0,0)& P_y(0,0)\\ Q_x(0,0)& Q_y(0,0)\end{array}\right]\left[\begin{array}{l}x\\ y\end{array}\right]$$

$x^{\prime }=Ax:A$ is the Jacobian.

Dynamics can be deduced from eigenvalues and eigenvectors of Jacobian

Name the four types of critical points and relate them to the eigenvalues of the coefficient matrix (linear system) or of the Jacobian (nonlinear system).

Node

• Two real eigenvalues of same sign
• Stable if both eigenvalues are negative
• Unstable if both eigenvalues are positive

Spiral

• Complex eigenvalues
• Stable if $\mathrm{Re}\lambda <0$ for both eigenvalues
• Unstable if $\mathrm{Re}\lambda >0$ for both eigenvalues

Saddle

• Two real eigenvalues, one positive and one negative

Center

• Two purely imaginary eigenvalues

Node

• Two real eigenvalues of same sign
• Stable if both eigenvalues are negative
• Unstable if both eigenvalues are positive

Spiral

• Complex eigenvalues
• Stable if $\mathrm{Re}\lambda <0$ for both eigenvalues
• Unstable if $\mathrm{Re}\lambda >0$ for both eigenvalues

Saddle

• Two real eigenvalues, one positive and one negative

Center

• Two purely imaginary eigenvalues
• Two real eigenvalues of same sign
• Stable if both eigenvalues are negative
• Unstable if both eigenvalues are positive

Spiral

• Two real eigenvalues, one positive and one negative

Center

• Two purely imaginary eigenvalues
• Two real eigenvalues of same sign
• Stable if both eigenvalues are negative
• Unstable if both eigenvalues are positive

Spiral

• Complex eigenvalues
• Stable if $\mathrm{Re}\lambda <0$ for both eigenvalues
• Unstable if $\mathrm{Re}\lambda >0$ for both eigenvalues

Saddle

• Two real eigenvalues, one positive and one negative

Center

• Two purely imaginary eigenvalues

What is a nullcline?

Explain how nullclines and vector fields are used to provide information on solution orbits

The nullclines of a system are the curves obtained when the derivatives $x^{\prime }$ and $y^{\prime }$ are set equal to zero:

$$\begin{array}{rl}& x^{\prime }=0\iff P(x,y)=0\\ & y^{\prime }=0\iff Q(x,y)=0\end{array}$$

What is a bifurcation and what is depicted in a bifurcation diagram?

• Assume there is a parameter (the bifurcation parameter) in our ODE system
• We have a bifurcation point at the value of this parameter where a new type of solution emerges
• A bifurcation diagram shows the possible steady state values (stable and unstable) of a system as a function of a bifurcation parameter
• I want- [ ] I want to know what is being plotted here on the bifurcation diagram

Existence and uniqueness theorem for IVPs: Define the integral form and the Picard iteration

Integral form of IVP:

• [?] IVP?
  • [i] Initial Value Problem

$$y(t)=y\left(t_0\right)+{\int }_{t_0}^{t}f(\tau ,y(\tau ))d\tau$$

The Picard iteration

Sequence of approximations to integral form:

$$y_{n+1}(t)=y_n(t)+{\int }_{t_0}^{t}f\left(\tau ,y_n(\tau )\right)d\tau ,y_0(t)=y\left(t_0\right)$$

• [?] No idea what this is supposed to mean. It’s the same formula but with some variable changed, but not explained

If the sequence converges to a limit function $y(t),y(t)$ is the solution of the IVP.

Existence and uniqueness theorem for IVPs: Define Cauchy sequence and completeness

Cauchy sequence

A sequence where the elements come arbitrarily close to each other when the sequence numbers are sufficiently large. Mathematically, for every positive $ϵ$, there is an integer $N$ such that for all $m,n>N$,

$$\left|x_m-x_n\right|<ϵ.$$

Completeness

A set $S$ is said to be complete provided that all Cauchy sequences of elements in the set converge to an element in $S$.

Thus, establishing that a sequence of elements belonging to a complete set $S$ is Cauchy, we can conclude that the sequence will converge to a limit in $S$.

Formulate the existence and uniqueness theorem

Given the IVP

$$y^{\prime }=f(t,y),y\left(t_0\right)=y_0$$

Providing that $f$ is continuous and differentiable ( $f\in C_1[E]$ ),

• What does this mean? Specifically C_1 and E there exists $a>0$ such that there exists a unique solution in $[-a,a]$. (Existence and uniqueness)

Give a sketch of the proof of existence

• SHIT this looks very very important
• Construct sequence of approximate solutions using Picard iteration.
• Show that this sequence is Cauchy within the space of continuous functions $C_1[E]$ (depends on the differentiability of $f$ ).
• Since $C_1[E]$ is complete, then the sequence must converge to an element within $C_1[E]$.
• This limit is the solution and existence is established.

Numerical methods for IVPs

Why do we need numerical methods to solve differential equations?

1. In general it is not possible to find an analytical solution
2. The differential equations need to be translated into equations that can be solved on a computer

How do we write an initial value problem on integral form?

Integrating the initial value problem (IVP)

$$y^{\prime }=f(t,y),y\left(t_0\right)=y_0$$

from $t$ to $t_0$ we obtain

$$y=y_0+{\int }_{t_0}^{t}f(t,y)dt$$

The integral is in general evaluated numerically to obtain a numerical solution of the IVP.

Use numerical integration to explain the accuracy of Euler’s method

• Use numerical integration to explain the accuracy of Euler’s method

$$y=y_0+{\int }_{t_0}^{t}f(t,y)dt$$

The integral is evaluated numerically (here $h=t_{n+1}-t_n$ ):

• Euler’s method: Riemann left sum (Order $h^2$ - method is order $h$ )
  • Not clear what “order h” means

$$y_{n+1}=y_n+{\int }_{t_n}^{t_{n+1}}f(t,y)dt\approx y_n+hf\left(t_n,y_n\right)$$

Euler’s steps

I think this is a iterative method. I’ve seen “Euler Steps” mentioned elsewhere. As in a number of timesteps to be integrated over piecewise.

Use numerical integration to explain the accuracy of Heun’s method

• Use numerical integration to explain the accuracy of Heun’s method

$$y=y_0+{\int }_{t_0}^{t}f(t,y)dt$$

The integral is evaluated numerically (here $h=t_{n+1}-t_n$ ):

• Heun’s method: Trapezoidal rule (Order $h^3$ - method is order $h^2$ )

$$y_{n+1}=y_n+{\int }_{t_n}^{t_{n+1}}f(t,y)dt\approx y_n+\frac{h}{2}\left(f\left(t_n,y_n\right)+f\left(t_{n+1},y_n+hf\left(t_n,y_n\right)\right)\right)$$

Use numerical integration to explain the accuracy of Runge-Kutta’s 4th order method

• Use numerical integration to explain the accuracy of Runge-Kutta’s 4th order method

$$y=y_0+{\int }_{t_0}^{t}f(t,y)dt$$

The integral is evaluated numerically (here $h=t_{n+1}-t_n$ ):

• Runge-Kutta’s 4th order method: Simpson’s rule (Order $h^5$ method is order $h^4$ )
• Estimate for solution at time step $t_{j+1}$ is obtained by stepwise estimation of slopes at time steps $t_j$ and $t_{j+1/2}=t_j+\Delta t/2$

PDEs

What is the difference between ODEs and PDEs ?

ODEs: solution is function of one variable

PDEs: solution is function of two or more variables and partial derivatives of more than one of the variables are part of the equation

Examples:

$$\begin{array}{c}{u}_t=k{u}_{xx}+f(u,x,t),\\ \\ u_{xt}=xt\\ \\ u_x+a{u}_t=0\end{array}$$

How is the diffusion PDE derived from a conservation law?

$u(x,t)$ : Quantity of a substance $J(x)$ : flux of substance (or heat flux) $f(t,x,u)$ : production rate Conservation law in integral form is obtained by requiring the following (on interval $[a,b]$ ):

Rate of change of amount of substance $=$ what comes in - what goes out + rate of production

This, in turn, gives the conservation law in differential form :

$$u_t=-J_x+f(x,t,u)$$

Diffusion (and heat): $J=-D{u}_x$ (Fick’s law) $\Rightarrow J_x=-D{u}_{xx}$.

$$u_t=D{u}_{xx}+f(x,t,u)$$

What do we mean by linearity and superposition of solutions?

PDEs are often written using operators. Example:

$$Lu=f(x,t),Lu=u_t-k{u}_{xx}$$

The operator $L$ is linear when

1. $L(u+v)=L(u)+L(v)$
2. $L(cu)=cL(u)$

For homogeneous problems $Lu=0$, super-positions of solutions are also solutions:

1. $w=c_1{u}_1+c_2{u}_2+\cdots +c_n{u}_n+\cdots$
2. $w={\int }_{I}c(\xi )u(x,t,\xi )d\xi$

What is an initial boundary value problem

(IBVP) for a PDE?

Example: 1D diffusion PDE with initial and boundary conditions

$$\begin{array}{cc}{u}_t=k{u}_{xx},t>0,& 0<x<L\\ u(0,t)=a(t),& t\ge 0\\ u(L,t)=b(t),& t\ge 0\\ u(x,0)=g(x),& 0\le x\le L\end{array}$$

Prove uniqueness of the IBVP

Assume two solutions $u_1$ and $u_2$. Then $w=u_1-u_2$ solves

$$\begin{array}{c}{w}_t=k{w}_{xx},t>0,0<x<L\\ w(0,t)=0,w(L,t)=0,t\ge 0\\ w(x,0)=0,0\le x\le L\end{array}$$

• Define the energy integral $E(t)={\int }_0^L{w}^2(x,t)dx$
• Show (without too much detail) by differentiation and integration by parts that $w$ must be zero
• This proves that $w\equiv 0$
• Uniqueness is proven!

Numerical methods for IBVPS

How are partial derivatives approximated by finite difference expressions?

Typical approximations:

$$\begin{array}{c}{u}_{xx}\approx \frac{u(x+\Delta x,t)-2u(x,t)+u(x-\Delta x,t)}{\Delta x^2}\\ u_t\approx \frac{u(x,t+\Delta t)-u(x,t)}{\Delta t}\end{array}$$

Let $v_{i,j}\approx u\left(x_i,t_j\right)$. Then,

$$\begin{array}{c}{u}_{xx}\approx \frac{v_{i-1,j}-2v_{i,j}+v_{i+1,j}}{\Delta x^2},\\ u_t\approx \frac{v_{i,j+1}-v_{i,j}}{\Delta t}\end{array}$$

Explain the explicit scheme

PDE: $u_t=u_{xx}$. Let $v_{i,j}\approx u\left(x_i,t_j\right)$. Then,

$$\begin{array}{rl}& \frac{v_{i,j+1}-v_{i,j}}{\Delta t}=\frac{v_{i-1,j}-2v_{i,j}+v_{i+1,j}}{\Delta x^2}\iff \\ & v_{i,j+1}=v_{i,j}+\alpha \left(v_{i-1,j}-2v_{i,j}+v_{i+1,j}\right)\end{array}$$

That is,

$${\mathbf{v}}_{j+1}=\hat{A}{\mathbf{v}}_j$$

where

$$\hat{A}=I-\alpha \hat{B},\hat{B}=\left[\begin{array}{cccc}2& -1& 0& 0\\ -1& 2& \ddots & 0\\ 0& \ddots & \ddots & -1\\ 0& 0& -1& 2\end{array}\right],{\mathbf{v}}_j=\left[\begin{array}{c}{v}_{1,j}\\ v_{2,j}\\ ⋮\\ v_{n,j}\end{array}\right]$$

Derive the stability criterion in terms of $\alpha =\Delta t/\Delta x^2$ for the explicit scheme

Explicit scheme:

$${\mathbf{v}}_{j+1}=\hat{A}{\mathbf{v}}_j$$

where

$$\hat{A}=I-\alpha \hat{B},\hat{B}=\left[\begin{array}{cccc}2& -1& 0& 0\\ -1& 2& \ddots & 0\\ 0& \ddots & \ddots & -1\\ 0& 0& -1& 2\end{array}\right],{\mathbf{v}}_j=\left[\begin{array}{c}{v}_{1,j}\\ v_{2,j}\\ ⋮\\ v_{n,j}\end{array}\right].$$

Scheme is stable if eigenvalues of $\hat{A}$ are less than 1 in absolute value. Eigenvalues are $1-\alpha \mu$ ( $\mu$ are eigenvalues of $\hat{B},0<\mu <4$ ).

$$-1<1-\alpha \mu <1\iff \alpha \mu <2$$

$\alpha \mu <2$ must hold for all values of $\mu$ including $\mu =4$ :

$$4\alpha <2\iff \alpha <\frac{1}{2}\iff \frac{\Delta t}{\Delta x^2}<\frac{1}{2}$$

What is the difference between explicit and implicit schemes?

• Explicit scheme:

$$\underbrace{\frac{v_{i,j+1}-v_{i,j}}{\Delta t}}=\frac{v_{i-1,j}-2v_{i,j}+v_{i+1,j}}{\Delta x^2}\Rightarrow {\mathbf{v}}_{j+1}=\hat{A}{\mathbf{v}}_j$$

forward time derivative (forward Euler)

• Implicit scheme:

$$\underbrace{\frac{v_{i,j}-v_{i,j-1}}{\Delta t}}=\frac{v_{i-1,j}-2v_{i,j}+v_{i+1,j}}{\Delta x^2}\Rightarrow {\mathbf{v}}_{j-1}=\hat{A}{\mathbf{v}}_j$$

backward time derivative (backward Euler)

Show that the implicit scheme is stable for all values of $\alpha =\Delta t/\Delta x^2$

• Implicit scheme is

$${\mathbf{v}}_{j-1}=\hat{A}{\mathbf{v}}_j\iff {\mathbf{v}}_{j+1}={\hat{A}}^{-1}{\mathbf{v}}_j(\hat{A}=I+\alpha \hat{B})$$

• Eigenvalues of $\hat{A}$ are $1+\alpha \mu >1$.

$$\text{ ⇒ eigenvalues of }{\hat{A}}^{-1}\text{ are }(1+\alpha \mu {)}^{-1}<1.$$

• Scheme is stable!

What is the difference between the Crank-Nicolson scheme and the previous schemes?

Same as implicit method, except that $u_{xx}$ is approximated with the average of approximations at time steps $t_j$ and $t_{j-1}$ :

$$u_{xx}\left(x_i,t_j\right)\approx \frac{1}{2}\frac{v_{i-1,j}-2v_{i,j}+v_{i+1,j}}{\Delta x^2}+\frac{1}{2}\frac{v_{i-1,j-1}-2v_{i,j-1}+v_{i+1,j-1}}{\Delta x^2}$$

Advantage: the Crank-Nicolson scheme is more accurate than the other schemes when evaluating at points $\left(x,t^{\prime }\right)=\left(x,t-\frac{\Delta t}{2}\right)$ : error at his point is $O\left(\Delta t^2\right)+O\left(\Delta x^2\right)$

Show that the Crank-Nicolson scheme is stable

• Crank-Nicolson scheme is of the form

$$(2\mathrm{I}+\alpha \hat{B}){\mathbf{v}}_j=(2\mathrm{I}-\alpha \hat{B}){\mathbf{v}}_{j-1}\iff {\mathbf{v}}_j={\left[(2I+\alpha \hat{B}{)}^{-1}(2\mathrm{I}-\alpha \hat{B})\right]}^j{\mathbf{v}}_{j-1}$$

• Eigenvalues of $(2I+\alpha \hat{B}{)}^{-1}(2\mathrm{I}-\alpha \hat{B})$ are $(2+\alpha \mu {)}^{-1}(2-\alpha \mu )$.

⇒ eigenvalues are $<1$.

• Scheme is stable!

$$i\hslash \frac{\partial \Psi }{\partial t}=\left(-\frac{{\hslash }^2}{2m}{\nabla }^2+V(x,y,z,t)\right)\Psi$$

The Schrödinger equation

$$i\hslash \frac{\partial \Psi }{\partial t}=-\frac{{\hslash }^2}{2m}\frac{d^2\Psi }{d{x}^2}$$

• How can the numerical methods previously considered be applied to solve this PDE ?
  • Quite straightforward, with some adjustments (e.g. solutions of the Schrödinger equation are generally complex)
• How can this PDE be rewritten as a coupled system of ODEs?
  • Approximate ${\psi }_k(t)\approx \Psi \left(x_k,t\right)$ and solve the resulting large ODE system
  • Explain how the solution is interpreted as a probability density
  • The integral ${\int }_a^b|\Psi (x,t){|}^{2}dx$ is probability of a particle being located in $[a,b]$