\documentclass[11pt]{scrartcl}

% vim: tw=50

\usepackage{notesheader}

\begin{document}

\title{Differential Equations}
\author{}
\date{\today}
\maketitle

\tableofcontents

\newpage
\setcounter{section}{-1}
\section{Introduction}

Colours for notes:
\begin{itemize}
    \item Text will be in black
    \item Displayed maths will be in blue
    \item Comments will be in green
    \item Examples will be in red
\end{itemize}

This is an \emph{applied} course --- we will
emphasise methods rather than proofs.
\newpage

\section{Basic Calculus}

\setcounter{subsection}{0}

\subsection{Differentiation}

\begin{definition*}[Derivative of a function]
    Define the \emph{derivative} of a function
    $f(x)$ with respect to $x$ as a function such
    that
    \[ \left. \dfrac{f}{x}\right|_{x_0}
    = \lim_{h \to 0} \frac{f(x_0 + h) - f(x)}{h} \]
\end{definition*}

\noindent The derivative can sort of be thought of
as the ``slope''. \\
For the derivative to exist at $x_0$, the
left-hand and right-hand limits must exist and be
equal, i.e.
\[ \lim_{h \to O^-} \frac{f(x_0 + h) - f(x)}{h}
= \lim_{h \to 0^+} \frac{f(x_0 + h) - f(x)}{h} \]

\begin{example*}[$f(x) = |x|$ at $x = 0$]
    Then
    \[ \lim_{h \to 0^-} = \frac{|h|}{h} = -1
    \qquad \text{but} \qquad \lim_{h \to 0^+}
    \frac{|h|}{h} = +1 \]
    so $|x|$ is not differentiable at $x = 0$ (but
    it is everywhere else).
\end{example*}

\begin{notation*}[Ways of writing a derivative]
    \[ \dfrac{f}{x} = f'(x) = \dot{f}(x) \]
    and we denote the $n$-th derivative by
    $f^{(n)}(x)$.
\end{notation*}

\subsection{$O$ and $o$ notation}

These are useful concepts to give comparative
scalings of a function as they approach some
limiting point

\begin{definition*}[$O$ and $o$ notation]
    $f(x)$ is $o(g(x))$ as $x \to x_0$ if
    \[ \lim_{x \to x_0} \frac{f(x)}{g(x)} = 0 .\]
    This is often written as $f(x) = o(g(x))$ as
    $x \to x_0$. Informally you can think of this
    as ``$f(x)$ is much smaller than $g(x)$ as
    $x \to x_0$.'' \\
    $f(x)$ is $O(g(x)$ as $x \to x_0$ if
    $\exists \delta > 0$ and $M > 0$ such that for
    all $x$ with $0 < |x - x_0| < \delta$, we have
    \[ |f(x)| \le M|g(x)| .\]
    So $\frac{f(x)}{g(x)}$ is bounded for $x$
    sufficiently close to $x_0$.
\end{definition*}

\begin{remark*}[Using $=$ for $O$ and $O$]
    This is an abuse of notation, because $O(g(x))$
    or $o(g(x))$ does note denote a function: it
    denotes a class of functions. This can be
    though of more as $f(x) = o(g(x))$ means
    $f(x) \in o(g(x))$.
\end{remark*}

The definition of $O$ can be extended to behaviour
at $\infty$: we say that $f(x)$ is $O(g(x))$ as $x
\to \infty$ if $\exists X > 0$ and $M > 0$ such
that for all $x > X$ we have
\[ |f(x)| \le M|g(x)| .\]

\begin{remark*}[$O \subset o$]
    $f(x) = o(g(x)) \implies f(x) = O(g(x))$, but
    not the other way around.
\end{remark*}

\begin{example*}
    In this case we have $f(x) = O(x)$ but $f(x)
    \neq o(x()$ as $x \to 0$ since the limit is 2
    which is not 0.
\end{example*}

\begin{example*}
    We have that $x^2 + x = O(x^2)$ as $x \to
    \infty$, since for $x > 1$, we have that
    $|x^2 + x| < 2|x^2|$. In general, when we
    focus on behaviour of a polynomial at
    infinity, we see that the highest power term
    dominates.
\end{example*}

\begin{example*}
    $\sin 2x = O(x)$ as $x \to 0$ since
    $\sin 2x \approx 2x$ for small $x$.
\end{example*}

Order parameters are useful to classify remainder
terms before taking limits.

\begin{example*}
    \[ \left. \dfrac{f}{x} \right|_{x_0}
    = \frac{f(x_0 + h) - f(x_0)}{h} +
    \frac{o(h)}{h} \qquad \text{as } h \to 0 \]
    (this equation can be understood by
    considering the limit as $h \to 0$.) It
    follows from this equation that
    \[ f(x_0 + h) = f(x_0) + h \left. \dfrac{f}{x}
    \right|_{x_0} +o(h) .\]
\end{example*}

\subsection{Rules for differentiation}

\begin{theorem*}[Chain rule]
    Given $f(x) = F(g(x))$ we have that
    \[ \dfrac{f}{x} = F'(g(x)) \times \dfrac{g}{x} =
    \dfrac{F}{g} \times \dfrac{g}{x} \]
\end{theorem*}

\begin{proof}
    See printed notes.
\end{proof}

\begin{example*}
    \[ \dfrac{}{x} \sin(x^2 + x - 2)
    = \underbrace{\cos(x^2 + x - 2)}_{F'(g(x))}
    \underbrace{(2x + 1)}_{g'(x)} \]
\end{example*}

\begin{theorem*}[Product rule]
    Given $f(x) = u(x)v(x)$ we have
    \[ \dfrac{f}{x} = v \dfrac{u}{x}
    + u \dfrac{v}{x} \]
\end{theorem*}

\begin{proof}
    Left as an exercise.
\end{proof}

\begin{example*}[Quotient rule]
    The Quotient Rule is a special case of product
    rule:
    \[ \dfrac{}{x} \left( \frac{u}{v} \right)
    = \frac{1}{v} u' + u \dfrac{}{x} \left(
    \frac{1}{v} \right) = \frac{u'}{v} -
    \frac{v'u}{v^2} = \frac{u'v + v'u}{v^2} \]
\end{example*}

\begin{remark*}
    For any real number $\lambda$, if a function
    $f(x) = o(g(x))$, then we also have that
    $\lambda f(x) = o(g(x))$.
\end{remark*}

\subsubsection*{Leibniz's rule}

Consider $f(x) = u(x) v(x)$; then
\begin{align*}
    f' &= u'x + uv' \\
    f'' &= u''v + \underbrace{u'v' + u'v'}_{2u'v'}
    + uv'' \\
    f''' &= u'''v + \underbrace{u''v' +
    2u''v'}_{3u''v'} + \underbrace{2u'v'' +
    u'v''}_{3u'v''} + uv'''
\end{align*}

\begin{theorem*}[Leibniz's rule]
    Given $f(x) = u(x)v(x)$ then
    \begin{align*}
        f^{(n)}(x) &= \sum_{r = 0}^n
        {n \choose r} u^{(n - r)}(x)v^{(r)}(x) \\
        &= u^{(n)}v + nu^{(n - 1)}v^{(1)} +
        \frac{n(n - 1)}{2!} u^{(n - 2)} v^{(2)} \\
        &\,\,\,\,\,\, + \cdots + \frac{n!}{m!(n-m)!}
        u^{(n - m)}v^{(m)} + \cdots + uv^{(n)}
    \end{align*}
\end{theorem*}

\begin{proof}
    Induction.
\end{proof}

\subsection{Taylor Series}

\subsubsection*{Taylor's Theorem}

Previously:
\[ f(x_0 + h) = f(x_0) + hf'(x_0) + o(h)
\qquad \text{as $h \to 0$} \]
This extends to the following theorem:
\begin{theorem*}[Taylor's Theorem]
    For $n$-times differentiable $f(x)$, then
    \[ f(x_0 + h) = f(x_0) + hf'(x_0) +
    \frac{h^2}{2!} f''(x_0) + \cdots
    + \frac{h^{(n)}}{n!}f^{(n)}(x_0) + E_n \]
    where $E_n$ is $o(h^n)$ as $h \to 0$.
\end{theorem*}

\noindent Stronger version if $f^{(n + 1)}(x)$ exists
$\forall x \in (x_0, x_0 + h)$ and $f^{(n)}$ is
continuous in this range:
\begin{align*}
    E_n &= O(h^{n + 1}) &&\text{as $H \to 0$} \\
    &= f^{(n + 1)}(x_n) \frac{h^{n + 1}}{(n + 1)!}
    &&\text{for some $x_0 \le x_n \le x_0 + h$}
\end{align*}

\begin{remark*}
    $E_n = O(h^{n + 1})$ is a stronger statement
    than $o(h^n)$. For example
    \[ h^{n + 1/2} = o(h(n)) \qquad \text{but}
    \qquad h^{n + 1/2} \neq O(h^{n + 1}) \quad
    \text{as $h \to 0$} \]
\end{remark*}

\subsubsection*{Taylor Polynomials}

With $x_0 + h = x$, Taylor's theorem gives
\[ f(x) = f(x_0) + hf'(x_0) + \cdots
+ \frac{(x - x_0)^n}{n!}f^{(n)}(x_0) + E_n .\]
This is called the $n$-th order Taylor polynomial
of $f(n)$ about $x_0$.
\begin{itemize}
    \item $n$-th order Taylor polynomial matches
        first $f$ derivatives of $f(x)$ at $x_0$.
    \item Provides local approximation to $f(x)$
        in vicinity of $x_0$, with error
        $E_n = O(h^{n + 1}$.
    \item If Taylor polynomial converges (as $x$
        varies) as $n \to \infty$, gives
        \emph{Taylor series} of $f(x)$ about
        $x_0$.
\end{itemize}

\subsection{L'H\^opital's Rule}

L'H\^opital's Rule is used to deal with
indeterminate forms in limits.

\begin{theorem*}[L'H\^opital's Rule]
    Let $f(n)$ and $g(n)$ be differentiable at
    $x_0$, with continuous first derivatives
    there, and
    \[ \lim_{x \to x_0} = f(x_0) = 0 \qquad
    \text{and} \qquad \lim_{x \to x_0} g(x) =
    g(x_0) = 0 \]
    Then if $g'(x_0) \neq 0$,
    \[ \lim{x \to x_0} = \frac{f(x)}{g(x)} =
    \lim_{x \to x_0} \frac{f'(x)}{g'(x)} \]
    \emph{provided} the limit on the right exists.
\end{theorem*}

\begin{proof}
    Note that
    \[ f(x) = \underbrace{f(x_0)}_{0} + (x - x_0)f'(x_0)
    + o(x - x_0) \qquad \text{as $x \to x_0$} \]
    and similarly
    \[ g(x) = \underbrace{g(x_0)}_{0} + (x - x_0)g'(x_0)
    + o(x - x_0) \qquad \text{as $x \to x_0$} .\]
    Then
    \begin{align*}
        \lim{x \to x_0} \frac{f(x)}{g(x)} &=
        \frac{f'(x_0) + \frac{o(x - x_0)}{x - x_0}}
        {\underbrace{g'(x_0)}_{\neq 0}
        + \frac{o(x - x_0)}{x - x_0}} \\
        &= \frac{f'(x_0)}{f'(x_0)} \\
        &= \lim_{x \to x_0} \frac{f'(x)}{g'(x)} 
    \end{align*}
\end{proof}

\begin{remark*}[Generalisation]
    The proof above also applies more generally,
    but is easier to understand for the special
    case mentioned in L'H\^opital's Rule. In
    particular
    \[ \lim_{x \to x_0} \frac{f(x)}{g(x)}
    = \lim_{x \to x_0} \frac{f''(x)}{g''(x)} \]
    provided the limit on the right side exists.
\end{remark*}

\begin{example*}
    Consider
    \begin{align*}
        f(x) &= 2\sin x - \sin 3x \\
        g(x) &= 2x - \sin 2x
    \end{align*}
    Notice that we have $f(0) = g(0) = 0$ and we
    have
    \begin{align*}
        f'(x) &= 3\cos x - 3\cos 3x \qquad &&
        (= 0 \text{ at } x = 0) \\
        g'(x) &= 2 - 2\cos 2x &&(= 0 \text{ at } x = 0)
    \end{align*}
    Note that
    \[ f''(0) = g''(0) = 0
    \qquad \text{ at } x = 0 \]
    However one can compute that
    \[ f'''(0) = 24 \qquad \text{and} \qquad
    g'''(0) = 8 \]
    From this it follows that
    \[ \lim_{x \to 0} \frac{f(x)}{g(x)}
    = \lim{x \to 0} \frac{f'''(x)}{g'''(x)}
    = \frac{24}{8} = 3 .\]
\end{example*}

The following is a non-examinable sketch-proof of
Taylor's Theorem:
\begin{proof}
    Start from the Fundamental Theorem of Calculus
    (proved in the next section)
    \begin{align*}
        \int_)^x f'(t) \dd t &= f(x) - f(0) \\
        \implies f(x) &= f(0) = \int_0^x f'(t)
        \dd t \\
        &= f(0) + \int_0^x \dfrac{}{t} (t -
        x)f'(t) \dd t \\
        &= f(0) + \left[ (t - x)f'(t) \right]_0^x
        - \int_0^x (t - x)f''(t) \dd t \\
        &= f(0) + xf'(0) - \half \int_)^x
        \dfrac{}{t} (T - x)^2 f''(t) \dd t \\
        &= f(0) + xf'(0)
        - \half \left[ (t - x)^2f''(t) \right]
        _0^x + \half \int_0^x (t - x)^2 f'''(t)
        \dd t \\
        &\hspace{15em}\vdots \\
        &= f(0) = xf'(0) + \frac{x^2}{2!}f''(0) +
        \cdots + \frac{x^n}{n!}f^{(n)}(0)
        + \frac{1}{n!} \int_0^x (x - t)^n f^{(n +
        1)} (t) \dd t
    \end{align*}
    And we can note that the integral on the last
    line is $E_n = o(h^n)$.
\end{proof}

\section{Integration}

\subsection{Integrals as Riemann Sums}

The aim of this section is to formalise ``area
under a curve''.

\begin{definition*}[Integral]
    The \emph{integral} of a suitably well behaved
    function $f(x)$ is the limit of a sum:
    \[ \in_a^b f(x) \dd x \equiv
    \lim_{\Delta x \to 0} \sum_{n = 0}^{N - 1}
    f(\underbrace{x_n}_{a + n\Delta x}) \underbrace{
    \Delta x}_{\frac{b - a}{N}} \]
    \begin{center}
    \begin{tsqx}
        ! size(5cm);
        (0,0)->(1.5,0)
        (0,0)->(0,1.5)
        label $x$ @ (1.6,0)
        label $a$ @ (0.2,-0.1)
        label $b$ @ (1.3,-0.1)
        label $x_1, x_2, \dots, x_{N - 1}$ @ (0.75,-0.1)
        ! real f(real x) {return 0.5*x*x + 0.4;}
        ! draw(graph(f,-0.1,1.4,operator ..), red);
        (0.2,-0.02)--(0.2,-0.05)
        (0.7,-0.02)--(0.7,f(0.6))
        ! for(int i; i < 5; ++i) {
            (0.2+0.1*i,f(0.2+0.1*i))--(0.2+0.1*i,-0.02)
            (0.2+0.1*i,f(0.2+0.1*i))--(0.3+0.1*i,f(0.2+0.1*i))
        ! }
        label $\cdots$ @ (0.95,0.4)
        label $f(x)$ @ (0.95, f(0.95)+0.2)
        (1.2,-0.02)--(1.2,f(1.2))--(1.3,f(1.2))--(1.3,-0.05)
    \end{tsqx}
    \end{center}
    \textit{Riemann sum: Limit should not depend
    on exact choice of rectangles (e.g. $\Delta x$
    does not have to be uniform). Relation to
    ``area under curve''?}\\
\end{definition*}
\noindent
Consider one rectangle for finite $N$:
\begin{center}
\begin{tsqx}
    ! size(5cm);
    (0,0)->(0,1.5)
    (0,0)->(1.5,0)
    ! real f(real x) {return 0.3*x*x + 0.7;}
    ! draw(graph(f,0,1.4,operator ..), red);
    label $x$ @ (1.6,0)
    label $x_n$ @ (0.4,-0.1)
    label $x_{n+1}$ @ (1.1,-0.1)
    (0.4,-0.05)--(0.4,f(0.4))--(1.1,f(0.4))--(1.1,-0.05)
    label $f(x)$ @ (1.2,f(1.2)+0.15)
    (0.4,f(0.4))--(0.4,f(1.1))--(1.1,f(1.1))--(1.1,f(0.4)) dashed
\end{tsqx}
\end{center}
One can see that the solid rectangle has area less
than the area under the curve, and the dashed one
has greater area. One might ask whether we could
pick a point to construct the rectangle from such
that the area of the rectangle and the area of the
curve are the same. The mean value theorem states
that this is possible.

\begin{theorem*}[Mean-value theorem]
    For $f(x)$ continuous, then if $A_n$ is the
    area under the curve between $x_n$ and $x_{n +
    1}$, then we have
    \[ A_n = (x_{n + 1} - x_n)f(c_n) \]
    for some $c_n$ satisfying $x_n \le c_n \le
    x_{n + 1}$. (See Analysis I for proof).
    If $f(x)$ is differentiable then
    \begin{align*}
        f(c_n) &= f(x_n) + O(c_n - x_n) \qquad
        \text{as $c_n - x_n \to 0$} \\
        &= f(x_n) + O(\Delta x) \qquad \text{since
        $c_n \le x_n + \Delta x$}
    \end{align*}
    It follows that
    \[ A_n = \Delta x f(x_n) + O(\Delta x^2)
    \qquad \text{as $\Delta x \to 0$} \]
    Total area between $a$ and $b$:
    \begin{align*}
        A &= \lim_{N \to \infty} \sum_{n = 0}^N A_n
        \\
        &= \underbrace{\lim_{N \to \infty} \sum_{n
        = 0}^N f(x_n) \Delta x}_{\text{Definition of
        integral}}
        + \underbrace{\lim_{N \to \infty} NO\left( \left(
        \frac{b - a}{N} \right)^2 \right)}
        _{\lim_{N \to \infty} O\left( \frac{(b -
        a)^2}{N} \right)} \\
        &= \int_a^b f(x) \dd x + 0
    \end{align*}
\end{theorem*}

\subsection{Fundamental Theorem of Calculus}

Formalise ``inverse of differentiation''

\begin{theorem*}[FTC]
    Let $F(x)$ be defined as
    \[ F(x) = \int_a^x f(Tt) \dd t \]
    Then
    \[ \dfrac{F}{x} = f(x) .\]
\end{theorem*}

\begin{proof}
    \begin{align*}
        \dfrac{F}{x} &= \lim_{h \to 0}
        \frac{F(x + h) - F(x)}{h} \\
        &= \lim_{h \to 0} \frac{1}{h} \left[ 
        \int_a^{x + h} f(t) \dd t - \int_a^{x}
        f(t) \dd t \right] \\
        &= \lim_{h \to 0} \frac{1}{h} \int_x^{x +
        h} f(t) \dd t \\
        &= \lim_{h \to 0} \frac{1}{h} \left[ f(x)h
        + O(h^2) \right] \qquad \text{(using Mean
        Value Theorem)} \\
        &= f(x) + \lim_{h \to 0} O(h) \\
        &= f(x)
    \end{align*}
\end{proof}

\begin{remark*}
    $F(x)$ is the solution of $\dfrac{F}{x} =
    f(x)$ with $F(a) = 0$.
\end{remark*}
\noindent
Corollaries:
\[ \dfrac{}{x} \int_a^b f(t) \dd t = -f(x) \]
\[ \dfrac{}{x} \int_a^{g(x)} f(t) \dd t = f(g(x))
\dfrac{g}{x} \]

\begin{notation*}[Indefinite Integrals]
    We may write indefinite integrals as
    \[ \int f(x) \dd x \qquad \text{or} \qquad
    \int^x f(t) \dd t \]
    (note that in the second case the lack of
    lower limit gives rise to integration
    constant.)
\end{notation*}

\subsection{Methods of integration}

\subsubsection*{Integration by subs}

Useful when integrand is ``function of a
function''

\begin{example*}
    Consider 
    \[ I = \int \frac{1 - 2x}{\sqrt{x - x^2}} \dd
    x \]
    Let $u = x - x^2 \implies \dfrac{u}{x} = 1 -
    2x$
    \[ I = \int \frac{1}{\sqrt{u}} = 2\sqrt{u} + c
    = 2\sqrt{x - x^2} + c \]
\end{example*}
\noindent
Here are some useful trig / hyperbolic subs:
\[ \sin^2 \theta + \cos^2\theta = 1 \qquad 
\cosh^2 u \sinh^2 u = 1 \]
\[ \tan^2\theta + 1 = \sec^2\theta \qquad 
1 - \tanh^2 u = \sech^2 u\]
Generally we use the subs $\sin\theta$ or $\tanh
u$, $\tan\theta$ or $\sinh u$, $\sec\theta$ or
$\cosh u$ for terms $1 - x^2$, $1 + x^2$, $x^2 -
1$ respectively.

\begin{example*}
    \[ I = \int \sqrt{2x - x^2} \dd x \]
    Try $x - 1 = \sin\theta \implies \dd x = \cos
    \theta \dd\theta$. Then
    \begin{align*}
        I &= \int \sqrt{1 - \sin^2\theta}
        \cos\theta \dd\theta \\
        &= \int \cos^2\theta \dd\theta \\
        &= \half \int 1 + \cos 2\theta \dd\theta
        \\
        &= \half(\theta + \half \sin 2\theta) + c
        \\
        &= \half(\theta + \sin\theta\cos\theta) +
        c \\
        &= \half \arcsin(x - 1) + \half(x -
        1)\sqrt{1 - (x - 1)^2} + 1
    \end{align*}
\end{example*}

\subsection{Integration by parts}

Recall the product rule: $uv' = (uv)' - u'v$.

\begin{theorem*}[Integration by parts]
    For functions $u$ and $v$
    \[ \int uv' \dd x = uv - \int u'v \dd x \]
\end{theorem*}

\begin{example*}
    \begin{align*}
        I &= \int_0^\infty xe^{-x} \dd x \\
          &= [ -xe^{-x} ]_0^\infty + \int_0^\infty
          e^{-x} \dd x \\
          &= 1
    \end{align*}
\end{example*}

\begin{example*}
    \begin{align*}
        I &= \int \ln x \dd x \\
        &= x \ln x - \int \frac{1}{x} x \dd x \\
        &= x \ln x - x + c
    \end{align*}
\end{example*}

\section{Partial Differentiation}

\subsection{Functions of Several Variables}

\begin{definition*}[Multivariate functions]
    A function is called \emph{multivariate} if it
    depends on more than one variable.
\end{definition*}
\noindent
Some examples of multivariate functions include
\begin{itemize}
    \item Height of terrain: $h(\text{latitude},
        \text{longitude})$.
    \item Density of air: $\rho(x, y, z, t)$.
\end{itemize}

\begin{center}
\begin{tsqx}
    (-0.2,0)--(1.5,0) EndArrow
    (0,-0.2)--(0,1.5) EndArrow
    label $y$ @ (1.5,-0.15)
    label $x$ @ (0.15,1.5)
    CR (0.7,0.7) 0.2
    CR (0.7,0.7) 0.4
    CR (0.7,0.7) 0.6
    label $A^*$ @ (0.7,0.4)
    label $f(x, y)$ @ (1.3,1.3)
\end{tsqx}
\end{center}
The slope of $f$ at $A$ depends on direction.

\subsection{Partial Derivatives}

\begin{definition*}[Partial derivative]
    Given some function of several variables, for
    example $f(x, y)$, the \emph{partial
    derivative} of $f$ with respect to $x$ at
    fixed $y$ is
    \[ \left. \pfrac{f}{x} \right|_y
    = \lim_{\delta x \to 0} \left[
    \frac{f(x + \delta x, y) - f(x, y)}{\delta x}
    \right] \]
\end{definition*}
\noindent
``Slope of $f$ when moving in $x$ direction.''

Note that $\left. \pfrac{f}{x} \right|_x$ is defined
similarly.

\begin{example*}
    Let $f$ be defined as
    \[ f(x, y) = x^2 + y^2 + e^{xy^2} \]
    Then we have that
    \[ \left. \pfrac{f}{x} \right|_y
    = 2x + 0 + y^2e^{xy^2} \]
    \[ \left. \pfrac{f}{x} \right|_x
    = 0 + 2y + 2xye^{xy^2} \]
    We can also compute second partial
    derivatives:
    \[ \left. \pfrac[2]{f}{x} \right|_y
    = 2 + y^2 \cdot y^2 \cdot e^{xy^2} \]
    \[ \left. \pfrac[2]{f}{x} \right|_x
    = 2 + 2xe^{xy^2} + 2xy \cdot 2xy e^{xy^2} \]
    and we can compute \emph{mixed} derivatives
    \[ \left. \pfrac{}{x} \left( \left.
    \pfrac{f}{y} \right|_x \right) \right|_y
    = 2ye^{xy^2} + 2xy \cdot y^2 \cdot e^{xy^2} \]
    \[ \left. \pfrac{}{y} \left( \left.
    \pfrac{f}{x} \right|_y \right) \right|_x
    = 2ye^{xy^2} + y^2 \cdot 2xy \cdot e^{xy^2} \]
\end{example*}

\begin{notation*}
    A useful convention that is often used is to
    omit the $|_y$ symbol if all other
    variables are being held fixed. For example
    \[ \frac{\partial^2 f}{\partial x \partial y}
    = \left. \pfrac{}{x}
    \left( \left. \pfrac{f}{y} \right|_x
    \right) \right|_y \]
\end{notation*}

Note that if a function as continuous (mixed)
second derivatives then
\[ \frac{\partial^2 f}{\partial x\partial y}
= \frac{\partial^2 f}{\partial y \partial x} \]
This is known as ``Schwarz's theorem''. (In other
words, partial derivatives commute). \\\\
Alternative notation:
\[ f_x = \pfrac{f}{x} \qquad;\qquad
f_{x, y} = \frac{\partial^2 f}{\partial y \partial
x} \]

\subsection{Multivariate Chain Rule}

Given path $x(t)$, $y(t)$ and $f(x, y)$, what is
$\dfrac{f}{t}$? \\\\
Consider change in $f$ under $(x, y) \to (x +
\delta x, y + \delta y)$:
\begin{align*}
    \delta f
    &= f(x + \delta x, y + \delta x) - f(x, y) \\
    &= [f(x + \delta x, y + \delta y) - f(x +
    \delta x, y)]
    + [f(x + \delta x, y) - f(x, y)]
\end{align*}
Then Taylor expanding we get
\[ f(x + \delta x, y) - f(x, y)
= f_x(f, y) \delta x + o(\delta x) \]
and
\[ f(x + \delta x, y + \delta y) - f(x + \delta x,
y) = f_y(x + \delta x, y) \delta y + o(\delta y) \]
We also Taylor expand $f_y(x + \delta x, y)$ and
find that
\[ f_y(x + \delta x, y) = f_y(x, y) + f_{y, x}(x,
y) \delta x + o(\delta x) \]
And hence
\[ \delta f = [f_y(x, y) + f_{y, x}(x, y)\delta x
+ o(\delta x)]\delta y + o(\delta y) + f_x(x, y)
\delta x + o(\delta x) \hspace{3em} (*) \]
Take limit as $\delta x, \delta y \to 0$ to get
\emph{differential} of $f$:
\[ \dd f = \lim_{\substack{\delta x \to 0 \\\delta
y \to 0}} \delta f \]

\begin{theorem*}[Chain Rule for Partial Derivatives]
    Differential $\dd f$ of $f(x, y)$ is related
    to differentials of arguments $\dd x$ and $\dd
    y$, as
    \[ \dd f = \pfrac{f}{x} \dd x + \pfrac{f}{x}
    \dd y \]
    This is because the remaining terms in $(*)$
    approach 0 faster than $\dd x$ or $\dd y$.
\end{theorem*}
\noindent
For the path $x(t), y(t)$ we now have
\[ \dfrac{}{t}f(x(t), y(t))
= \pfrac{f}{x} \dfrac{x}{t} + \pfrac{f}{y}
\dfrac{y}{t} \]
If path parametrised as $y(x)$, we have $f(x,
y(x))$ and
\[ \dfrac{f}{x} = \pfrac{f}{x}
\ub{\dfrac{x}{f}}_{=1} + \pfrac{f}{y} \dfrac{y}{x}
= \pfrac{f}{x} + \pfrac{f}{y} + \dfrac{y}{x} .\]

\subsubsection*{Integral Form of the Chain Rule}

\[ \Delta f = \int \dd f = \int \left(
\pfrac{f}{x} \dd x = \pfrac{f}{y} \dd y \right) \]
Note that $\Delta f$ denotes the change in $f$
between endpoints.
For a path $x(t), y(t)$ connecting the endpoints,
have $f(x(t), y(t))$
\[ \Delta f = \int \left( \pfrac{f}{x} \dfrac{x}{t}
+ \pfrac{f}{y} \dfrac{y}{t} \right) \dd t \]
Note that the result doesn't depend on path chosen
between two given endpoints.

\subsection{Applications of Multivariate Chain Rule}

\subsubsection*{Change of Variables}

\begin{example*}
    Cartesian coordinates $(x, y)$ to plane-polar
    $(r, \theta)$.
    \begin{center}
    \begin{tsqx}
        (-0.2,0)--(1.5,0) EndArrow
        (0,-0.2)--(0,1.5) EndArrow
        label $x$ @ (1.5,-0.15)
        label $y$ @ (0.15,1.5)
        P .= (1.2,0.9)
        label $(x, y)$ @ P+(0.15,0)
        label $r$ @ 0.5*P+(-0.1,0.1)
        anglemark4 (1,0) origin P 3
        label $\theta$ @ 0.2*dir(20)
        origin--P red
        label $x = r\cos\theta$ @ (2,1)
        label $y = r\sin\theta$ @ (2,0.5)
    \end{tsqx}
    \end{center}
    Original function $f(x, y)$ can be thought of
    as a function of $r$, $\theta$:
    \[ f(x, y) = f(x(r, \theta), y(r, \theta)) \]
    \[ \left. \pfrac{f}{r} \right|_y
    = \left. \pfrac{f}{x} \right|_y \left.
    \pfrac{x}{r} \right|_\theta
    + \left. \frac{f}{y} \right|_x \left.
    \pfrac{y}{r} \right|_y \]
    \[ \left. \pfrac{f}{\theta} \right|_r
    = \left. \pfrac{f}{x} \right|_y \left.
    \pfrac{x}{\theta} \right|_r + \left.
    \pfrac{f}{y} \right|_x \left.
    \pfrac{y}{\theta} \right|_r \]
\end{example*}

\subsubsection*{Implicit Differentiation}

Consider $f(x, y, z) = z$ which is a surface in
3D space. Implicitly define
\[ z = z(x, y), \qquad x = x(y, z), \qquad y =
y(x, z) \]
(Note that we may not be able to find these
explicitly but can still determine partial
derivatives).

\begin{example*}
    Consider
    \[ xy + y^2z + z^5 = 1 \tag{$*$} \]
    Take derivative of $(*)$ with respect to $x$
    holding $y$ fixed
    \[ y + y^2 \left. \pfrac{z}{x} \right|_y
    + 5z^4 \left. \pfrac{z}{x} \right|_y = 0 \]
    \[ \implies \left. \pfrac{z}{x} \right|_y
    = \frac{-y}{y^2 + 5z^4} .\]
\end{example*}
\noindent
Generally given $f(x, y, z) = c$, MVC gives
\[ 0 = \dd f = \left. \pfrac{f}{x} \right|_{y, z}
\dd x
+ \left. \pfrac{f}{y} \right|_{x, z} \dd y
+ \left. \pfrac{f}{z} \right|_{x, y} \dd z \]
Note that we can't vary $x$, $y$, $z$
independently \emph{and} stay in the surface.

Find $\left. \pfrac{z}{x} \right|_y$ by
\[ \left. \pfrac{f}{x} \right|_{y, z}
\left. \pfrac{x}{x} \right|_y
+ \left. \pfrac{f}{y} \right|_{x, z}
\left. \pfrac{y}{x} \right|_y
+ \left. \pfrac{f}{z} \right|_{x, y}
\left. \pfrac{z}{x} \right|_y \]
\[ \left. \pfrac{z}{x} \right|_y
= \frac{- \partial f / \partial x \big|_{y, z}}
{\partial f / \partial x \big|_{x, y}}
\tag{$\cross$} \]
Similarly,
\[ \left. \pfrac{x}{y} \right|_z
= \frac{-\partial f / \partial y \big|_{x, z}}
{\partial f / \partial y \big|_{y, z}}\]
and
\[ \left. \pfrac{y}{z} \right|_x
= \frac{-\partial f / \partial z \big|_{x, z}}
{\partial f / \partial y \big|_{x, z}} \]
Follows that
\[ \left. \pfrac{x}{y} \right|_z
\left. \pfrac{y}{z} \right|_x
\left. \pfrac{z}{x} \right|_y = -1 \]

\begin{theorem*}[Reciprocal Rule]
    Applies to partial derivatives provided some
    variables held fixed: \\
    Similarly to $(\cross)$
    \[ \left. \pfrac{x}{y} \right|_z
    = \frac{-\partial f / \partial z \big|_{x, z}}
    {\partial f / \partial x \big|_{y, z}} \]
    so that
    \[ \left. \pfrac{z}{x} \right|_y
    = \frac{1}{\partial x / \partial z \big|_y} .\]
\end{theorem*}

\subsubsection*{Differentiation of an integral with respect to a parameter}

Consider a family of functions $f(x; \alpha)$.
Define the integral
\[ I(\alpha) = \int_{a(\alpha)}^{b(\alpha)} f(x;
\alpha) \dd x \]
What is $\dd I / \dd \alpha$?

\begin{theorem*}[Differentiation of integral with respect to a parameter]
    \[ \dfrac{I}{\alpha} =
    \int_{a(\alpha)}^{b(\alpha)} \pfrac{f}{\alpha}
    (x; \alpha) \dd x
    + f(b(\alpha); \alpha) \dfrac{b}{\alpha}
    -f(a(\alpha); \alpha) \dfrac{a}{\alpha} \]
\end{theorem*}

\begin{proof}
    \begin{align*}
        \dfrac{I}{\alpha}
        &=
        \lim_{\delta \alpha \to 0}
        \frac{1}{\delta
        \alpha} \left[ \int_{a(\alpha + \delta
        a)}^{b(\alpha + \delta \alpha)}
        f(\alpha; \alpha + \delta \alpha) \dd x
        - \int_{a(\alpha)}^{b(\alpha)} f(x;
        \alpha) dd x \right]  \\
        &= \lim_{\delta x \to 0}
        \frac{1}{\delta \alpha}
        \left[ \int_{a(\alpha)}^{b(\alpha)}
        f(x; \alpha + \delta \alpha) - f(x; \alpha)
        \dd x \right] \\
        &\,\,\,\,\,\, + \lim_{\delta \alpha \to 0}
        \frac{1}{\delta \alpha} \int_{b(\alpha)}
        ^{b(\alpha + \delta\alpha)} f(\alpha;
        \alpha + \delta \alpha) \dd x
        - \lim_{\delta \alpha \to 0}
        \frac{1}{\delta \alpha}
        \int_{a(\alpha)}^{b(\alpha +
        \delta\alpha)} f(x; \alpha + \delta\alpha)
        \dd x\\
        &\,\,\,\,\,\, + \lim_{\delta \alpha \to 0}
        \frac{1}{\delta\alpha}
        \int_{a(\alpha)}^{a(\alpha +
        \delta\alpha)} f(\alpha; \alpha + \delta
        \alpha) \dd x \\
        &= \int_{a(\alpha)}^{b(\alpha)}
        \pfrac{f}{\alpha} (x; \alpha) \dd x
        + f(b(\alpha); \alpha) \dfrac{b}{\alpha}
        - f(a(\alpha); \alpha) \frac{a}{\alpha}
    \end{align*}
\end{proof}

\begin{example*}
    \[ I(\lambda) = \int_0^\lambda e^{-\lambda
    x^2} \dd x \]
    \[ \dfrac{I}{x} = \int_0^\lambda -x^2
    e^{-\lambda x^2} \dd x +
    \dfrac{\lambda}{\lambda} e^{-\lambda^2} \]
\end{example*}

\begin{example*}
    Suppose we want to evaluate
    \[ \int_0^\infty x^n e^{-x} \dd x \qquad
    \text{($n$ an integer)} \]
    Let
    \[ I(\lambda) = \int_0^\infty e^{-\lambda x}
    \dd x = \frac{1}{\lambda} \]
    \[ \dfrac[n]{}{\lambda} I = \int_0^\infty
    (-x)^n e^{-\lambda n} \dd x = (-1)^n
    \frac{n!}{\lambda^{n + 1}} \]
    Now set $\lambda = 1$ to get
    \[ \int_0^\infty x^n e^{-x} \dd x = n! \]
\end{example*}

\newpage
\begin{center}
    \LARGE \bfseries
    TOPIC II
\end{center}

\newpage
\setcounter{section}{0}
\section{First Order Linear Differential Equations}

\begin{definition*}[Exponential function]
    Exponential function is defined by the
    infinite series
    \[ \exp(x) = 1 + x + \frac{x^2}{2!}
    + \frac{x^3}{3!} + \cdots
    = \sum_{r = 0}^\infty \frac{x^r}{n!} \]
    It can also be written as
    \[ \exp(x) = \lim_{k \to \infty}
    \left( 1 + \frac{x}{k} \right)^k
    = \lim_{k \to \infty} \left( 1
    + k\frac{x}{k} + \frac{k(k - 1)}{2!}
    \left( \frac{x}{k} \right)^2 + \cdots \right) \]
\end{definition*}
\noindent
Differentiating the series gives
\[ \dfrac{}{x} \exp(x) =
1 + 2 \frac{x}{2!} + 3 \frac{x^2}{3!} + \cdots
= \exp(x) \]
As $\exp(0) = 1$, can think of $\exp(n)$ as the
solution of ODE
\[ \dfrac{f}{x} = f \qquad \text{and} \qquad f(0)
= 1 \qquad \implies \int_1^{\exp(x)} \dfrac{t}{t}
= x .\]
Key Property: $\exp(x_1) \exp(x_2) = \exp(x_1 +
x_2)$, which can be verified by expanding.

This allows us to write
\[ \exp(x) = e^x \]
where
\[ e = \exp(1) = \lim_{k \to \infty} \left( 1 +
\frac{1}{k} \right)^k = 2.718\dots \]
We can also write the inverse function, $\ln$ such
that
\[ \exp(\ln x) = x \]
\begin{notation*}
    Sometimes $\ln$ is written as $\log_e$.
\end{notation*}
It follows that
\[ a^x = (a^{\ln a})^x = e^{x \ln a} \]
\[ \dfrac{}{x} a^x = \ln a e^{x \ln a} = (\ln a)
a^x \]

\begin{definition*}[Eigenfunction]
    An \emph{eigenfunction} of the derivative
    operator is a function that is unchanged up to
    a multiplicative scaling by the
    \emph{eigenvalue} under the action of the
    operator
    \[ \dfrac{f}{x} =
    \ub{\lambda}_{\text{eigenvalue}}
    \ub{f(x)}_{\text{eigenfunction}} \]
    Eigenfunctions of $\dfrac{}{x}$ are $\alpha
    e^{\lambda x}$.
\end{definition*}

\newpage
\section{First-order Linear ODEs}

\begin{itemize}
    \item \emph{LINEAR}: dependent variable ($y$,
        say) and its derivatives only appear
        linearly.
    \item \emph{FIRST ORDER}: highest derivative
        of $y$ that appears in $\dfrac{y}{x}$.
\end{itemize}

\subsection{Homogeneous First-order Linear ODEs}

\begin{itemize}
    \item \emph{HOMOGENEOUS}: an ODE in which all
        terms involve dependent variable or its
        derivatives, which implies that $y = 0$ is
        always a solution.
    \item \emph{CONSTANT COEFFICIENTS}:
        independent variable $(x)$ does not appear
        explicitly.
\end{itemize}

\begin{example*}
    \[ 5 \dfrac{y}{x} - 3y = 0 \]
    Try $y = Ae^{\lambda x}$:
    \[ \dfrac{y}{x} = \lambda A e^{\lambda x} =
    \lambda y \]
    So ODE $\implies (5 \lambda - 3) = 0$, so
    $\lambda = \frac{3}{5}$. Since linear,
    homogeneous equation, solution holds for all A.
\end{example*}

\begin{itemize}
    \item Generally, for any linear, homogeneous
        ODE (doesn't have to be 1st order), any
        constant multiple of a solution is a
        solution.
    \item An $n$-th order linear ODE (highest
        derivative is $\dfrac[n]{y}{x}$ has
        exactly $n$ independent solutions, so $y =
        A e^{3x/5}$ is \emph{the} general
        solution.
    \item To specify unique solution requires
        giving suitable boundary conditions; e.g.
        $y(0)$ determines $A$.
\end{itemize}

\subsubsection*{Discrete Equations}

Consider
\[ 5 \dfrac{y}{x} - 3y = 0 \]
with $y(0) = y_0$. Then
\[ y(x) = y_0 e^{3x/5} \]
Now approximate the equation by discrete form at
$\{x_n\}$ where $x_n = nh$.
\begin{center}
\begin{tsqx}
    ! size(5cm);
    (-0.1,0)->(1.5,0)
    (0,-0.05)->(0,1.5)
    label $x$ @ (1.5,-0.15)
    label $y$ @ (0.15,1.5)
    x_0 = (0,-0.1) 0R
    x_1 = (0.3,-0.1) 0R
    x_2 = (0.6,-0.1) 0R
    (0.3,0.05)--(0.3,-0.05)
    (0.6,0.05)--(0.6,-0.05)
    label $\cdots$ @ (0.9,-0.1)
    y_0 ;= (0,0.3) R180
    y_1 ;= (0.3,0.6) R135
    y_2 ;= (0.6,1.1) R90
    y_1--y_2
\end{tsqx}
\end{center}
To do the approximation, we set
\[ \left. \dfrac{y}{x} \right|_{x_n}
= \frac{y_{n + 1} - y_n}{h} \]
so
\[ 5\left( \frac{y_{n + 1} - y_n}{h} \right)  -
3y_n = 0 \]
\[ \implies y_{n + !} = \left( 1 + \frac{3h}{5}
\right) y_n \]
which is an example of a recurrence relation.
\[ y_n = \left( 1 + \frac{3h}{5} \right) y_{n - 1}
= \left( 1 + \frac{3h}{5} \right)^2 y_{n - 2}
= \left( 1 + \frac{3h}{5} \right)^n y_0 \]
It follows that
\[ y_n = y_0 \left( 1 + \frac{3x_n}{5n} \right)^n \]
If we take $x_n \to x$, i.e. take $n$ steps from
$x_0$ to $x$, as $n \to \infty$ we get
\[ \lim_{n \to \infty} y_n = \lim_{n \to \infty}
y_0 \left( 1 + \frac{3x}{5n} \right)^n
= y_0 \exp\left( \frac{3x}{5} \right) \]
(note that this agrees with the continuous case)

\subsubsection*{Series Solution}

\[ y(x) = \sum_{n = 0}^\infty a_n x^n \]
where we determine $a_n$ by substituting into the
ODE.

\begin{example*}
    \[ 5 \dfrac{y}{x} - 3y = 0 \]
    \begin{align*}
        \dfrac{y}{x}
        &= \sum_{n = 0}^\infty a_n n x^{n - 1} \\
        &= \sum_{n = 1}^\infty a_n n x^{n - 1}
    \end{align*}
    so
    \[ x\dfrac{y}{x} = \sum_{n = 1}^\infty a_n n
    x^n \]
    and
    \begin{align*}
        xy
        &= \sum_{n = 0}^\infty a_n x^{n + 1} \\
        &= \sum_{n = 1}^\infty a_{n - 1} x^n
    \end{align*}
    since $5x \dfrac{y}{x} - 3xy = 0$, we have
    \[ \sum_{n = 1}^\infty (5na_n - 3a_{n - 1})
    x^n = 0 \]
    \[ \implies a_n = \frac{3}{5n} a_{n - 1}
    \qquad (n \ge 1) \]
    It follows that
    \begin{align*}
        a_n
        &= \frac{3}{5n} a_{n - 1} \\
        &= \left( \frac{3}{5} \right)^2
        \frac{1}{n(n - 1)} a_{n - 2} \\
        &= \cdots \\
        &= \left( \frac{3}{5} \right)^n
        \frac{a_0}{n!}
    \end{align*}
    and
    \[ y(x) = a_0 \sum_{n = 0}^\infty \frac{1}{n!}
    \left( \frac{3x}{5} \right)^n =
    a_0 \exp\left( \frac{3x}{5} \right) .\]
\end{example*}

\subsection{Forced (inhomogeneous) ODEs}

Include terms involving only the independent
variable so $y = 0$ no longer a (trivial)
solution.

\subsubsection*{Constant forcing}

\begin{example*}
    \phantom{.}\\[-2\baselineskip]
    \[ 5\dfrac{y}{x} - 3y = 10 \]
\end{example*}
\noindent
General method of solution:
\begin{enumerate}[(1)]
    \item Find \emph{any} solution of the forced
        equation $\to$ \emph{particular integral}
        $y_p(x)$. In the example, we could choose
        $y_p(x) = -\frac{10}{3}$, for example.
    \item Write general solution of ODE as
        \[ y(x) = y_p(x) + y_c(x) \]
        where $y_c$ is a complementary function.
        Then $y(x)$ satisfies the ODE if and only
        if $y_c$ satisfies the homogeneous
        version (this only follows because the
        differential equation is \emph{linear}).
        So in the example, this means we need
        \[ 5 \dfrac{y_c}{x} - 3y_c = 0 \]
        \[ \implies y_c(x) = A e^{3x/5} \]
    \item Combine $y_c$ and $y_p$ to get full
        general solution.
        \[ y(x_ = A e^{3x/5} - \frac{10}{3} \]
        ($A$ could be determined by a boundary
        condition if given). This method is
        general for \emph{linear} ODEs.
\end{enumerate}

\subsubsection*{Eigenfunction Forcing}

Forcing term is an eigenfunction of differential
operator on LHS.

\begin{example*}[radioactive decay]
    \phantom{.}
    \begin{center}
    \begin{tsqx}
        ! size(8cm);
        X := (-1,0)
        X+(-0.3,0.2)--X+(0.3,0.2)--X+(0.3,-0.2)
        X+(0.3,-0.2)--X+(-0.3,-0.2)--X+(-0.3,0.2)
        label $A$ @ X
        label $a(t)$ @ X+(0,-0.3)
        Y := (0,0)
        Y+(-0.3,0.2)--Y+(0.3,0.2)--Y+(0.3,-0.2)
        Y+(0.3,-0.2)--Y+(-0.3,-0.2)--Y+(-0.3,0.2)
        label $B$ @ Y
        label $b(t)$ @ Y+(0,-0.3)
        Z := (1,0)
        Z+(-0.3,0.2)--Z+(0.3,0.2)--Z+(0.3,-0.2)
        Z+(0.3,-0.2)--Z+(-0.3,-0.2)--Z+(-0.3,0.2)
        label $C$ @ Z
        label $c(t)$ @ Z+(0,-0.3)
        (-0.65,0)->>(-0.35,0)
        (0.35,0)->>(0.65,0)
        label $k_a a$ @ (-0.5,0.15)
        label $k_b b$ @ (0.5,0.15)
    \end{tsqx}
    \end{center}
    \[ \dfrac{a}{y} = -k_a a \qquad \text{and}
    \qquad \dfrac{b}{t} = k_a a - k_b b \]
    Have $a(t) = a_0 e^{-k_a t}$
    \[ \dfrac{b}{t} + k_b b = k_a a_0 e^{-k_a t}
    \tag{$*$} \]
    Try particular integral
    \[ b_p(t) = C e^{-k_a t} \]
    \[ (*)\qquad\qquad (k_s - k_a) C e^{-k_a t}
    = k_a a_0 e^{-k_a t} \]
    \[ \implies C = \frac{k_a a_0}{k_s - k_a}
    \qquad (k_a \neq k_s) \]
    General solution $b(t) = b_c(t) + b_p(t)$
    where
    \[ \dfrac{b_c}{t} + k_s b_c = 0 \implies
    b_c(t) = D e^{-k_s t} \]
    Follows that
    \[ b(t) = D e^{-k_s t} + \frac{k_a}{k_s - k_a}
    a_0 e^{-k_a t} \]
    In the case that $b(0) = 0$:
    \[ b(t) = \frac{b_a a_0}{k_s - k_a} (e^{-k_a
    t} - e^{k_s t}) \]
    \begin{center}
    \begin{tsqx}
        ! size(5cm);
        (-0.1,0)->>(1.5,0)
        (0,-0.1)->>(0,1.5)
        label $t$ @ (1.6,0)
        label $a_0$ @ (-0.1,1.2)
        label $a(t)$ @ (1,0.6) +col:blue
        label $b(t)$ @ (0.3,0.1) +col:red
        ! real f(real x) {return 1.5*(0.5**x) - 0.3;}
        ! draw(graph(f,0,1.4,operator ..), blue);
        ! real g(real x) {return x / (1 + 100*x**6);}
        ! draw(graph(g,0,1.4,operator ..), red);
    \end{tsqx}
    \end{center}
    If $k_a = k_s$, need a different particular
    integral.
\end{example*}

\subsection{Non-constant coefficients}

General linear first-order ODE
\[ a(X) \dfrac{y}{x} + b(x) y = c(x) \]
``Standard form''
\[ \dfrac{y}{x}  + p(x)y = f(x) \]
Solve using \emph{integrating factor} $\mu(x)$;
multiply ODE by $\mu(x)$:
\[ \ub{\mu y' + (\mu p) y)}_{\dfrac{}{x} \text{if
$\mu' \mu p$}} = \mu f(x) \tag{$*$} \]
Requires
\[ \ub{\int \frac{1}{\mu} \dfrac{\mu}{x} \dd x}
_{\ln \mu}
= \int p \dd x \]
\[ \implies \mu(x) = \exp\left( \int^x p(u) \dd u \right) \]
($\mu$ unique up to an irrelevant constant factor)

\[ (*) \implies \dfrac{}{x}(\mu y) = \mu f(x) \]
so
\[ \mu y = \int^x \mu f \dd u \]

\begin{example*}
    \[ x \dfrac{y}{x} + (1 - x)y = 1 \]
    Standard form:
    \[ \dfrac{y}{x} + \left( \frac{1}{x} - 1
    \right) y = \frac{1}{x} \]
    Integrating factor $\mu$:
    \begin{align*}
        \ln \mu(x)
        &= \int^x \left( \frac{1}{u} - 1 \right)
        \dd u \\
        &= \ln x - x \\
        &= \ln(x e^{-x})
    \end{align*}
    \[ \implies \mu(X) = xe^{-x} \]
    Have
    \[ \dfrac{}{x} \left( xe^{-x} y \right)
    = \frac{x}{x} e^{-x} = e^{-x} \]
    Integrate to find
    \[ xe^{-x} y(x) = -e^{-x} + C \]
    \[ y(x) = -\frac{1}{x} + \frac{c}{x} e^x \]
    (if for example we wanted $y$ finite at $x \to
    0$, then we would set $C = 1$.)
\end{example*}

\subsubsection*{Radioactivity example revisited}

Recall
\[ \dfrac{b}{t} + k_b b = k_a a_0 e^{-k_a t} \]
Integrating factor $\mu$:
\[ \ln\mu = \int^t k_b \dd t' = k_b t \]
\[ \implies \mu = e^{k_b t} \]
Have
\[ \dfrac{}{t} ( e^{k_s t} b) = k_a a_0 e^{-k_a t}
e^{k_b t} \]
Consider $k_a = k_s (= k)$
Have
\[ \dfrac{}{t} (e^{kt} b) = k a_0 \]
\[ \implies e^{kt} b(t) = k a_0 t + D \]
\[ \implies b(t) = D e^{-kt} + k a_0 t e^{-kt} \]
Note: $t e^{-kt}$ rather than $e^{-kt}$.

\newpage
\begin{center}
    \LARGE\bfseries
    TOPIC III: NON-LINEAR, FIRST-ORDER ODEs
\end{center}

\noindent
General form
\[ Q(x, y) \dfrac{y}{x} + P(x, y) = 0 \tag{$*$} \]
[Could be non-linear in $\dfrac{y}{x}$ but not
considered here.]

\newpage
\setcounter{section}{0}
\section{Separable Equations}

\begin{definition*}[Separable Equation]
    First-order ODE is separable if can be written
    as
    \[ q(y) \dd y = p(x) \dd x \]
    Solve by direct integration.
    \[ \int q(y) \dd y = \int p(x) \dd x \]
\end{definition*}

\begin{example*}
    \[ (x^2y - 3y) \dfrac{y}{x} - 2xy^2 = 4x \]
    Rearrange:
    \[ y(x^3 - 3) \dfrac{y}{x} = 2x(2 + y^2) \]
    \[ \implies \frac{y}{2 + y^2} \dd y
    = \frac{2x}{x^2 - 3} \dd x \]
    \[ \implies \half \ln(2 + y^2) = \ln(x^2 - 3)
    + C \]
    \[ \implies (y^2 + 2)^\half = A(x^2 - 3) \]
\end{example*}

\section{Exact Equations}

\begin{definition*}[Exact Equation]
    ODE $(*)$ is \emph{exact} if $Q(x, y) \dd y +
    P(x, y) \dd y$ is an exact differential i.e.,
    $\exists f(x, y)$ such that
    \[ \dd f = P(x, y) \dd x + Q(x, y) \dd y \]
    If $(*)$ is exact, then if $\dd f = 0$, then
    $f(x, y) = c$ is a solution. \\
    If $P(x, y) \dd x + Q(x, y) \dd y$ is exact,
    use multivariate chain rule
    \[ \dd f = \pfrac{f}{x} \dd x + \pfrac{f}{y}
    \dd y \]
    so $\exists f(x, y)$ such that
    \[ \pfrac{f}{x} = P \qquad\text{and}\qquad
    \pfrac{f}{y} = Q \]
    Partial derivatives commute so
    \[ \pfrac{}{y} \left( \pfrac{f}{x} \right) =
    \pfrac{P}{y} \qquad\text{and}\qquad
    \pfrac{}{x} \left( \pfrac{f}{y} \right)
    = \pfrac{Q}{x} \]
    are equal. Therefore we find the following
    \emph{necessary} (but not sufficient)
    condition for exactness:
    \[ \text{ODE exact} \implies \pfrac{P}{y} =
    \pfrac{Q}{x} \]
\end{definition*}

\begin{theorem*}
    If $\pfrac{y}{x} = \pfrac{Q}{x}$ throughout a
    \emph{simply-connected} domain $\mathcal{R}$,
    then $P\dd x + Q \dd y$ is an exact
    differential of a single-valued function $f$
    in $\mathcal{R}$.
\end{theorem*}

\begin{definition*}[Simply connected]
    A domain $\mathcal{R}$ is \emph{simply
    connected} if it is path connected (every
    pair of points can be connected by a path in
    $\mathcal{R}$) and any closed curve can be
    shrunk continuously to a point in
    $\mathcal{R}$ without leaving $\mathcal{R}$.
    \begin{center}
    \begin{tsqx}
        CR (-1,0) 0.8 lightblue / black
        CR (1,0) 0.8 lightblue / black
        CR (1,0) 0.4 white / black
        label $\mathcal{R}$ @ (-1,0)
        label $\mathcal{R}$ @ (1,0)
        label Simply connected @ (-1,-1)
        label Not simply connected @ (1,-1)
    \end{tsqx}
    \end{center}
    [see notes for example of $f$ in non-simply
    connected case.]
\end{definition*}

\begin{example*}
    \[ 6y(y - x) \dfrac{y}{x} + 2x - 3y^2 = 0 \]
    Follows that
    \[ \ub{(2x - 3y^2)}_P \dd x + \ub{6y(y - x)}_Q
    \dd y = 0 \]
    Have
    \[ \pfrac{P}{y} = -6y,\qquad \pfrac{Q}{x} =
    -6y \]
    hence \emph{exact} in any simply connected
    domain. Must have
    \[ P = \left. \pfrac{f}{x} \right|_y = 2x -
    3y^2 \]
    \[ Q = \left. \pfrac{f}{y} \right|_x = 6y^2 -
    6xy \]
    \[ \left. \pfrac{f}{x} \right|_y = 2x - 3y^2
    \implies x^2 - 3xy^2 + h(y) = f(x, y) \]
    Follows that
    \[ \left. \pfrac{f}{y} \right|_x = -6xy +
    \dfrac{h}{y} = Q = 6y^2 - 6xy \]
    \[ \implies \dfrac{h}{y} = 6y^2 \implies h(y)
    = 2y^2 + c \]
    Solution of ODE is
    \[ f(x, y) = x^2 - 3y^2x + 2y^2 = C .\]
\end{example*}

\section{Solution Curves and Isoclines}

Graphical methods for ``solving'' of ODEs

\subsection{Solution Curves}

Consider
\[ \dfrac{y}{t} = f(t, y) \]
Each initial condition (e.g. $y(0) = y_0$)
generates a distinct solution curve (trajectory).
\begin{center}
\begin{tsqx}
    ! size(5cm);
    (-0.2,0)->>(1.5,0)
    (0,-0.2)->>(0,1.5)
    label $x$ @ (1.5,-0.15)
    label $y$ @ (0.15,1.5)    
    ! real f(real x) {return 0.7 + 0.6*x*x;}
    ! draw(graph(f,0,1,operator ..), green);
    ! real g(real x) {return 0.3 + 0.4*x*x;}
    ! draw(graph(g,0,1.4,operator ..), red);
    ! dot((0,0.7), green);
    ! dot((0,0.3), red);
\end{tsqx}
\end{center}
Can we sketch these solution curves without
solving ODE?

\begin{example*}
    \[ \dfrac{y}{t} = t(1 - y^2) \]
    [separable, so can solve directly!]
    Separable:
    \[ \frac{1}{1 - y^2} \dd = t \dd t \]
    \[ \implies \half \left( \frac{1}{1 + y} +
    \frac{1}{1 - y} \right) \dd y = t \dd t \]
    \[ \implies \half \left( \ln(1 + y) - \ln(1 -
    y)\right) = \half t^2 + C \]
    \[ \implies \half\ln\left| \frac{1 + y}{1 - y}
    \right| = \half t^2 + C \]
    \[ y = \frac{A - e^{-t^2}}{A + e^{-t^2}} \]
    $A$ parametrises the solution curves. \\
    Initial condition $y(0) = y_0 \implies A =
    \frac{1 + y_0}{1 - y_0}$
    \begin{center}
    \begin{tsqx}
        (-3.5,0)->>(3.5,0)
        (0,-3.5)->>(0,3.5)
        label $y_0$ @ (3.5,-0.3)
        label $A$ @ (0.3,3.5)
        ! real f(real x) {return (1 + x)/(1 - x);}
        ! draw(graph(f,-3.5,0.55,operator, ..), red);
        ! draw(graph(f,1.8,3.5,operator, ..), red);
        (3.5,-1)--(-3.5,-1) dashed
        (1,3.5)--(1,-3.5) dashed
        label $1$ @ (1,0)+(-0.3,0.3)
        label $1$ @ (0,1)+(-0.3,0.3)
        label $-1$ @ (-1,0)+(-0.3,0.3)
        label $-1$ @ (0,-1)+(-0.3,0.3)
    \end{tsqx}
    \end{center}
\end{example*}

\subsection{Slope fields and Isoclines}

\[ \dfrac{y}{t} = f(t, y) \implies
\text{gradient of solution curve through any point} \]
\emph{Slope field} represents gradient field by
short ``sticks'', one centre on each point.

\begin{example*}
    \[ \dfrac{y}{t} = t(1 - y^2) \]
    For $t > 0$, $\dot y > 0$ for $|y| < 1$ and
    $\dot y < 0$ for $|y| > 1$.
\end{example*}

\noindent
\emph{Isoclines} sometimes useful: curves along
which $\dot{y} = \text{const}$.

\begin{example*}
    \[ \dfrac{y}{t} = t(1 - y^2) \]
    Isosclines:
    \[ t(1 - y^2 = D \]
    \[ y^2 = 1 - \frac{D}{t} \]
\end{example*}

Slope field is tangent to solution curves.
[If $f(t, y)$ is single-valued, solution curves
cannot cross.]

\section{Fixed (Equilibrium) Points and Stability}

\begin{definition*}[Fixed / equilibrium points]
    Fixed point of ODE $\dfrac{y}{t} = f(t, y)$ is
    a constant solution $y = c$,
    i.e.
    \[ \dfrac{y}{t} = 0 \,\,\forall t \text{  at
    $y = C$.} \]
\end{definition*}

\begin{example*}
    \[ \dfrac{y}{t} = t(1 - y^2) \]
    Fixed points at $y = \pm 1$.
\end{example*}

\begin{definition*}[Stability of fixed points]
    Fixed point $y = C$ is \emph{stable}
    (\emph{unstable}) if whenever $y$ deviates
    slightly from $C$, $y$ converges (diverges) to
    (from) $y = C$ as $t \to \infty$.
\end{definition*}

\subsection{Perturbation analysis and stability}

Suppose $y = C$ is a fixed point of
\[ \dfrac{y}{t} = f(t, y) \]
[so $f(t, c) = 0\,\,\forall t$] \\
Consider some small perturbation from $y = C$:
\[ y(t) = C + \epsilon(t) \]
so
\begin{align*}
    ODE \implies \dfrac{\epsilon}{t}
    &= f(t, C + \epsilon) \\
    &= f(t, C) + \epsilon \left. \pfrac{f}{y}
    \right|_t (t, C) + O(\epsilon^2)
\end{align*}
Linear for small $\epsilon$:
\[ \dfrac{\epsilon}{t} \approx \epsilon
\pfrac{f}{y} (t, C) \]
(note that this is a linear ODE so should be
easier to solve!) [If $\pfrac{f}{y} = 0$ at the
fixed point, then we need higher-order terms in
Taylor expansion to determine stability.]

\begin{example*}
    \[ \dfrac{y}{t} = t(1 - y^2) \]
    Fixed points at $y = \pm 1$
    \[ \pfrac{f}{y} = -2ty = \begin{cases}
        -2t & \text{at $y = +1$} \\
        2t & \text{at $y = -1$}
    \end{cases} \]
    Near $y = +1$,
    \[ \dfrac{\epsilon}{t} = -2t\epsilon \]
    \[ \implies \epsilon(t) = \epsilon_0 e^{-t^2} \]
    As $t \to \infty$, $\epsilon \to 0 \implies$
    \emph{stable} fixed point. \\
    Near $y = -1$,
    \[ \dfrac{\epsilon}{t} = 2t\epsilon \implies
    \epsilon = \epsilon_0 e^{t^2} \]
    Now $|\epsilon(t)| \to \infty$ as $t \to
    \infty$ so trajectories diverge from $y = -1
    \implies$ \emph{unstable} fixed point.
\end{example*}

\subsection{Autonomous systems and phase portraits}

Special case:
\[ \dfrac{y}{t} = f(y) \]
(so no explicit $t$ dependence) then near a fixed
point $y = C$ ($f(C) = 0$) we have
\[ \dfrac{\epsilon}{t} = \epsilon \dfrac{f}{y}(C)
= k\epsilon \]
\[ \implies \epsilon(t) = \epsilon_0 e^{kt} \]
For autonomous ODE:
\[ f'(c) < 0 \implies \text{stable fixed point} \]
\[ f'(c) > 0 \implies \text{unstable fixed point} \]

\begin{example*}[Chemical kinetics]
    \[ A + B \to C + D \]
    \[ a(t) = a_0 - c(t) \qquad b(t) = b_0 - c(t) \]
    Assume rate of reaction $\alpha$ $a$ $b$ (e.g.
    dilute gas)
    \begin{align*}
        \dfrac{c}{t}
        &= \lambda a(t)b(t) \\
        &= \ub{\lambda [a_0 - c(t)][b_0 - c(t)]}
        _{f(c)}
    \end{align*}
    Fixed points are $C = a_0$ and $C = b_0$.
    Assume $a_0 < b_0$: $c = b_0$ is unphysical.
    \\
    Perturbation analysis:
    \begin{align*}
        \dfrac{f}{c}
        &= \lambda (2c - a_0 - b_0) \\
        &= \begin{cases}
            \lambda (a_0 - b_0) & \text{at $c =
            a_0$} \\
            \lambda (b_0 - a_0) & \text{at $c
            = b_0$}
        \end{cases}
    \end{align*}
    For $a_0 < b_0$, $C = a_0$ is a \emph{stable}
    fixed point (and $C = b_0$ is unstable) \\
    Phase portrait:
    \begin{center}
    \begin{tsqx}
        (-2,0)->>(-1.5,0)
        (-1.5,0)--(-0.5,0)
        (1,0)->>(-0.5,0)
        (1,0)->>(1.5,0)
        (1.5,0)--(2,0)
        ! dot((-1,0), red);
        ! dot((1,0), red);
        label $a_0$ @ (-1,-0.1) +col:red
        label $b_0$ @ (1,-0.1) +col:red
        label $c(t)$ @ (2.2,0)
    \end{tsqx}
    \end{center}
\end{example*}

Since there is no time dependence, we can shift
any solution in time and it will remain a
solution. This is because:
\[ \dfrac{y}{t} = f(y) \implies \int^y \frac{\dd
u}{f(u)} = t + t_0 \]
so if $y(T)$ is a solution, so is $y(t - t_0)
\,\,\forall t_0$.

\begin{example*}[Population dynamics]
    ``Logistic eq'' \\
    Population size $y(t)$ \\
    Birth rate $\alpha y(t)$ with $\alpha > 0$. \\
    Death rate $\beta y + \gamma y^2$
    \begin{align*}
        \dfrac{y}{t}
        &= (\alpha - \beta) y - \gamma y^2 \\
        &= \lambda y (1 - y/\gamma) = f(y)
    \end{align*}
    Autonomous ODE fixed points: $y = 0$ or $y =
    \gamma$.
    \[ \dfrac{f}{y} = \lambda (1 - 2y / \gamma) =
    \begin{cases}
        \lambda & \text{at $y = 0$} \\
        -\lambda & \text{at $y = \gamma$}
    \end{cases} \]
    For $\lambda > 0$, $y = 0$ is unstable fixed
    point; $y = \gamma$ is stable.
    \begin{center}
    \begin{tsqx}
        (0,0)->>(1.5,0)
        (0,0)->>(0,1)
        label $y$ @ (1.6,0)
        label $f(y)$ @ (0,1.1)
        label $0$ @ (0,-0.1)
        label $\gamma$ @ (1.2,-0.1)
        label unstable @ (0,-0.3)
        label stable @ (1.2,-0.3)
        ! real f(real x) {return -2*(x - 0)*(x - 1.2);}
        ! draw(graph(f,0,1.2,operator ..), red);
        (-0.1,-1)->>(0.1,-1)
        (0.1,-1)->>(1.1,-1)
        (1.1,-1)--(1.3,-1)
        (2,-1)->>(1.3,-1)
        label $0$ @ (0,-1.1)
        label $\gamma$ @ (1.2,-1.1)
        ! dot((0,-1), red);
        ! dot((1.2,-1), red);
        label $y(t)$ @ (2.1,-1)
    \end{tsqx}
    \end{center}
    \begin{center}
        \includegraphics[width=0.6\linewidth]
        {images/7153ab9e38a211ec.png}
    \end{center}
\end{example*}

\subsection{Fixed points in discrete equations}

Consider \emph{first order} (links $x_{n + 1}$
only to $x_n$) discrete equation
\[ x_{n + 1} = f(x_n) \]

\begin{definition*}[Fixed point of discrete equation]
    A \emph{fixed point} of a first-order discrete
    equation is a value of $x_n$ such that $x_{n +
    1} = x_n$, i.e.,
    \[ f(x_n) = x_n \]
    Stability now perturbation analysis.
\end{definition*}

Let $x_f$ be a fixed point, write
\[ x_n = x_f + \epsilon_n \]
\begin{align*}
    x_f + \epsilon_{n + 1}
    &= f(x_f + \epsilon_n) \\
    &= f(x_f) + \epsilon_n \dfrac{f}{x} (x_f) +
    O(\epsilon_n^2) \\
    \epsilon_{n = 1}
    &\approx \epsilon_n \dfrac{f}{x} (x_f)
\end{align*}
Follows that
\begin{align*}
    |f'(x_f)| < 1
    &\implies \text{ stable F.P.} \\
    |f'(x_f)| > 1
    &\implies \text{ unstable F.P.}
\end{align*}

\begin{example*}[Logistic map]
    \[ x_{n + 1} = \ub{r x_n(1 - x_n)}_{f(x_n)} \]
    ``Discrete logistic equation'' or ``logistic
    map'' \\
    Complimentary function:
    \[ \dfrac{x}{t} = \lambda x \left(1 -
    \frac{x}{X}\right) \]
    Interested in $x_n \ge 0$. If $r < 5$,
    $\{x_n\}$ stay in range $0 \le x_n \le 1$:
    \begin{center}
    \begin{tsqx}
        (0,0)->>(1.5,0)
        (0,0)->>(0,1)
        label $x_n$ @ (1.6,0)
        label $f(x_n)$ @ (0,1.1)
        ! real f(real x) {return -2*(x - 0)*(x - 1.2);}
        ! draw(graph(f,0,1.2,operator ..), red);
        label $1$ @ (1.2,-0.1)
        label $0$ @ (0,-0.1)
        label $\frac{r}{4}$ @ (-0.1,0.72)
        (-0.05,0.72)--(0.05,0.72)
    \end{tsqx}
    \end{center}
    Fixed points: $x_n = rx_n (1 - x_n)$
    \[ \implies x_n = 0 \wedge x_n = 1 -
    \frac{1}{r} \]
    (note that the second solution only works for
    $x_n > 0$ if $r > 1$) \\
    Stability?
    \[ f'(x) = r(1 - 2\alpha) = \begin{cases}
        r & \text{at $x = 0$} \\
        2 - r & \text{at $x = 1 - \frac{1}{r}$}
    \end{cases} \]
    \emph{Fixed point at $x = 0$}: stable for $0 <
    r < 1$; unstable for $r > 1$. \\
    \emph{Fixed point at $x = 1 - \frac{1}{r}$}:
    stable for $1 < r < 3$; unstable for $r > 3$.
\end{example*}

\noindent
Illustrate with ``cobweb diagrams'' (Details are
non-examinable) \\
$0 < r < 1$:
\begin{center}
    \includegraphics[width=0.6\linewidth]
    {images/2e7f72fa38a511ec.png}
\end{center}
Converge to $x_n = 0$ for any starting points $0
\le x_1 \le 1$.

\newpage

\begin{center}
    \LARGE\bfseries
    TOPIC IV: HIGHER-ORDER LINEAR ODEs
\end{center}
\setcounter{section}{0}

\newpage
\section{Second order ODEs}

General form
\[ \ub{a \dfrac[2]{y}{x} + b \dfrac{y}{x} +
    cy}_{\mathcal{D} y} = f(x)
\tag{$*$} \]
where $\mathcal{D}$ is a \emph{differential
operator} defined by
\[ \mathcal{D} \equiv a \dfrac[2]{}{x} + b
\dfrac{}{x} + c \]
Note that this is \emph{linear}.

\begin{definition*}[Linear operator]
    $\mathcal{D}$ is \emph{linear} if for any
    $y_1(x)$ and $y_2(x)$, and constants $\alpha$
    and $\beta$,
    \[ \mathcal{D}(\alpha y_1 + \beta y_2) =
    \alpha \mathcal{D} y_1 + \beta \mathcal{D} y_2 \]
\end{definition*}

Exploit linearity to solve ODE in 2 steps:
\begin{enumerate}[(1)]
    \item Find complimentary functions that satisfy
        homogeneous equation
        \[ a \dfrac[2]{y_c}{x} + b \dfrac{y_c}{x} +
        cy_c = 0 \]
    \item Find \emph{any} particular integral that satisfies the
        full equation.
\end{enumerate}
Then the solution to full equation from $y_c +
y_p$:
\[ \mathcal{D}(y_c + y_p) = \mathcal{D} (y_c) +
\mathcal{D} (y_p) = 0 + f(x) \]
If $y_{c_1}$ and $y_{c_2}$ are linearly
independent complimentary functions, then
$y_{c_1} + y_p$ and $y_{c_2} + y_p$ are linearly
independent solutions of the full equation.

\begin{definition*}[Linear independence of functions]
    Set of $N$ functions $\{f_i(x)\}$ is
    \emph{linearly dependent} if
    \[ \sum_{i = 1}^N c_i f_i(x) = 0 \]
    for $N$ constants $c_i$, not all of which are
    zero. Otherwise the functions are
    \emph{linearly independent}. [Equivalently, if
    any function can be written as a linear
    combination of the others they are linearly
    dependent.]
\end{definition*}

\subsection{Complementary functions}

Recall
\[ \dfrac{}{x} e^{\lambda x} = \lambda e^{\lambda
x} \]
($\lambda$ is an eigenvalue and $e^{\lambda x}$ is
a eigenfunction) Also $e^{\lambda x}$is an
eigenfunction of $\mathcal{D}$:
\[ \mathcal{D} (e^{\lambda x}) = a \dfrac[2]{}{x}
e^{\lambda x} + b \dfrac{}{x} e^{\lambda x} +
ce^{\lambda x} = \ub{(a\lambda^2 + b\lambda
c)}_{\text{eigenvalue}}
e^{\lambda x} \]
Complimentary functions of $(*)$ satisfy
$\mathcal{D} y_c = 0$: eigenfunctions of
$\mathcal{D}$ with eigenvalue zero:
\[ y_c = A e^{\lambda x} \]
provided $a\lambda^2 + b\lambda + c = 0$.
(\emph{Characteristic equation} of $ay'' + by' +
cy = 0$) \\
Two roots $\lambda_1$ and $\lambda_2$.
\begin{itemize}
    \item If $\lambda_1 \neq \lambda_2$ then we
        have two linearly independent
        complimentary functions
        \[ y_{c_1} \propto e^{\lambda_1 x} \qquad
        \text{and} \qquad y_{c_2} \propto
        e^{\lambda_2 x} \]
        Most general complimentary function is a
        linear combination of $y_{c_1}$ and
        $y_{c_2}$:
        \[ y_c = c_1 y_{c_1}(x) + c_2 y_{c_2} (x) \]
        (so $y_{c_1}$ and $y_{c_2}$ form a basis
        for the solution subspace of the
        homogeneous ODE)
        \begin{note*}
            Roots may be complex $\to$
            oscillations.
        \end{note*}
    \item $\lambda_1 = \lambda_2$: (degenerate
        case). Only generates one linearly
        independent complimentary function of form
        $e^{\lambda_1 x}$. How to find second
        complimentary function? See example below.
\end{itemize}

\begin{example*}[non-degenerate, real roots]
    \[ \dfrac[2]{y}{x} - 5 \dfrac{y}{x} + 6y = 0 \]
    Characteristic equation:
    \[ \lambda^2 - 5\lambda + 6 = 0 \]
    \[ \implies (\lambda - 2)(\lambda - 3) = 0 \]
    Complimentary function:
    \[ y_c(x) = A e^{2x} + B e^{3x} \]
\end{example*}

\begin{example*}[complex roots]
    \[ \dfrac[2]{y}{x} + 4y = 0 \]
    Characteristic equation:
    \[ \lambda^2 + 4 = 0 \implies \lambda = \pm 2i \]
    General complimentary function:
    \begin{align*}
        y_c
        &= A \ub{e^{2i x}}_{\cos 2x + i \sin 2x}
        + B e^{-2i x} \\
        &= \ub{(A + B)}_{\alpha} \cos 2x + \ub{i(A
        - B)}_{\beta} \sin 2x
    \end{align*}
\end{example*}

\begin{example*}[degeneracy and detuning]
    \[ \dfrac[2]{y}{x} - 4 \dfrac{y}{x} + 4y = 0 \]
    Characteristic equation:
    \[ \lambda^2 - 4\lambda + 4 = 0 \]
    \[ \implies (\lambda - 2)^2 = 0 \]
    Degenerate: $\lambda = 2, 2$ so only
    \emph{one} linearly independent solution of
    the form $e^{\lambda x}$. Now \emph{detuning}
    is removing degeneracy by considering
    ``detuned'' (modified) equation:
    \[ \dfrac[2]{y}{x} - 4 \dfrac{y}{x} + (4 -
    \epsilon^2) y = 0 \qquad (\epsilon \ll 1) \]
    in $\lim_{\epsilon \to 0}$, detuned equation
    $\to$ one we really want to solve. Try $y =
    e^{\lambda x}$:
    \[ \lambda^2 - 4 \lambda + (4 - \epsilon^2) =
    0 \]
    \[ \implies \lambda = 2 \pm \epsilon \]
    General solution of detuned equation is
    \begin{align*}
        y
        &= A e^{(2 + \epsilon) x} + B e^{(2 -
        \epsilon) x} \\
        &= e^{2x} [ A e^{\epsilon x} + B
        e^{-\epsilon x} ] \\
        &= e^{2x} [(A + B) + \epsilon (A - B) + (A
        \epsilon^2) + O(B \epsilon^2)]
    \end{align*}
    Consider applying initial conditions
    \[ y(0) = C \qquad \text{and} \qquad y'(0) = D \]
    (to original detuned equation)
    Have
    \[ A + B = C \qquad \text{and} \qquad 2(A + B)
    + \epsilon (A - B) = D \]
    so
    \[ y \to e^{2x} [C + (D - 2C)x + O(A
    \epsilon^2) + O(B\epsilon^2)] \]
    We have
    \begin{align*}
        A
        &= \half \left( C + \frac{D -
        2C}{\epsilon} \right) \\
        \implies A \epsilon^2
        &= \half \left( C\epsilon^2 + (D -
        2C)\epsilon \right) \\
        \implies O(A\epsilon^2)
        &= O(\epsilon)
    \end{align*}
    It follows that the general solution of
    original equation is
    \[ y = e^{2x} [\alpha + \beta x] \]
    for arbitrary constants $\alpha$ and $\beta$.
\end{example*}

General rule: if $y_{c_1}(x)$ is a degenerate
complimentary function of a linear ODE with
constant coefficients, then
\[ y_{c_2}(x) = x y_{c_1}(x) \]
is a second, linearly independent complimentary
function.

\newpage
\section{Homogeneous 2nd-order ODEs with non-constant coefficients}

General form
\[ y'' + p(x) y' + q(x) y = 0 \tag{$*$} \]

\subsection{Second CF --- reduction of order}

Aim: given one solution of $(*)$, $y_2(x)$, find
another, $y_2(x)$. \\
Look for a solution
\[ y_2 = v(x) y_1(x) \]
So we have
\begin{align*}
    y_2
    &= v'y_2 + vy_1' \\
    y_2''
    &= v''y_1 + 2v'y_1' + vy_1''
\end{align*}
If $y_2$ satisfies $(*)$ then
\[ v''y_1 + v'(2y_1 + py_1) + v(\ub{y_1'' + py_1' +
qy_1}_{\text{$= 0$ by $y_1$ in $(*)$}}) = 0 \]
Follows that
\[ v''y_1 + v'(2y_1' + py_1) = 0 \]
\[ u'y_1 + u(2y_1' + py_1) = 0 \qquad \text{for $u
= v'$} \]
(note that this is a first order ODE in $u$!)
Separable:
\[ \frac{u'}{u} = -\frac{y_1'}{y_1} - p(x) \]
\[ \implies \ln u = -2 \ln y_1 - \int_0^x p(w) \dd w
+ \ln A \]
\[ \implies u(x) = \frac{A}{[y_1(x)]^2} \exp \left[
-\int_0^x p(w) \dd w \right] \]
Integrate again to find $v(x)$.

\begin{example*}
    \[ y'' \ub{- 4}_{p(x)}y' \ub{+ 4}_{q(x)}y = 0 \]
    Have a degenerate solution $y_1 = e^{2x}$.
    Reduction of order:
    \[ \frac{u'}{u} = -2 \times 2 + 4 = 0 \]
    \[ \implies v' = \text{const} \implies v = Ax
    + B .\]
    Have $y_2(x) = (Ax + B)e^{2x}$. So $xe^{2x}$
    is a second linearly independent solution.
\end{example*}

\subsection{Phase space}

For $n$-th order linear ODE
\[ p(x) y^{(n)} + q(x) y^{(n - 1)} + \cdots +
r(x)y = f(x) \]
From this equation, $y^{(n)}(x)$ is determined by
$y(x)$, $y'(x)$, \dots, $y^{(n - 1)}(x)$. By
differentiating the equation, we can also find
that higher derivatives are determined. Hence we
can compute all the derivatives at a point, and
therefore construct the Taylor series about $x_0$
if $y(x_0)$, $y'(x_0$, \dots, $y^{(n - 1)}(x)$ are
specified. So solution is determined by these
initial conditions. \\
State of system fully specified (can predict
subsequence evalutation) at any $x$ by
\emph{solution vector}
\[ \bf{Y}(x) =
\begin{pmatrix}
    y(x) \\
    y'(x) \\
    \vdots \\
    y^{(n - 1)}(X)
\end{pmatrix}
\]
At any $x$, $\bf{Y}(x)$ defines a point in $nD$
\emph{phase space}.

\begin{example*}
    \[ y'' = 4y = 0 \tag{SHM} \]
    \[ y_1 = \cos 2x; \qquad y_2 = \sin 2x \]
    \[ \bf{Y}_1(x) =
    \begin{pmatrix}
        y_1 \\
        y_1'
    \end{pmatrix}
    =
    \begin{pmatrix}
        \cos 2x \\
        -2\sin 2x
    \end{pmatrix}
    \]
    \[ \bf{Y}_2(x) =
    \begin{pmatrix}
        y_2 \\
        y_2'
    \end{pmatrix}
    =
    \begin{pmatrix}
        \sin 2x \\
        2 \cos 2x
    \end{pmatrix}
    \]
    Phase space
    \begin{center}
    \begin{tsqx}
        ! size(5cm);
        (-1.5,0)->>(1.5,0)
        (0,-1.5)->>(0,1.5)
        label $y$ @ (1.6,0)
        label $y'$ @ (0,1.6)
        unitcircle red
        y_2 := dir(70)
        y_1 := dir(-30)
        origin->>y_2 blue
        origin->>y_1 blue
        label $\mathbf{Y}_2(x)$ @ 1.2*y_2 +col:blue
        label $\mathbf{Y}_1(x)$ @ 1.3*y_1 +col:blue
    \end{tsqx}
    \end{center}
\end{example*}

\begin{note*}
    $\bf{Y}_1$ and $\bf{Y}_2$ are linearly
    independent vectors $\forall x$ (in linear
    algebra sense).
\end{note*}

\noindent
Can use $\bf{Y}_1$ and $\bf{Y}_2$ at any $x$ as a
\emph{basis} for phase space.

\subsection{Wronskian and linear dependence}

Recall, $y_i(x)$ ($i = 1, \dots, n$) are linearly
depdendent if
\[ \sum_{i = 1}^n c_i y_i(x) = 0 \qquad \forall x \]
must hold for all $x$ hence differentiable $n - 1$
times:
\[ \sum_{i = 1}^n c_i \bf{Y}_i(x) = \bf{0} \]
See that
\[ \{y_i(x)\} \text{ linearly dependent} \implies
\{\bf{Y}_i\} \text{ linearly depdendent $\forall
x$} \]
\emph{Fundamental matrix}: $(\bf{Y}_1, \bf{Y}_2,
\dots, \bf{Y}_n)$ has zero determinant if the
$\{y_i\}$ are linearly dependent.

\begin{definition*}[Wronskian]
    Wronskian $W(x)$ of $n$ functions $y_i(x)$ ($i
    = 1, \dots, n$) is
    \[ W(x) = |\bf{Y}_1, \bf{Y}_2, \bf{Y}_3,
    \dots, \bf{Y}_n| =
    \begin{pmatrix}
        y_1 & y_2 & \cdots & y_n \\
        y_1' & y_2 & \cdots & y_n' \\
        \vdots & \vdots & \ddots & & \\
        y_1^{(n - 1)} & y_2^{(n - 1)} & & y_n^{(n - 1)}
    \end{pmatrix}
    \]
\end{definition*}
We have
\[ \{y_i(x)\} \text{ linearly dependent}
\implies W(x) = 0 \,\, \forall x \]
(test for linear dependence of $n$ functions)
\begin{remark*}
    $W(x) = 0$ does not necessarily imply
    linear dependence of the $y_i(x)$.
\end{remark*}

\begin{example*}
    \[ y'' + 4y = 0 \implies y_1 = \cos 2x, y_2 =
    \sin 2x \]
    \[ W(x) =
    \left| \begin{matrix}
        y_1 & y_2 \\
        y_1' & y_2'
    \end{matrix} \right|
    =
    \left| \begin{matrix}
        \cos 2x & \sin 2x \\
        -2\sin 2x & 2\cos 2x
    \end{matrix} \right|
    = 2 \]
    \[ W \neq 0 \implies \text{$y_1$ and $y_2$ are
    linearly independent.} \]
\end{example*}

\subsection{Abel's Theorem}

\begin{theorem*}[Abel's Theorem]
    Given any two solutions of
    \[  y'' + p(x)y + q(x)y = 0 \]
    If $p(x)$ and $q(x)$ are continuous on an
    interval $I$, then either $W(x) = 0
    \,\,\forall x \in I$ or $W(x) \neq 0
    \,\,\forall x \in I$.
\end{theorem*}

\begin{proof}[(sketch)]
    \begin{align*}
        W
        &= y_1y_2' - y_2y_1' \\
        \implies W'
        &= y_1y_2'' - y_2y_1'' \\
        &= y_2(py_1' + qy_1) - y_1(py_2' + qy_2)
        \\
        &= y_1py_1' - y_1py_2' \\
        &= -p(x)W
    \end{align*}
    Separable ODE for $W$:
    \[ W(x) = W(x_0) \ub{\exp \left( -\int_{x_0}^x
    p(u) \dd u \right)}_{\text{never zero}} \]
    (This is known as ``Abel's identity'') If
    $W(x)$ vanished at some point $x_0$, it is
    zero $\forall x$; else is is never zero.
\end{proof}

\begin{note*}
    If $p(x) = 0$ (no $y'$ term in ODE), then the
    Wronskian is constant.
\end{note*}

\begin{example*}[Bessel's Equation]
    \[ x^2 y'' + xy' + (x^2 - n^2)y = 0 \]
    generally no closed form solutions (but for
    example for some half integer values there
    exists a closed form solution)
    \[ y'' + \frac{1}{x} y' + \left( 1 -
    \frac{n^2}{x^2} \right) y = 0 \]
    Abel:
    \begin{align*}
        W(x)
        &= W(x_0) \exp \left( -\int_{x_0}^x \frac{\dd
        u}{u} \right) \\
        &= W(x_0) \frac{x_0}{x}
    \end{align*}
    Note: Abel's identity determines form of $W$
    without having to solve ODE.
\end{example*}

Geometric interpretation: (phase space)
\begin{center}
\begin{tsqx}
    ! size(5cm);
    (-1.5,0)->>(1.5,0)
    (0,-1.5)->>(0,1.5)
    label $y$ @ (1.6,0)
    label $y'$ @ (0,1.6)
    yscale(1.3)*unitcircle red
    origin->>yscale(1.3)*dir(-80) blue
    label $\mathbf{Y}_1$ @ yscale(1.3)*1.5*dir(-80) +col:blue
    origin->>yscale(1.3)*dir(60) blue
    label $\mathbf{Y}_2$ @ yscale(1.3)*1.4*dir(60) +col:blue
\end{tsqx}
\end{center}
Solution vectors always collinear or never
collinear as $x$ varies. [Basis at some $x \to$
basis $\forall x$.]

\subsubsection*{Applications of Abel's Theorem}

\[ y_1y_2' - y_2y_1' = W(x_0) \exp \left(
-\int_{x_0}^x p(w) \dd w \right) \]
If $y_1$ known, can solve this ODE to find $y_2$:
\[ \dfrac{}{x} \left( \frac{y_2}{y_1} \right) =
\frac{W(x_0)}{y_1^2(x)} \exp \left( -\int_{x_0}^x
p(u) \dd u \right) \]
Exactly as for ``reduction of order''.

\subsubsection*{Generalisation}

Abel's Theorem also holds for solutions of
$n$th-order linear, homogeneous ODEs. \\
Such an equation can always be written as a system
of first-order equations:
\[ \bf{Y}' + \ul{\bf{A}} \bf{Y} = \bf{0} \]

\begin{example*}
    \[ y''' + p(x)y'' + q(x)y' + r(x)y = 0 \]
    \[
    \begin{pmatrix}
        y \\
        y' \\
        y''
    \end{pmatrix}
    +
    \ub{
    \begin{pmatrix}
        0 & -1 & 0 \\
        0 & 0 & -1 \\
        r & q & p
    \end{pmatrix}
    }_{\ul{\bf{A}}}
    \begin{pmatrix}
        y \\
        y' \\
        y''
    \end{pmatrix}
    = \bf{0} \]
\end{example*}

Can be shown (Q7, Ex 3) that
\begin{align*}
    W' + \mathrm{Tr}(\ul{\bf{A}}) W
    &= 0 \\
    \implies W(x)
    &= W(x_0) \exp \left( -\int_{x_0}^x
    \mathrm{Tr} [\ul{\bf{A}}(u)] \dd u \right)
\end{align*}
so Abel's theorem holds.

\subsection{Linear equidimensional ODEs}

\begin{definition*}[Equidimensional ODE]
    A linear, second-order ODE is
    \emph{equidimensional} if it has the following
    form
    \[ ax^2y'' + bxy' + cy = f(x) \tag{$*$} \]
    (where $a$, $b$ and $c$ are constants)
\end{definition*}

\noindent
Why call this equidimensional? \\
Let $y(x)$ be a complimentary function of $(*)$.
Consider
\[ \phi(x) = y(\lambda x) \]
for some arbitrary real $\lambda$. Then
\[ x\dfrac{\phi}{x} = x\lambda y'(\lambda x) \]
\[ x^2 \dfrac[2]{\phi}{x} = x^2 \lambda^2 y''
(\lambda x) \]
\[ \implies ax^2 \dfrac[2]{\phi}{x} + bx
\dfrac{\phi}{x} + c\phi = \text{LHS of $(*)$
evaluated at $\lambda x$} = 0 \]

\subsubsection*{Solving by Eigenfunctions}

$y = x^k$ is an eigenfunction of $x \dfrac{}{x}$,
with eigenvalue $k$:
\[ x \dfrac{}{x} x^k = kx^{k} \]
Look for complimentary functions of $(*)$: $y_c
= x^k$
\[ ak(k - 1) + bk + c = 0 \]
\[ \implies ak^2 + (b - a)k + c = 0 \]
Solve quadratic equation to determine $k_1$ and
$k_2$. If $k_1 \neq k_2$ (non-degenerate) general
complimentary function:
\[ y_c = A x^{k_1} + B x^{k_2} \]

\subsubsection*{Solving by substitution}

Switch independent variable from $x$ to
\[ z \equiv \ln x \]
For $y(x) = y(e^z)$
\[ \dfrac{y}{z} = \ub{e^z y'(e^z)}_{xy'(x)} \]
\[ \dfrac[2]{y}{z} = \ub{e^z y'(e^z)}_{x'(x)} +
\ub{e^{2z} y''(e^z)}_{x^2 y''(x)} \]
Follows that $(*)$ becomes
\[ \ub{a \dfrac[2]{y}{z} + (b - a) \dfrac{y}{z} +
cy}_{\text{constant coefficients}} = f(e^z) \]
Complimentary function $y_c \propto e^{\lambda z}$
\[ a\lambda^2 + (b - a)\lambda + c = 0 \]
(Some characteristic equation as for $k_1$ and
$k_2$)
General complimentary function ($k_1 \neq k_2$):
\[ y_c = A e^{k_1 z} + B e^{k_2 z} = A x^{k_1} +
B x^{k_2} \]
Degenerate case ($k_1 = k_2 = k$)
\[ y_c = A e^{kz} + Bx e^{kz} = A x^k + B \ln(x)x^k \]

\newpage
\section{Inhomogeneous (Forced) 2nd-order ODEs}

\subsection{Particular Integrals of equations with constant coefficients}

\[ a \dfrac[2]{y}{x} + b \dfrac{y}{x} + cy = f(x) \]
\begin{center}
    \begin{tabular}{c|c}
        $f(x)$ & $y_p(x)$ \\
        \hline
        $e^{mx}$ & $A e^{mx}$ \\
        $\sin kx$ or $\cos kx$ & $A \sin kx + B
        \cos kx$ \\
        Polynomial $p_n(x)$ ($n$-th degree) &
        Polynomial $q_n(x)$ ($a_nx^n + \cdots +
        a_1x + a_0$)
    \end{tabular}
\end{center}

\begin{example*}
    \[ y'' - 5y' + 6y = \ub{2x + e^{4x}}_{f(x)} \]
    Try
    \begin{align*}
        y_p
        &= \ob{(Ax + B)}^{\text{for $2x$}} +
        \ob{ce^{4x}}^{\text{for $e^{4x}$}} \\
        y_p'
        &= A + 4ce^{4x} \\
        y_p''
        &= 16ce^{4x}
    \end{align*}
    \[ \ub{(16C - 20C + 6C)}_{\to 1}e^{4x} +
    \ub{(6A)}_{\to 2}x + \ub{(-5A + 6B)}_{\to
    0} = 2x + e^{4x} \]
    \[ C = \half,\quad A = \frac{1}{3},\quad B =
    \frac{5}{18} \]
    Complementary function: $\alpha e^{3x} + \beta
    e^{2x}$. General solution:
    \[ \alpha e^{3x} + \beta e^{2x} + \half e^{4x}
    + \frac{1}{3}x + \frac{5}{18} \]
\end{example*}

\subsubsection*{Resonance}

What if forcing term involves a complimentary
function? $\to$ Detuning.

\begin{example*}
    \[ \ddot{y} + \omega_0^2 y = \sin \omega_0 t
    \tag{$*$} \]
    Complimentary function:
    \[ y_c(t) = A \cos \omega_0t + B \sin
    \omega_0t \tag{SHM} \]
    Oscillator driven ``resonantly''. \\
    ``Detuned equation'':
    \[ \ddot{y} + \omega_0^2 y = \sin\omega t
    \qquad (\omega \neq \omega_0) \tag{\cross} \]
    Try particular integral:
    \[ y_p(t) = C \sin \omega t \]
    Substitute into (\cross):
    \begin{align*}
        C(\omega_0^2 - \omega^2)
        &= 1 \\
        \implies C
        &= \frac{1}{\omega_0^2 - \omega^2}
    \end{align*}
    Limit as $\omega \to \omega_0$ does not exist.
    Add in a complimentary function to regularise
    the limit $\omega \to \omega_0$:\
    \[ y_p(t) = \frac{1}{\omega_0^2 - \omega^2}
    \left( \sin \omega t - \ub{\sin \omega_0
    t}_{\text{CF}} \right) \]
    Evaluate indeterminate limit with
    L'H\^opital's rule:
    \[ \lim_{\omega \to \omega_0} y_p(t) =
    \lim_{\omega\to \omega_0} \left( \frac{t
    \cos \omega t}{-2\omega} \right) =
    -\frac{t \cos \omega_0 t}{2\omega_0} \]
    Particular integral of original equation ($*$)
    is
    \[ y_p(t) = -\frac{t}{2\omega_0} \cos\omega_0
    t \]
\end{example*}

\noindent
General rule: if forcing is a linear combination
of complimentary functions, particular integral is
of the form
\[ y_p(t) = t \times \ub{(\text{non-resonant
PI})}_{\text{$\cos\omega_0 t$ above}} \]

\subsubsection*{Resonance in Equidimensional ODEs}

Complimentary function for equidimensional ODE
(non-degenerate):
\[ y_c = Ax^{k_1} + B x^{k_2} \]
If forcing $f(x) \propto x^{k_1}$ (or $x^{k_2}$),
particular integral is
\[ y_p(x) \propto \ln x \times x^{k_1} \]
Follows from transforming to $z = \ln x$:
\[ y_p(z) \propto z e^{k_1z} = \ln x \,\,x^{k_1} \]

\subsection{Variation of Parameters}

Systematic way of finding a particular integral
given two linearly independent complimentary
functions ($y_1$ and $y_2$).

\bigskip

\noindent
Consider
\[ y'' + p(x) y' + q(x) y = f(x) \tag{$*$} \]
with linearly independent complimentary functions
$y_1$ and $y_2$. Solution vectors $\bf{Y}_1(x)$
and $\bf{Y}_2(x)$. We will use these as a
\emph{basis} in phase space at any argument $x$ to
write the solution vector of the particular
integral as
\[ \bf{Y}_p(x) = u(x)\bf{Y}_1(x) + v(x)\bf{Y}_2(x) \]
In components:
\[ y_p = u(x) y_1 = v(x) y_2 \tag{\dag} \]
\[ y_p' = u(x)y_1' + v(x) y_2' \tag{\ddag} \]
\[ \dfrac{}{x} y_p' = y_p''
\ub{=}_{(\text{\ddag})} uy_1'' + u'y_1' + vy_2'' =
v'y_2' \]
Follows that
\[ f(x) = uy_1'' + u'y_1' + vy_2'' + v'y_2'
+ p(x)[uy_1' + vy_2'] + q(x)[uy_1 + vy_2] \]
Now using the fact that $y_1$ and $y_2$ are
complimentary functions, we get that
\[ f(x) = u'y_1' + v'y_2' \]
However (\ddag) has to be consistent with
$\dfrac{}{x}$(\dag), so
\[ \dfrac{}{x} y_p = u'y_1 + uy_1' + v'y_2 + vy_2' \]
Compare to (\ddag):
\[ u'y_1 + v'y_2 = 0 \]
Combining with $u'y_1' + v'y_2' = f(x)$, we have
\[
\begin{pmatrix}
    y_1 & y_2 \\
    y_1' & y_2'
\end{pmatrix}
\begin{pmatrix}
    u' \\
    v'
\end{pmatrix}
=
\begin{pmatrix}
    0 \\
    f
\end{pmatrix}
\]
\[ \implies
\begin{pmatrix}
    u' \\
    v'
\end{pmatrix}
= \frac{1}{\ub{W}_{y_1y_2' - y_2y_1'}}
\begin{pmatrix}
    y_2' & -y_2 \\
    -y_1' & y_1
\end{pmatrix}
\begin{pmatrix}
    0 \\
    f
\end{pmatrix}
\]
Hence
\[ u' = -\frac{y_2}{W}f;\quad v' = \frac{y_1}{W}f \]
Integrate to find
\[ y_p(x) = y_2(x) \int^x \frac{1}{W(t)}f y_1(t)
f(t) \dd t - y_1(x) \int^x \frac{1}{W(t)} y_2(t)
f(t) \dd t \]
(changing the lower limit in these integrals just
adds a multiple of a complimentary function).

\begin{example*}
    \[ y'' + 4y = \ub{\sin 2x}_{f(x)} \]
    Complimentary functions $y_1 = \sin 2x$, $y_2
    = \cos 2x$. (Note that forcing term is
    resonant).
    \[ W(x) = -2 \]
    so
    \[ y_p(x) = \cos 2x \int^x -\half \ub{\sin 2t
    \sin 2t}_{\half(1 - \cos 4t)} \dd t
    - \sin 2x \int^x - \half \ub{\cos 2t \sin
    2t}_{\half \sin 4t} \dd t \]
    \begin{align*}
        y_p(x)
        &= - \frac{1}{4} \cos 2x \left[ x -
        \frac{1}{4} \sin 4x \right]
        + \frac{1}{4} \sin 2x \left[ -\frac{1}{4}
        \cos 4x \right] \\
        &= \frac{1}{16} \sin 2x - \frac{1}{4} x
        \cos 2x
    \end{align*}
\end{example*}

\newpage
\section{Forced Oscillating Systems: Transients and Damping}

\begin{center}
    \includegraphics[width=0.6\linewidth]
    {images/260251aa748811ec.png}
\end{center}
Newton II:
\begin{align*}
    m \ddot{y}
    &= \sum \text{forces} \\
    &= -ky - b \dot{y} + F(t)
\end{align*}
\[ m \ddot{y} + b \dot{y} + ky = F(t) \]
($b$, $k$ positive constants). For $b = 0$ and
$F(T) = 0$, simple harmonic motion at
\[ \omega_0 = \sqrt{\frac{k}{m}} \]
Dimensionless $t$ coordinate:
\[ \tau = \omega_0 t \qquad \left( \dfrac{}{\tau}
= 1 \right) \]
Divide ODE by $k$:
\[ y'' + \ub{2 \kappa y'}_{\frac{b}{m\omega_0}} +
y = \ub{f(\tau)}_{\frac{F(t)}{k}} \]
Unforced system described by one dimensionless
parameter $\kappa$.

\subsection{Free (unforced or natural) response}

$f = 0$:
\[ y'' + 2\kappa y' + y = 0 \]
Look for solutions $y \propto e^{\lambda x}$:
\[ \lambda^2 + \kappa \lambda + 1 = 0 \]
(characteristic equation).
\[ \lambda_1, \lambda_2 = -\kappa \pm
\sqrt{\kappa^2 - 1} \]

\subsubsection*{Light damping (underdamping): $\kappa < 1$}

$\lambda_1$ and $\lambda_2$ complex:
\[ \lambda_1, \lambda_2 = -\kappa \pm i \sqrt{1 -
\kappa^2} \]
General solution:
\[ y(\tau) = e^{-\kappa \tau} [A \sin(\sqrt{1 -
\kappa^2} \tau)] + B \cos(\sqrt{1 - \kappa^2}
\tau) \]
($A$ and $B$ constants)
\begin{center}
    \includegraphics[width=0.6\linewidth]
    {images/72a7c14c748911ec.png}
\end{center}

\subsubsection*{Critical damping: $\kappa = 1$}

Degenerate case:
\[ \lambda_1 = \lambda_2 = -\kappa \]
General solution:
\[ y(\tau) = e^{-\kappa \tau}(A + B\tau) \]
($A$ and $B$ constants)
\begin{center}
    \includegraphics[width=0.6\linewidth]
    {images/b25fc0a0748911ec.png}
\end{center}

\subsubsection*{Heavy damping (overdamping): $\kappa > 1$}

$\lambda_1$ and $\lambda_2$ are real: take
$|\lambda_2| > |\lambda_1|$ without loss of
generality.
\[ \lambda_1 = \ub{-\kappa + \sqrt{\kappa^2 -
1}}_{\text{negative}} \qquad
\lambda_2 = \ub{-\kappa - \sqrt{\kappa^2 -
1}}_{\text{positive}} \]
General solution:
\[ y(\tau) = A e^{-|\lambda_1|\tau} + B
e^{-|\lambda_2|\tau} \]
(first term dominates long-term motion if
present)
\begin{center}
    \includegraphics[width=0.6\linewidth]
    {images/2eff4de2748a11ec.png}
\end{center}
Note: unforced response decays eventually in all
cases.

\subsection{Forced response}

Initial ``transient'' response from PI + CF, but
CF decays leaving ``steady-state'' response (PI).

\begin{example*}
    \[ \ddot{y} + \mu \dot{y} + \omega_0^2 y =
    \frac{F_0}{m} \sin \omega t \tag{$*$} \]
    Assuming light damping ($\mu < 2\omega_0$).
    Complimentary function:
    \[ y_c(t) = e^{-\frac{\mu t}{2}} [A \sin
    \omega_{\text{free}} t + B \cos
    \omega_{\text{free}} t] \]
    ($\kappa = \frac{\mu}{2\omega_0}$) where
    \[ \omega_{\text{free}} = \sqrt{\omega_0^2 -
    \frac{\mu^2}{4}} \]
    For PI: try
    \[ y_p(t) = \frac{F_0}{m} (C \sin \omega t + D
    \cos \omega t) \]
    Substitute into ($*$):
    \[ \sin \omega t: \qquad -\omega^2 C -
    \mu\omega D + \omega_0^2 C = 1 \]
    \[ \cos \omega t: \qquad \omega^2 D + \mu
    \omega C + \omega_0^2 D = 0 \]
    Gives
    \[ y_p(t) = \frac{F_0/m}{(\omega_0^2 -
    \omega^2)^2 + (\mu \omega)^2} [(\omega_0^2 -
    \omega^2) \sin \omega t - \mu \omega \cos
    \omega t] \]
    \begin{center}
        \includegraphics[width=0.6\linewidth]
        {images/5d45dc4c748b11ec.png}
    \end{center}
    Amplitude of $y_p$:
    \[ \frac{F_0/m}{\sqrt{(\omega_0^2 -
    \omega^2)^2 + (\mu \omega)^2}} \]
    ($\frac{F_0}{m}\sqrt{c^2 + p^2}$)
    \begin{center}
        \includegraphics[width=0.6\linewidth]
        {images/110670c0748c11ec.png}
    \end{center}
\end{example*}

\newpage
\section{Impulses and Point Forces}

Consider sudden force $F(t)$ applied between $t =
T - \epsilon$ and $T + \epsilon$.
\begin{center}
\begin{tsqx}
    ! size(5cm);
    (-0.2,0)->>(1.5,0)
    (0,-0.2)->>(0,1.5)
    ! real f(real x) {return 1-16*(x-0.75)*(x-0.75);}
    ! draw(graph(f,0.5,1,operator..), red);
    (0.25,0)--(0.5,0) red
    (1,0)--(1.25,0) red
    label $T$ @ (0.75,-0.1)
    label $T - \epsilon$ @ (0.5,-0.1)
    label $T + \epsilon$ @ (1,-0.1)
    (0.75,0.05)--(0.75,-0.01)
    label $t$ @ (1.6,0)
    label $F(t)$ @ (0,1.6)
\end{tsqx}
\end{center}
For example strike on oscillator. \\
As $\epsilon \to 0$, only the \emph{impulse}
\[ I = \int_{T - \epsilon}^{T + \epsilon} F(t) \dd
t \]
matters for subsequent motion $\to$
(mathematically) convenient to consider limit of a
sudden impulse.

\bigskip

\noindent
Forced, damped oscillator
\[ m \ddot{y} + b\dot{Y} + ky = F(t) \]
Integrate ODE from $T - \epsilon$ to $T +
\epsilon$
\[ \lim_{\epsilon \to 0} \left(
m[\dot{y}]_{T - \epsilon}^{T + \epsilon} +
\ub{b[y]_{T - \epsilon}^{T + \epsilon}}_{\text{0
if $y$ continuous}} + \ub{k \int_{T -
\epsilon}^{T + \epsilon} y \dd t}_{\text{0 if $y$
remains finite}} \right) = I \]
here $I$ denotes impulse. Follows that
\[ \lim_{\epsilon \to 0} m[\hat{y}]_{T -
\epsilon}^{T + \epsilon} = I \]
velocity is discontinuous.
\begin{center}
\begin{tsqx}
    (-0.1,0)--(0,0)
    (0.5,0)->>(1.5,0)
    (0,-0.75)->>(0,0.75)
    label $t$ @ (1.6,0)
    label $y(t)$ @ (0,0.85)
    label $T$ @ (0.5,-0.1)
    (0.5,-0.05)--(0.5,0.05)
    (0,0)--(0.5,0) red
    ! real f(real x) {return 0.25*Sin(1000*(x-0.5))*(1/x)*(1/x);}
    ! draw(graph(f,0.5,1.5,operator..), red);
\end{tsqx}
\end{center}

\subsection{Dirac delta function}

Consider a family of functions $D(t; \epsilon)$
with
\begin{enumerate}[1.]
    \item $\lim_{\epsilon\to 0} D(t; \epsilon) = 0
        \qquad \forall t \neq 0$.
    \item $\int_{-\infty}^\infty D(t, \epsilon)
        \dd t$.
\end{enumerate}
Impulse force considered scalar $F(t) = ID(t - T,
\epsilon)$. \\
Example family:
\[ D(t; \epsilon) =
\frac{e^{-t^2/\epsilon^2}}{\epsilon\sqrt{\pi}} \]
(Q14 on example sheet 1)
\begin{center}
\begin{tsqx}
    (-1.5,0)->>(1.5,0)
    (0,-0.1)->>(0,1.5)
    label $t$ @ (1.6,0)
    label $D(t; \epsilon)$ @ (0,1.6)
    ! real f(real x) {return 1.3*(10)^(-x**2);}
    ! draw(graph(f,-1.5,1.5,operator..), red);
    ! real g(real x) {return 0.8*(10)^(-0.3*x**2);}
    ! draw(graph(g,-1.5,1.5,operator..), blue);
\end{tsqx}
\end{center}
Family not unique, but for any such family limit
$\epsilon \to 0$ yields \emph{Dirac delta
function}.

\begin{definition*}[Dirac delta function]
    The Dirac delta function is defined by
    \[ \delta(t) \equiv \lim_{\epsilon \to 0}
    D(t; \epsilon) .\]
\end{definition*}
Should only use inside integrals. \\
Key properties:
\begin{enumerate}[1.]
    \item $\delta(x) = 0 \qquad \forall x \neq 0$.
    \item $\int_{-\infty}^\infty \delta(x) \dd x =
        1$.
    \item For all $g(x)$ continuous at $x = 0$
        \begin{align*}
            \int_{-\infty}^\infty g(x) \delta(x)
            \dd x \\
            &= \lim_{\epsilon \to 0}
            \int_{-\infty}^\infty g(x)D(x;
            \epsilon) \dd x \\
            &= g(0)
        \end{align*}
        Generalises to ($b > a$)
        \[ \int_a^b g(x) \delta(x - x_0) \dd x =
        \begin{cases}
            g(x_0) & \text{if $a < x_0 < b$} \\
            0 & \text{otherwise}
        \end{cases} \]
\end{enumerate}

\subsection{Delta-function forming}

Consider
\[ y'' + p(x)y' + q(x)y = \delta(x) \]
For $x < 0$ or $x > 0$ we have
\[ y'' + p(x)y' + q(x)y = 0 \]
However, discontinuity in $y'$ at $x = 0$:
integrate ODE from $-\epsilon$ to $+\epsilon$
\[ \lim_{\epsilon \to 0} [y']_{-\epsilon}^\epsilon + p(0)
\ub{\lim_{\epsilon \to 0}
[y]_{-\epsilon}^\epsilon}_{\text{0 if $y$
continuous}} +
\ub{\lim_{\epsilon \to 0} \int_{-\epsilon}^\epsilon
q(x) y \dd x}_{\text{0 if $y$ is finite}} = 1 \]
\[ \implies \lim_{\epsilon \to 0}
[y']_{-\epsilon}^\epsilon = 1 \]
``Jump condition''

\begin{note*}
    Continuity of $y$ at $x = 0$ is required to
    avoid $y' \sim \delta(x)$ around $x = 0$ and
    $y''$ even worse behaved.
\end{note*}
\noindent
General rule: highest-order derivative term in the
ODE addresses the delta-function forcing.

\begin{example*}
    \[ y'' - y = 3\delta \left( x - \frac{\pi}{2} \right) \]
    boundary conditions $t = 0$ at $x = 0$ and $y
    = 0$ at $x = \pi$. \\
    For $0 \le x < \frac{\pi}{2}$: $y'' - y = 0$
    \[ y = A\sinh(x) \qquad [y(0) = 0] \]
    For $\frac{\pi}{2} < x \le \pi$: $y'' - y = 0$
    \[ y = C\sinh(\pi - x) \qquad [y(\pi) = 0] \]
    Join up at $x = \frac{\pi}{2}$ with
    \begin{enumerate}[1.]
        \item We have $y$ continuous $\implies A = C$.
        \item Now ``Jump condition'':
            \[ \lim_{\epsilon \to 0}
            [y']_{\frac{\pi}{2}-\epsilon}^{\frac{\pi}{2}
            + \epsilon} = 3 \]
            \[ \implies -C \cosh \left(
            \frac{\pi}{2} \right) - A \cosh \left(
            \frac{\pi}{2} \right) = 3 \]
            hence
            \[ A = C = -\frac{3}{2\cosh \left(
            \frac{\pi}{2} \right) } \]
            \begin{center}
            \begin{tsqx}
                (-0.2,0)->>(1.5,0)
                label $x$ @ (1.6,0)
                (0,-1)->>(0,0.2)
                label $y$ @ (0,0.3)
                label $\frac{\pi}{2}$ @ (0.7,0.1)
                label $\pi$ @ (1.4,0.1)
                (0.7,-0.05)--(0.7,0.05)
                (1.4,-0.05)--(1.4,0.05)
                ! real f(real x) {return (0-1.3)*x**2;}
                ! draw(graph(f,0,0.7,operator..), red);
                ! real g(real x) {return (0-1.3)*(x-1.4)**2;}
                ! draw(graph(g,0.7,1.4,operator..), red);
            \end{tsqx}
            \end{center}
            (discontinuity in $y'$ at the
            ``spike'')
    \end{enumerate}
\end{example*}

\subsection{Heaviside Step Function}

\begin{definition*}[Heaviside function]
    The \emph{Heaviside step function} is
    \[ H(x) = \int_{-\infty}^x \delta(t) \dd t \]
    Follows that
    \[ t(x) = \begin{cases}
        0 & \text{for $x < 0$} \\
        1 & \text{for $x > 0$} \\
        \text{undefined} & \text{for $x = 0$}
    \end{cases} \]
\end{definition*}

\begin{center}
\begin{tsqx}
    (0,0)->>(1,0)
    label $x$ @ (1.1,0)
    (0,-0.1)->>(0,1.2)
    label $H(x)$ @ (0,1.3)
    (0,1)--(1,1) red
    (-1,0)--(0,0) red
    ! dot((0,1), red, filltype=FillDraw(fillpen=white));
    ! dot((0,0), red, filltype=FillDraw(fillpen=white));
\end{tsqx}
\end{center}
Fundamental theorem of calculus implies
\[ \dfrac{H}{x} = \delta(x) \]

\subsubsection*{Forcing with $H(x)$}

Consider
\[ y'' + p(x)y' + q(x)y = H(x) \]
Have
\[ y'' + p(x)y' + q(x)y = \begin{cases}
    0 & \text{for $x < 0$} \\
    1 & \text{for $x > 0$}
\end{cases} \]
Follows that
\[ \lim_{\epsilon \to 0}
[y'']_{-\epsilon}^\epsilon + p(0) \lim_{\epsilon
\to 0} [y']_{-\epsilon}^\epsilon + q(0)
[y]_{-\epsilon}^\epsilon = 1 \tag{$*$} \]
If $y'' \sim H(x)$ around $x = 0$, $y'$ and $y$
are continuous and ($*$) is satisfied.
\[ \lim_{\epsilon \to 0} [y']_{-\epsilon}^\epsilon
= 0 \quad \text{and} \quad \lim_{\epsilon \to 0}
[y]_{-\epsilon}^\epsilon = 0 \]
``Jump'' conditions for Heaviside forcing.

\newpage
\section{Higher Order Discrete Equations}

Consider a linear, discrete second-order equation
\[ a y_{n + 2} + by_{n + 1} + cy_n = f_n \]
where $a$, $b$ and $c$ are constants and $f_n$ is
a forcing sequence.

\begin{note*}
    Might arise from discretising a second order
    ODE with constant coefficients.
    \[ \left. \dfrac[2]{y}{x} \right|_{x_n} 
    \simeq \frac{\ob{y(x_n + h)}^{y_{n+1}} +
    \ob{y(x_n - h)}^{y_{n-1}} -
    \ob{2y(x_n)}^{y_n}}{h^2} \]
\end{note*}

\noindent
Solve with similar methods to ODEs:
\[ y_n = y_n^{(c)} + y_n^{(p)} \]
Complementary function satisfies
\[ ay_{n+2}^{(c)} + by_{n+1}^{(c)} + cy_n^{(c)} =
0 \tag{$*$} \]
For ODE with constant coefficients, $y$ or
$e^{\lambda x} \to e^{\lambda nk} \o k^n$. \\
Try $y_n^{(c)} \propto k^n$ in ($*$):
\begin{align*}
    ak^{n + 2} + bk^{n + 1} + ck^n
    &= 0 \\
    \implies ak^{2} + bk + c
    &= 0 \qquad k \neq 0
\end{align*}
``Characteristic equation''. \\
Two roots of characteristic equation in general,
$k = k_1$ and $k = k_2$. So general complimentary
function:
\[ y_n^{(c)} = \begin{cases}
    Ak_1^n + Bk_2^n & \text{if $k_1 \neq k_2$} \\
    (A + Bn)k_1^n & \text{if $k_1 = k_2$}
\end{cases} \]
(Note that this is similar to differential
equations, where we get $A + Bx$).

\bigskip

\noindent
Now we will look at common forcing functions, and
suitable particular integrals to try in these
cases.
\begin{center}
    \begin{tabular}{c|c}
        $f_n$ & $y_n^{(c)}$ \\
        \hline
        $k^n$ & $Ak^n \text{if $k \not\in \{k_1, k_2\}$}$ \\
        $k_1^n$ & $Ank_1^n$ \\
        $n^p \text{($p$ a non-negative integer)}$ & $An^p + Bn^{p - 1} + \cdots + Cn + D$
    \end{tabular}
\end{center}

\begin{example*}[Fibonacci sequence]
    \[ y_n = y_{n - 1} + y_{n - 2} \]
    with initial conditions $y_0 = 1$ and $y_1 =
    1$. So the first few terms are $1, 1, 2, 3,
    5, 8$. We have
    \[ y_{n + 2} - y_{n + 1} - y_n = 0 \tag{$*$} \]
    Try $y_n = k^n$:
    \[ \implies k^2 - k - 1 = 0 \]
    \[ \implies k = \frac{1 \pm \sqrt{5}}{2} \]
    ``Golden ratio'' $\phi_1 = \frac{1 +
    \sqrt{5}}{2}$ and $\frac{-1}{\phi_1} = \frac{1
    - \sqrt{5}}{2}$.
    General solution of ($*$) is
    \[ y_n = A\phi_1^n + B\phi_2^n \]
    for some constants $A$ and $B$, which are
    determined by the initial conditions:
    \[ y_0 = 1 \implies A + B = 1 \]
    \[ y_1 = 1 \implies A\phi_1 + B\phi_2 = 1 \]
    \[ \implies A = \frac{\phi_1}{\sqrt{5}}, B =
    \frac{-\phi_2}{\sqrt{5}} =
    \frac{1}{\sqrt{5}\phi_1} \]
    Follows that
    \[ y_n = \frac{\phi_1^{n + 1} - \phi_2^{n + 1}}{\sqrt{5}}
    = \frac{\phi_1^{n + 1} - (-1 / \phi_1)^{n +
1}}{\sqrt{5}} \]
    Note that $y_n$ is a sequence of integers, but
    $\phi_1$ is irrational, so it's almost magical
    that the formula always gives an integer. \\
    Since $\phi_1 > 1$, $y_n \to \phi_1^{n + 1} /
    \sqrt{5}$ as $n \to \infty$, so
    \[ \lim_{n \to \infty} \frac{y_{n + 1}}{y_n} =
    \phi_1 .\]
\end{example*}

\newpage
\section{Series Solutions}

Series solutions for linear, homogeneous,
second-order ODEs.
\[ p(x)y'' + q(x)y' + r(x)y = 0 \tag{$*$} \]
Seek a series solution around $x_0$. Feasibility
depends on nature of $p(x)$, $q(x)$ and $r(x)$
around $x_0$.

\subsection{Classification of Singular Points}

\begin{definition*}[Ordinary and Singular Points]
    Point $x_0$ is an \emph{ordinary point} if
    both $q(x) / p(x)$ \emph{and} $r(x) / p(x)$
    have Taylor series around $x_0$. (``analytic''
    there). Otherwise, $x_0$ is a \emph{singular
    point}.
\end{definition*}

\noindent
Types of singular point: if $x_0$ is a singular
point of ($*$), but ($*$) can be rewritten as
\[ \ub{P(x)(x - x_0)^2}_{p(x)} y'' + \ub{Q(x) (x -
x_0)}_{q(x)} y' + \ub{R(x)}_{r(x)} y = 0 .\]
with $\frac{Q(x)}{P(x)}$ and $\frac{R(x)}{P(x)}$
analytic at $x_0$, then $x_0$ is a \emph{regular
singular point}. Otherwise $x_0$ is an
\emph{irregular singular point}. \\
Equivalently: if
\[ (x - x_0) \frac{q(x)}{p(x)} \left[ =
\frac{Q}{P} \right] \]
and
\[ (x - x_0)^2 \frac{r(x)}{p(x)} \left[ =
\frac{R}{P} \right] \]
are analytic at $x_0$, then $x_0$ is a regular
singular point. \\
``No non-singular and equidimensional equation.''

\begin{example*}
    \[ \ub{(1 - x^2)}_{p}y'' - \ub{2x}_{q}y' +
    \ub{2}_{r}y = 0 \]
    \[ \frac{q(x)}{p(x)} = \frac{-2x}{1 - x^2} =
    \frac{2x}{(x - 1)(x + 1)} \]
    \[ \frac{r(x)}{p(x)} = \frac{2}{1 - x^2} =
    -\frac{2}{(x - 1)(x + 1)} \]
    $x = \pm 1$ are \emph{singular} points. They
    are \emph{regular} since, for example
    \[ (x - 1) \frac{q(x)}{p(x)} = \frac{2x}{x + 1} \]
    analytic at $x = 0$.
\end{example*}

\begin{example*}
    \[ \ub{\sin x}_{p} y'' + \ub{\cos x}_{q} y' +
    \ub{2}_{r} y = 0 \]
    We have
    \[ \frac{q}{p} = \frac{\cos x}{\sin x} \]
    \[ \frac{r}{p} = \frac{2}{\sin x} \]
    $x = n\pi$ are singular points ($n$ integer).
    However
    \begin{align*}
        \ub{(x - n\pi)}_{\eps} \frac{q}{p}
        &= \eps \frac{\omega s(n\pi + \eps)}
        {\sin(n\pi + \eps)} \\
        &= \eps \frac{\cos \eps}{\sin \eps}
    \end{align*}
    (Has a Taylor series about $\eps = 0$). Also
    \[ (x - n\pi)^2 \frac{r}{p} =
    \frac{2\eps^2}{(-1)^n \sin \eps} \]
    which has a Taylor series about $\eps = 0$.
    Hence $x = n\pi$ are \emph{regular singular
    points}.
\end{example*}

\begin{example*}
    \[ \ub{(1 + \sqrt{x})}_{p} y'' \ub{-2}_{q} y'
    + \ub{2}_{r} y = 0 \]
    \[ x \frac{q}{p} = - \frac{2x^2}{1 + \sqrt{x}} \]
    does not have a Taylor series about $x = 0$,
    so $x = 0$ is an \emph{irregular singular
    point}.
\end{example*}

\subsection{Method of Frobenius}

\[ p(x) y'' + q(x) y' + r(x) y = 0 \tag{$*$} \]

\begin{theorem*}[Fuch's Theorem]
    \phantom{.}
    \begin{enumerate}[1.]
        \item If $x = x_0$ is an ordinary point of
            ($*$), there are two linearly
            independent solutions of the form
            \[ y = \sum_{n = 0}^\infty a_n (x -
            x_0)^n \]
            (Taylor series convergent in some
            neighbourhood of $x_0$).
        \item If $x_0$ is a regular singular point
            of ($*$) there is at least one
            solution of the form
            \begin{align*}
                y
                &= \sum_{n = 0}^\infty a_n (x -
                x_0)^{n + \sigma} \\
                &= (x - x_0)^\sigma \sum_{n =
                0}^\infty a_n(x - x_0)^n
            \end{align*}
            for some real $\sigma$ with $a_0 \neq
            0$ (this is known as a ``Frobenius
            series''). No guarantee that we obtain
            two linearly independent solutions in
            this case.
    \end{enumerate}
    \noindent
    Method may fail completely for irregular
    singular points.
\end{theorem*}

% \begin{example*}[ordinary point]
\noindent\textbf{Example of ordinary point}
    \[ (1 - x^2) y'' - 2x y' + 2y = 0 \]
    Expand around $x = 0$ (ordinary point).
    [$x = \pm 1$ are regular singular points].
    Try
    \begin{align*}
        y
        &= \sum_{n = 0}^\infty a_n x^n \\
        y'
        &= \sum_{n = 1}^\infty na_n x^{n - 1} \\
        y''
        &= \sum_{n = 2}^\infty n(n - 1)a_n x^{n - 2}
    \end{align*}
    Convenient to multiply ($*$) by $x^2$:
    \[ (1 - x^2)(x^2 y'') - 2x^2(xy') + 2x^2(y) =
    0 \]
    Substituting:
    \[ \sum_{n = 2}^\infty a_n [(1 - x)^2n(n - 1)]
    x^n - 2 \sum_{n = 1}^\infty a_n x^2n x^n + 2
    \sum_{n = 0}^\infty a_n x^{n + 2} = 0 \]
    which can be rewritten as
    \[ \sum_{n = 2}^\infty a_n[(1 - x^2)n(n -
    1)]x^n - 2 \sum_{n = 2}^\infty a_{n - 2} (n -
    2) x^n  + 2 \sum_{n = 0}^\infty a_n x^{n + 2}
    = 0 \]
    Equate coefficients of $x^n$ for $n \ge 2$:
    \[ a_n n(n - 1) - a_{n - 2}(n - 2)(n - 3) -
    2a_{n - 2}(n - 2) + 2a_{n - 2} = 0 \tag{$n \ge
    2$} \]
    \[ \implies n(n - 1) a_n = n(n - 3)a_{n - 2} \]
    \[ \implies a_n = \frac{n - 3}{n - 1}a_{n - 2} \]
    (recursive relation). $a_0$ and $a_1$ are not
    fixed (arbitrary constants). Since $a_3 = 0
    \implies a_5 = 0 \implies a_7 = 0$ etc, there
    is one odd solution:
    \[ y = a_1 x \]
    other, even solution, for $n$ even:
    \begin{align*}
        a_n
        &= \frac{n - 3}{n - 1} a_{n - 2} \\
        &= \left( \frac{n - 3}{n - 1} \right)
        \left( \frac{n - 5}{n - 3} \right) a_{n -
        4} \\
        &= \cdots \\
        &= -\frac{1}{n - 1} a_0
    \end{align*}
    \[ \implies y = a_0 \left[ 1 - x^2 -
    \frac{x^4}{3} - \frac{x^4}{5} - \cdots \right] \]
    Note:
    \[ \ln(1 \pm x) = \pm x - \frac{x^2}{2} \pm
    \frac{x^3}{3} - \cdots \]
    so
    \[ y = a_0 \left[ 1 - \frac{x}{2} \ln \left(
    \frac{1 + x}{1 - x} \right) \right] \]
    Closed-form solution - note behaviour near $x
    = \pm 1$ (regular singular points); see below.
% \end{example*}

\bigskip

% \begin{example*}[regular singular point]
\noindent\textbf{Example of Regular Singular Point}
    \[ 4xy'' + (1 - x^2)y' - xy = 0 \tag{$*$} \]
    $x = 0$: regular singular point \\
    Try
    \[ y(x) = \sum_{n = 0}^\infty a_n x^{n +
    \sigma} \qquad (a_0 \neq 0) \]
    Multiply ($*$) by $x$:
    \[ 4(x^2y'') + 2(1 - x^2)(xy') + x^2y = 0 \]
    Hence,
    \[ \sum_{n = 0}^\infty a_n x^{n + \sigma}(4(n +
    \sigma)(n + \sigma - 1) + 2(1 - x^2)(n +
    \sigma) - x^2) = 0 \]
    $\sigma$ determined by lowest power of $x$. \\
    \ul{$x^\sigma$}: $a_0(4\sigma(\sigma - 1) +
    2\sigma) = 0$ (note $a_0 \neq 0$ by
    assumption). Hence $2\sigma(2\sigma - 1) = 0$
    so $\sigma = 0$ or $\sigma = \half$. \\
    Next lowest power of $x$: \\
    \ul{$x^{\sigma + 1}$}: $a_1(4(\sigma +
    1)\sigma + 2(1 + \sigma)) = 0$ hence
    $a_1(\sigma + 1)(2\sigma + 1) = 0$. Since we
    know that $\sigma = 0$ or $\half$, we must
    have $a_1 = 0$. \\
    \ul{$x^{n + \sigma}$} ($n \ge 2$): $a_n(4(n +
    \sigma)(n + \sigma - 1) + 2(1 + \sigma)) -
    2a_{n - 2}(n + \sigma - 2) - a_{n - 2} = 0$ so
    \[ 2(n + \sigma)(2n + 2\sigma - 1)a_n = (2n +
    2\sigma - 3)a_{n - 2} \tag{\dag} \]
    Consider roots separately: \\
    \ul{$\sigma = 0$}: (\dag) gives $2n(2n - 1)a_n
    = (2n - 3)a_{n - 2}$
    \[ \implies a_n = \frac{2n - 3}{2n(2n -
    1)}a_{n - 2} \]
    \[ a_2 = \frac{1}{4 \times 3} a_0; \qquad a_4 =
    \frac{5}{8 \times 7} a_2 = \frac{5 \times 1}{8
    \times 7 \times 4 \times 3} a_0 \]
    ($a_1 = 0$, so $a_3 = 0$ etc) \\
    Solution:
    \[ y = a_0 \left( 1 + \frac{x^2}{4 \times 3} +
    \frac{5x^4}{8 \times 7 \times 4 \times 3} +
    \cdots \right) \]
    (Note that this is a Taylor series since the
    root $\sigma$ was an integer) \\
    \ul{$\sigma = \half$} (\dag) gives $(2n +
    1)2na_0 = (2n - 2)a_{n - 2}$ for $n \ge 2$,
    hence
    \[ a_n = \frac{n - 1}{n(2n + 1)}a_{n - 2} \]
    \[ a_2 = \frac{1}{2 \times 5}; \qquad a_4 =
    \frac{3 \times 1}{4 \times 9 \times 2 \times 5} \]
    Solution:
    \[ y = b_0 x^{\half} \left( 1 + \frac{x^2}{2
    \times 5} + \frac{3x^4}{4 \times 9 \times 2
    \times 5} + \cdots \right) \]
    (Frobenius series).
% \end{example*}

\noindent
Note that we have determined 2 linearly
independent solutions, but this is not
generally the case.

\subsection{Second solutions}

Expansion around a regular singular point
guarantees one Frobenius solution. Whether we get
a second depends on roots $\sigma_1$ and
$\sigma_"$ of indicial equation.
\begin{enumerate}[1.]
    \item $\sigma_1 - \sigma_2$ \emph{not} an
        integer, then we get two linearly
        independent solutions:
        \[ y_1 = (x - x_0)^{\sigma_1} \sum_{n =
        0}^\infty a_n(x - x_0)^n \]
        \[ y_2 = (x - x_0)^{\sigma_2} \sum_{n =
        0}^\infty b_n(x - x_0)^n \]
        Note: as $x \to x_0$, general solution
        $\sim a_0(x - x_0)^{\sigma_1}$ if
        $\sigma_1 < \sigma_2$.
    \item $\sigma_1 - \sigma_2$ is a non-zero
        integer. One series solution involving
        \emph{larger} root (say, $\sigma_2$).
        \[ y_1 = (x - x)6)^{\sigma_2} \sum_{n =
        0}^\infty a_n(x - x_0)^n \]
        Second solution of ODE is generally of the
        form
        \[ y_2 = (x - x_0)^{\sigma_1} \sum_{n =
        0}^\infty b_n(x - x_0)^n +cy_1(x) \ln(x -
        x_0) \]
        ($c$ may or may not be zero) \\
        Constant $c$ is determined by $a_1$ (in
        $y_1$) and $b_0 \to$ general solution
        depends on two arbitrary constants as
        required.
    \item $\sigma_1 = \sigma_2 = \sigma$ log term
        is \emph{always} required, $c \neq 0$.
        \[ y_1 = (x - x_0)^\sigma \sum_{n =
        0}^\infty a_n(x - x_0)^n \]
        \[ y_2 = (x - x_0)^\sigma \sum_{n =
        0}^\infty b_n(x - x_0)^n + cy_1(x) \ln(x -
        x_0) \]
\end{enumerate}

\begin{example*}[case 2]
    \[ x^2y'' - xy = 0 \tag{$*$} \]
    $x = 0$ is a regular singular point. \\
    Look for
    \[ y = \sum_{n = 0}^\infty a_n x^{n + \sigma}
    \qquad (a_0 \neq 0) \]
    ($*$) implies
    \[ \sum_{n = 0}^\infty (a_n(n + \sigma)(n +
    \sigma - 1)x^{n + \sigma} - a_n x^{n + \sigma
    + 1} ) = 0 \]
    Lowest power is $x^\sigma$: $a_0\sigma(\sigma
    - 1) = 0$ hence $\sigma_1 = 0$ and $\sigma_2 =
    1$ (since $a_0 \neq 0$). \\
    Since $\sigma_1 - \sigma_2$ is a non-zero
    integer, this is case 2. \\
    For $n \ge 1$.
    \[ a_n (n + \sigma)(n + \sigma - 1) = a_{n -
    1} \]
    \ul{$\sigma = 1$} (larger root):
    \[ a_n = \frac{a_{n - 1}}{n(n + 1)} \implies
    a_n = \frac{a_0}{n!(n + 1)!} \]
    Solution
    \[ y_1 = a_0 x \left( 1 + \frac{x}{2} +
    \frac{x^2}{12} + \cdots \right) \]
    \ul{$\sigma = 0$} (smaller root):
    \[ n(n - 1)a_n = a_{n - 1} \]
    but $n = 1$ gives $a_0 = 0$, which
    contradicts initial assumption. So we need a
    log term:
    \[ y_2(x) = \ub{\sum_{n = 0}^\infty b_n
    x^n}_{\text{$\sigma_1 = 0$ series}} +
    cy_1(x)\ln x \]
    Determine $\{b_n\}$ and $c$ by direct
    substitution in ODE or by method of reduction
    or order. (See non-examinable section of
    lecture notes)
\end{example*}

\newpage

\begin{center}
    \bfseries\LARGE
    TOPIC V: MULTIVARIATE FUNCTIONS - APPLICATIONS
\end{center}
\setcounter{section}{0}

\begin{itemize}
    \item Directional derivatives
    \item Extrema
    \item Coupled systems of first-order ODEs
    \item Partial differential equations
\end{itemize}

\newpage

\section{Directional Derivative}

Consider function $f(x, y)$ and vector
displacement $\bf{\dd s} = (\dd x, \dd y)$.
Infinitesimal change in $f$:
\begin{align*}
    \dd f
    &= \pfrac{f}{x} \dd x + \pfrac{f}{y} \dd y \\
    &= (\dd x, \dd y) \cdot \left( \pfrac{f}{x},
    \pfrac{f}{y} \right) \\
    &= \bf{\dd s} \cdot \nabla f
\end{align*}
$\nabla f$ denotes the ``gradient of $f$''.

\begin{definition*}[Gradient of $f$]
    \[ \nabla f = \left( \pfrac{f}{x},
    \pfrac{f}{y} \right) \tag{$*$} \]
\end{definition*}

\begin{remark*}
    This is vector valued, and sometimes instead
    denoted $\mathrm{grad}\,\,f$.
\end{remark*}

\noindent
If we write
\[ \bf{\dd s} = \dd s \hat{\bf{s}} \]
where $\dd s$ is distance displacement and
$\hat{\bf{s}}$ is a unit vector ($|\hat{\bf{s}}| = 1$).

\begin{definition*}[Directional derivative]
    The \emph{directional derivative} of $f$ in
    direction $\hat{\bf{s}}$ is
    \[ \dfrac{f}{s} = \hat{\bf{s}} \cdot (\nabla f) \]
    Rate of change of $f$ with distance $s$ along
    direction $\hat{\bf{s}}$.
\end{definition*}

Can then \emph{define} gradient vector $\nabla
f$ geometrically as that vector such that
\[ \dfrac{f}{s} = \hat{\bf{s}} \cdot (\nabla t) \]

\bigskip

\noindent
\textbf{Properties of gradient vector}

\begin{center}
\begin{tsqx}
    ! size(5cm);
    f := (1,0.8)
    s := (0.2,0.9)
    label $\hat{\mathbf{s}}$ @ 1.1*s +col:red
    label $\nabla f$ @ 1.1*f +col:blue
    anglemark4 f origin s 5
    origin->>f blue
    origin->>s red
    label $\theta$ @ 0.2*dir(60)
\end{tsqx}
\end{center}

\[ \dfrac{f}{s} = \hat{\bf{s}} \cdot \nabla f =
\cos\theta |\nabla f| \]

\begin{enumerate}[1.]
    \item Direction $\nabla f$ is that in which
        $f$ \emph{increases} most rapidly.
    \item Magnitude of $\nabla f$ is maximum rate
        of change of $f$:
        \[ |\nabla f| = \max \left( \dfrac{f}{s} \right) \]
    \item If $\hat{\bf{s}}$ is parallel to
        contours of $f$, then
        \[ 0 = \dfrac{f}{s} = \hat{\bf{s}} \cdot
        (\nabla f) \]
        Hence, $\nabla f$ is $\perp$ contours of
        $f(x, y)$.
\end{enumerate}

\newpage
\section{Stationary Points}

Always at least one direction for which
$\dfrac{f}{s} = 0$ (parallel to contours of $f$).
\\
\emph{Stationary points}
\[ \dfrac{f}{s} = 0 \forall \hat{\bf{s}} \]
Since
\[ \dfrac{f}{s} = \hat{\bf{s}} \cdot (\nabla f) \]
we must have $\nabla f$ at stationary points.

\bigskip

\noindent
\textbf{Types of stationary points} \\
Local maxima:
\begin{center}
    \includegraphics[width=0.6\linewidth]
    {images/a6de09384b8711ec.png}
\end{center}
Contours of $f$ are locally elliptical. \\
Local minima:
\begin{center}
    \includegraphics[width=0.6\linewidth]
    {images/f39313e04b8711ec.png}
\end{center}
Contours are locally elliptical. \\
Saddle points:
\begin{center}
    \includegraphics[width=0.6\linewidth]
    {images/57e56b864b8811ec.png}
\end{center}
Contours locally hyperbolic. Contours cross at and
only at saddle points.

\newpage
\section{Classification of Stationary Points}

Consider how $f(x, y)$ varies along line
\[ \bf{x}(s) = \bf{x}_0 + s \hat{\bf{s}} \]
\begin{center}
    \includegraphics[width=0.6\linewidth]
    {images/0a37b54a4b8a11ec.png}
\end{center}
Can think of $f(x(s), y(s))$ as a function of $s$
so usual single-variable Taylor series holds:
\begin{align*}
    f(\bf{x}_0 + s \hat{\bf{s}})
    &= f(\bf{x}_0) + s \left. \dfrac{f}{s}
    \right|_{\bf{x}_0} + \half s^2 \left.
    \dfrac[2]{f}{s} \right|_{\bf{x}_0} + \cdots \\
    &= f(\bf{x}_0) + \ub{s \hat{\bf{s}} \cdot \nabla f
    |_{\bf{x}_0}}_{1} + \ub{\half s^2 (\hat{\bf{s}} \cdot
    \nabla) (\hat{\bf{s}} \cdot \nabla f)}_{2} + \cdots
\end{align*}
Define vector displacement
\[ \delta \bf{x} = s \hat{\bf{s}} \]
components
\[ \delta x = x(s) - x_0, \qquad \delta y = y(s) -
y_0 \]
Now we focus on terms 1 and 2:
\begin{enumerate}[1.]
    \item \eqnoskip
        \[ s \hat{\bf{s}} \cdot (\nabla f) =
        \bf{\delta x} \cdot \nabla  = \delta x
        \dfrac{f}{x} + \delta y \dfrac{f}{y} \]
    \item \eqnoskip \begin{align*}
        s^2(\hat{\bf{s}} \cdot
        \ol{\bf{v}})(\hat{\bf{s}} \cdot \nabla f)
        &= (\bf{\delta x} \cdot \nabla f) \\
        &= \left( \delta x \pfrac{}{x} + \delta y
        \pfrac{}{y} \right) \left( \delta x
        \pfrac{f}{x} + \delta y \pfrac{f}{y}
        \right) \\
        &= (\delta x)^2 \pfrac[2]{f}{x} + \delta x
        \delta y \frac{\partial^2 f}{\partial x
        \partial y} + \delta y \delta x
        \frac{\partial^2 f}{\partial x \partial y}
        + (\delta y)^2 \pfrac[2]{f}{y} \\
        &= (\delta x, \delta y)
        \begin{pmatrix}
            f_{xx} & f_{xy} \\
            f_{yx} & f_{yy}
        \end{pmatrix}
        \begin{pmatrix}
            \delta x \\
            \delta y
        \end{pmatrix}
    \end{align*}
\end{enumerate}

\begin{definition*}[Hessian matrix]
    The \emph{Hessian matrix} is
    \[ \ul{\bf{H}} =
    \begin{pmatrix}
        f_{xx} & f_{xy} \\
        f_{yx} & f_{yy}
    \end{pmatrix}
    = \nabla \nabla f \]
    Symmetric since $f_{xy} = f_{yx}$.
\end{definition*}

Multivariate Taylor series
\begin{align*}
    f(x_0 + \delta x, y_0 + \delta y)
    &= f(x_0, y_0) + \left. \left( \delta x \pfrac{f}{x}
    + \delta y \pfrac{f}{y} \right) \right|_{x_0, y_0} 
    \\
    &\,\,\,+ \left. \half \left( (\delta x)^2
    \pfrac[2]{f}{x} + 2 \delta x \delta y
    \frac{\partial^2}{\partial x
    \partial y} + (\delta y)^2 \frac{\partial^2
    f}{\partial y^2} \right) \right|_{x_0, y_0} +
    \cdots
\end{align*}
Coordinate free form:
\[ f(\bf{x}_0 + \bf{\delta x}) + f(\bf{x}_0) +
\delta x \cdot \nabla f |_{x_0} + \half \bf{\delta
x} (\nabla \nabla f)|_{x_0} \bf{\delta x}^2 +
\cdots \]

\subsection{Nature of stationary points and the Hessian}

Let $\bf{x}_0$ be a stationary point ($\bf{\nabla}
f \mid_{\bf{x}_0} = \bf{0}$). Around $\bf{x}_0$:
\[ f(\bf{x}) \approx f(\bf{x}_0) + \half
\ub{\bf{\delta x}}_{\bf{x} - \bf{x}_0}
\ub{\ul{\bf{H}}(\bf{x}_0)}_{\text{Hessian}}
\bf{\delta x}^\top \]
Nature of stationary point depends on properties
of $\ul{\bf{H}}(\bf{x}_0)$.

\begin{definition*}[Positive-definite and negative definite matrices]
    A (real) symmetric matrix $\ul{\bf{H}}$ is
    positive definite if
    \[ \bf{x} \ul{\bf{H}} \bf{x}^\top > 0 \]
    for all real $\bf{x}$ ($\bf{x} \neq \bf{0}$).
    It is negative definite if
    \[ \bf{x} \ul{\bf{H}} \bf{x}^\top < 0 \]
    for all real $\bf{x}$ ($\bf{x} \neq \bf{0}$).
    Otherwise, indefinite.
\end{definition*}

If $\ul{\bf{H}}(\bf{x}_0)$ is positive definite,
then $\bf{\delta x} \ul{\bf{H}} \bf{\delta
x}^\top > 0$ for all $\bf{\delta x}$ so $f(\bf{x})
> f(\bf{x}_0)$ in vicinity of $\bf{x}_0$:
\begin{center}
    $\ul{\bf{H}}$ positive definite $\implies$
    local minimum
\end{center}
Similarly,
\begin{center}
    $\ul{\bf{H}}$ negative definite $\implies$
    local maximum
\end{center}
If indefinite, may be a maximum, minimum or a
saddle.

\subsubsection*{Definiteness and eigenvalues}

$\ul{\bf{H}}$ is real symmetric $\to$ diagonalise
$\ul{\bf{H}}$ with an orthogonal transformation.
In $n$D, use coordinates along principal axes.
\[ \bf{\delta x} \ul{\bf{H}} \bf{\delta x}^\top =
(\delta x_1, \delta x_2, \dots, \delta x_n)
\begin{pmatrix}
    \lambda_1 & 0 & \cdots & 0 \\
    0 & \lambda_2 & \cdots & 0 \\
    \vdots & \vdots & \ddots & \vdots \\
    0 & 0 & \cdots & \lambda_n
\end{pmatrix}
\begin{pmatrix}
    \delta x_1 \\
    \delta x_2 \\
    \vdots \\
    \delta x_n
\end{pmatrix}
\]
Follows that
\begin{center}
    $\ul{\bf{H}}$ positive definite $\iff$ all
    $\lambda_i > 0$
\end{center}
\begin{center}
    $\ul{\bf{H}}$ negative definite $\iff$ all
    $\lambda_i < 0$
\end{center}
If all $\lambda_i$ are non-zero, but of mixed
sign, have a saddle point. \\
If ant eigenvalues are 0, need higher terms in
Taylor expansion.

\begin{example*}
    \[ f(x, y) = x^2 + y^4 \]
    Global \emph{minimum} at $(0, 0)$
    \[ \ul{\bf{H}} =
    \begin{pmatrix}
        2 & 0 \\
        0 & 12y^2
    \end{pmatrix}
    =
    \begin{pmatrix}
        2 & 0 \\
        0 & 0
    \end{pmatrix}
    \qquad \text{at $(0, 0)$} \]
    \[ \implies \lambda_1 = 2, \lambda_2 = 0 \]
\end{example*}

\subsubsection*{Definiteness and signature of $\ul{\bf{H}}$}

Alternative to establish if positive or negative
definite.

\begin{definition*}[Signature of $\ul{\bf{H}}$]
    The \emph{signature} of $\ul{\bf{H}}$ is
    pattern of signs of the ordered determinants
    of the leading principal minors of
    $\ul{\bf{H}}$: for example for $f(x_1, x_2,
    \dots, x_n)$
    \[ \ub{| f_{x_1x_1} |}_{|\ul{\bf{H}}_1|},
    \ub{\left|
    \begin{matrix}
        f_{x_1x_1} & f_{x_1x_2} \\
        f_{x_2x_1} & f_{x_2x_2}
    \end{matrix}
    \right|}_{|\ul{\bf{H}}_2|}, \cdots,
    \ub{\left|
    \begin{matrix}
        f_{x_1x_1} & \cdots & f_{x_1x_n} \\
        \vdots & \ddots & \vdots \\
        f_{x_nx_1} & \cdots & f_{x_nx_n}
    \end{matrix}
    \right|}_{|\ul{\bf{H}}| = |\ul{\bf{H}}_n|} \]
\end{definition*}

\noindent
\textbf{Sylvestor's criterion} \\
\begin{center}
    $\ul{\bf{H}}$ is positive definite $\iff$
    signature is $+, +, \cdots , +$
\end{center}
\begin{center}
    $\ul{\bf{H}}$ is negative definite $\iff$
    signature is $-, +, -, \cdots, (-)^n$
\end{center}
[$\bf{x} \ul{\bf{H}} \bf{x}^\top$ with $\bf{x} =
(x_1, x_2, 0, 0, \dots, 0)$]

\subsection{Contours near stationary points}

Let $f(x, y)$ have a stationary point at $(x_0,
y_0)$. Adopt coordinates along principal axes of
$\ul{\bf{H}}$ at $\bf{x}_0$:
\[ \ul{\bf{H}} =
\begin{pmatrix}
    \lambda_1 & 0 \\
    0 & \lambda_2
\end{pmatrix}
\]
Write
\[ \bf{x} = \bf{x}_0 + (\xi, \zeta) \]
then around $\bf{x}_0$
\[ f(\bf{x}) \approx f(\bf{x}_0) + \half
(\lambda_1 \xi^2 + \lambda_2 \zeta^2) \]
Contours of $f$ close to $\bf{x}_0$ satisfy
\[ \lambda_1 \xi^2 + \lambda_2 \zeta^2 =
\text{constant} \]
Maximum or minimum: $\lambda_1$ and $\lambda_2$
same sign $\to$ elliptical contours \\
Saddle point: $\lambda_1$ and $\lambda_2$ opposite
sign $\to$ hyperbolic contours.

\begin{example*}
    \[ f(x, y) = 4x^3 - 12xy + y^2 + 10y - 6 \]
    \[ f_x = 12x^2 - 12y \]
    \[ f_y = -12x + 2y + 10 \]
    Stationary points:
    \[ f_x = 0 \implies x^2 = y \]
    \[ f_y = 0 \implies x^2 - 6x + 5 = 0 \implies
    x = 1, 5 \]
    Stationary points $(1, 1)$ and $(5, 25)$.
    \[ f_{xx} = 24x; f_{xy} = 12; f_{yy} = 2 \]
    \ul{$(1, 1)$}:
    \[ \ul{\bf{H}} =
    \begin{pmatrix}
        24 & -12 \\
        -12 & 2
    \end{pmatrix}
    \]
    \[ |\ul{\bf{H}}_1| = 24, |\ul{\bf{H}}_2| = -96 \]
    Signature is $+, -$ so \emph{neither} positive
    definite nor negative definite. $|\ul{\bf{H}}|
    < 0 \implies \lambda_1$ and $\lambda_2$ are
    non-zero and opposite signs $\to$
    \emph{saddle point}. \\
    \ul{$(5, 25)$}:
    \[ \ul{\bf{H}} =
    \begin{pmatrix}
        120 & -12 \\
        -12 & 2
    \end{pmatrix}
    \]
    $|\ul{\bf{H}}_1| = 120$, $|\ul{\bf{H}}_2| = 96$
    $\implies$ signature is $+, +$ so a
    \emph{minimum}. Near saddle point, $(x, y) =
    (1, 1) + (\delta x, \delta y)$. Contours
    locally have
    \[ f_{xx}(\delta x)^2 + 2 f_{xy}(\delta
    x)(\delta y) + f_{yy}(\delta y)^2 =
    \text{constant} \]
    \[ \implies 12(\delta x)^2 - 12 \delta x
    \delta y + (\delta y)^2 = \text{constant} \]
    Intersecting contours through $(1, 1)$ are
    asymptotes of hyperbola:
    \[ 12(\delta x)^2 - 12 \delta x \delta y +
    (\delta y)^2 = 0 \]
    \[ \implies \delta y = (6 \pm 2 \sqrt{6})
    \delta x .\]
\end{example*}

\newpage
\section{Systems of Linear Equations}

Consider two functions $y_1(t)$ and $y_2(t)$ with
\[ \dot{y}_1 = a y_1 + b y_2 + f_1(t) \]
\[ \dot{y}_2 = c y_1 + d y_2 + f_2(t) \]
($a, b, c, d$ constants). Vector form:
\[ \bf{\dot{Y}} = \ub{\ul{\bf{M}}}_{\ul{\bf{M}} =
\begin{pmatrix}
    a & b \\
    c & d
\end{pmatrix}
} \ub{\bf{Y}}_{\bf{Y} =
\begin{pmatrix}
    y_1 \\
    y_2
\end{pmatrix}
} + \ub{\bf{F}}_{
\begin{pmatrix}
    f_1 \\
    f_2
\end{pmatrix}
} \]
Two ways to solve.

\begin{enumerate}[(1)]
    \item Convert to higher-order ODE in one
        variable
        \begin{align*}
            \ddot{y}_1
            &= a \dot{y}_1 + b \dot{y}_2 +
            \dot{f}_1 \\
            &= a \dot{y}_1 + b(cy_1 + dy_2 + f_2)
            + \dot{f}_1 \\
            &= a \dot{y}_1 + bcy_1 + d(\dot{y}_1 -
            ay_2 - f_1) + bf_2 + \dot{f}_1
        \end{align*}
        \[ \implies \ddot{y}_1 - (a + d)\dot{y}_1
        + (ad - bc) y_1 = bf_2 - df_1 + \dot{f}_1 \]
        (Linear, second-order ODE with constant
        coefficients)
    \item Solve first-order ODEs directly with
        matrix methods. \\
        Sometimes convenient to deal with higher
        order ODEs as systems of first-order
        ODEs. For example
        \[ \ddot{y} + \alpha \dot{y} + \beta y = f \]
        Define:
        \[ y_1 \equiv y \qquad \text{and} \qquad
        y_2 \equiv \dot{y} \]
        so $\dot{y}_1 = y_2$ and $\dot{y}_2 =
        \ddot{y} = -\alpha y_2 - \beta y_1 + f$
        \[ \bf{\dot{Y}} =
        \begin{pmatrix}
            0 & 1 \\
            -\beta & -\alpha
        \end{pmatrix}
        \ub{\bf{Y}}_{
        \begin{pmatrix}
            y_1 \\
            y_2
        \end{pmatrix}
        } +
        \begin{pmatrix}
            0 \\
            f
        \end{pmatrix}
        \]
\end{enumerate}

\subsection{Matrix methods}

To solve
\[ \bf{\dot{Y}} = \ul{\bf{M}} \bf{Y} + \bf{F}(t) \]
\begin{enumerate}[(1)]
    \item Write
        \[ \bf{Y} = \bf{Y}_c + \bf{Y}_p \]
        \[ \bf{\dot{Y}}_c = \ul{\bf{M}} \bf{Y}_c \]
    \item Look for $\bf{Y}_c$ in the form
        \[ \bf{Y}_c = \bf{v} e^{\lambda t} \]
        Have
        \[ \bf{\dot{Y}}_c = \lambda \bf{Y}_c =
        \ul{\bf{M}} \bf{Y}_c \]
        \[ \implies \ul{\bf{M}} \bf{v} = \lambda
        \bf{v} \]
        ($\bf{v}$ must be an eigenvector of
        $\ul{\bf{M}}$ and $\lambda$ is the
        eigenvalue)
        For systems of $n$ equations with $n$ distinct
        eigenvalues $\lambda$, $n$ such complementary
        solutions which we can add to get general
        $\bf{Y}_c$.
    \item Find a particular solution $\bf{Y}_p$
        that satisfies the full system of ODEs.
\end{enumerate}

\begin{example*}
    \[ \bf{\dot{Y}} = \ub{
    \begin{pmatrix}
        -4 & 24 \\
        1 & -2
    \end{pmatrix}
    }_{\ul{\bf{M}}}
    \bf{Y} +
    \begin{pmatrix}
        4 \\
        1
    \end{pmatrix}
    e^t \]
    Complementary solutions:
    \[ \bf{Y}_c = \bf{v} e^{\lambda t} \]
    \[ \implies \ul{\bf{M}} \bf{v} = \lambda
    \bf{v} \implies |\ul{\bf{M}} - \lambda
    \ul{\bf{I}}| = 0 \]
    \[ \lambda_1 = 2 \qquad \text{or} \qquad
    \lambda_2 = -8 \]
    \[ \bf{v}_1 =
    \begin{pmatrix}
        4 \\
        1
    \end{pmatrix}
    \qquad \,\,\,\,\,\,\,\, \qquad \bf{v}_2 =
    \begin{pmatrix}
        -6 \\
        1
    \end{pmatrix}
    \]
    General complementary solution:
    \[ \bf{Y}_c = A
    \begin{pmatrix}
        4 \\
        1
    \end{pmatrix}
    e^{2t} + B
    \begin{pmatrix}
        -6 \\
        1
    \end{pmatrix}
    e^{-8t} \]
    For particular solution, try
    \[ \bf{Y}_p = \bf{u} e^t \]
    Requires
    \[ \bf{u} = \ul{\bf{M}} \bf{u} +
    \begin{pmatrix}
        4 \\
        1
    \end{pmatrix}
    \]
    \[ \implies \bf{u} = -(\ul{\bf{M}} -
    \ul{\bf{I}})^{-1}
    \begin{pmatrix}
        4 \\
        1
    \end{pmatrix}
    \]
    (inverse exists since 1 is not an eigenvalue
    of $\ul{\bf{M}}$)
    \[ \bf{u} = -
    \begin{pmatrix}
        4 \\
        1
    \end{pmatrix}
    \]
    General solution
    \[ \bf{Y} = A
    \begin{pmatrix}
        4 \\
        1
    \end{pmatrix}
    e^{2t} + B
    \begin{pmatrix}
        -6 \\
        1
    \end{pmatrix}
    e^{-8t} -
    \begin{pmatrix}
        4 \\
        1
    \end{pmatrix}
    e^t \]
    (Note: if forcing $\propto e^{\lambda t}$ with
    $\lambda$ an eigenvalue of $\ul{\bf{M}}$, try
    $\bf{Y}_p = \bf{u} t e^{\lambda t}$.)
\end{example*}

\subsection{Non-degenerate phase portraits}

$n$ first-order ODEs has phase space with points
\[ \bf{Y} = (y_1, y_2, \dots, y_n)^\top \]
\ul{Phase portrait}: solution trajectories in
phase space. \\
Consider
\[ \bf{\dot{Y}} = \ul{\bf{M}} \bf{Y} \]
(homogeneous). Fixed point $\bf{Y} = \bf{0}$. For
$n = 2$, and $\lambda_1 \neq \lambda_2$
(non-degenerate), general solution:
\[ \bf{Y}(t) = \bf{v}_1 e^{\lambda_1 t} + \bf{v}_2
e^{\lambda_2 t} \]
($\lambda_1, \lambda_2$ eigenvalues and $\bf{v}_1,
\bf{v}_2$ eigenvectors). Consider phase portraits
for $\lambda_1 \neq 0$, $\lambda_2 \neq 0$ (and
$\lambda_1 \neq \lambda_2$).

\bigskip

\noindent
\ul{Case 1}: $\lambda_1$ and $\lambda_2$ real and
opposite signs. WLOG $\lambda_2 < 0 < \lambda_1$.
$\bf{v}_1$ and $\bf{v}_2$ can be chosen to be
real.
\begin{center}
    \includegraphics[width=0.6\linewidth]
    {images/508022ce74a411ec.png}
\end{center}

\noindent
\ul{Case 2}: $\lambda_1$ and $\lambda_2$ real and
same sign. WLOG $|\lambda_1| > |\lambda_2|$ and
$\bf{v}_1$ and $\bf{v}_2$ real.
\begin{center}
    \includegraphics[width=0.6\linewidth]
    {images/b259f5d874a411ec.png}
\end{center}
For $\lambda_1, \lambda_2 < 0$ - sable node (all
arrows reversed). \\\\

\noindent
\ul{Case 3}: $\lambda_1$ and $\lambda_2$ complex
conjugate pairs, $\lambda_1 = \ol{\lambda_2}$,
$\bf{v}_1 = \ol{\bf{v}_2}$
\begin{align*}
    \bf{Y}(t)
    &= c \bf{v}_1 e^{\lambda_1 t} + c
    \cdot c_0 \\
    &= c \bf{v}_1 e^{\mathrm{Re}(\lambda_1)t} e^{i
    \mathrm{Im}(\lambda_1)t} + c \cdot c_0
\end{align*}
$\mathrm{Re}(\lambda_1) < 0 \implies$ stable
spiral.
\begin{center}
    \includegraphics[width=0.6\linewidth]
    {images/4cf2760674a511ec.png}
\end{center}
$\mathrm{Re}(\lambda_1) > 0 \implies$
unstable spiral.
\begin{center}
    \includegraphics[width=0.6\linewidth]
    {images/6dd9201874a511ec.png}
\end{center}

$\mathcal{Re}(\lambda_1) = 0 \to$ \emph{centre} \\
$\bf{Y}(t)$ is periodic $\to$ closed elliptical
trajectories in phase space
\begin{center}
    \includegraphics[width=0.6\linewidth]
    {images/fc72ee8274a811ec.png}
\end{center}
Sense of rotation? \\
Determine $\bf{\dot{Y}}$ at one point:
\begin{center}
    \includegraphics[width=0.6\linewidth]
    {images/31fe051e74a911ec.png}
\end{center}

\newpage
\section{Non-linear Dynamical Systems}

Consider an \emph{autonomous system} of two,
non-linear, first-order ODEs:
\[ \begin{cases}
    \dot{x} = f(x, y) \\
    \dot{y} = g(x, y)
\end{cases} \tag{$*$} \]
($f, g$ general non-linear functions of $x$, $y$
which are independent of $t$).

\subsection{Equilibrium points}

\begin{definition*}[Equilibrium point]
    An \emph{equilibrium point} (fixed point) of
    system ($*$) is a point at which $\dot{x} = 0$
    and $\dot{y} = 0$.
\end{definition*}

\noindent
If $(x_0, y_0)$ is a fixed point of ($*$):
\[ f(x_0, y_0) = 0 = g(x_0, y_0) \]
Solve these to locate the fixed points. Stability
via perturbation analysis:
\[ x(T) = x_0 + \xi(t); y(t) = y_0 + \zeta(t) \]
($\xi, \zeta$ small perturbations about $(x_0,
y_0)$)
\[ \dot{x} = f(x, y) \implies \dot{\xi} = f(x_0 +
\xi, y_0 + \zeta) \]
\[ \implies \dot{\xi} \approx \ub{f(x_0,
y_0)}_{=0\text{ since fixed point}} + \xi
\pfrac{f}{x} (x_0, y_0) + \zeta \pfrac{f}{y} (x_0,
y_0) \]
Similarly,
\[ \dot{\zeta} = \xi \pfrac{g}{x} (x_0, y_0) +
\zeta \pfrac{g}{y} (x_0, y_0) \]
\[ \implies
\begin{pmatrix}
    \dot{\xi} \\
    \dot{\zeta}
\end{pmatrix}
=
\ub{
\begin{pmatrix}
    f_x & f_y \\
    g_x & g_y
\end{pmatrix}
}_{\ul{\bf{M}}}
\begin{pmatrix}
    \xi \\
    \zeta
\end{pmatrix}
\]
Linear system of homogeneous, first-order ODEs
$\to$ eigenvalues of $\ul{\bf{M}}$ determine
stability of equilibrium points.

\bigskip

% \begin{example*}[predator-prey system]
\noindent\textbf{Example (predator-prey
system)}. \\
    $x(t)$: Number of prey at time $t$ \\
    $y(t)$: Number of predators at time $t$. \\
    Dynamics of prey:
    \[ \dot{x} = \alpha x - \beta x^2 - \gamma xy \]
    Where $\alpha x$ represents the excess births
    over natural deaths, $\beta x^2$ represents
    competition over scarce resource and $\gamma
    xy$ represents prey being killed - predators
    consume all prey they encounter. \\
    Dynamics of predators:
    \[ \dot{y} = \eps xy - \delta y \]
    where $\eps xy$ represents birth of predators
    and $\delta y$ represents natural death rate.
    \\
    Note $\alpha, \beta, \gamma, \eps, \delta$ are
    all positive constants. \\
    Specific case:
    \[ \dot{x} = \ub{8x - 2x^2 - 2xy}_{f(x, y)} \]
    \[ \dot{y} = \ub{xy - y}_{g(x, y)} \]
    Equilibrium points:
    \[ 2x(4 - x - y) = 0 \]
    \[ y(x - 1) = 0 \]
    First of these requires: $x = 0$ \emph{or} $x
    = 4 - y$. \\
    Second requires: $y = 0$ \emph{or} $y(3 - y) =
    0$ \\
    Have equilibrium points:
    \[ (0, 0), (1, 3), (4, 0) \]
    \[
    \begin{pmatrix}
        f_x & f_y \\
        g_x & g_y
    \end{pmatrix}
    =
    \begin{pmatrix}
        8 - 4x - 2y & -2x \\
        y & x - 1
    \end{pmatrix}
    \]
    \ul{$(0, 0)$}
    \[
    \begin{pmatrix}
        \dot{\xi} \\
        \dot{\zeta}
    \end{pmatrix}
    =
    \begin{pmatrix}
        8 & 0 \\
        0 & -1
    \end{pmatrix}
    \begin{pmatrix}
        \xi \\
        \zeta
    \end{pmatrix}
    \]
    Eigenvalues / eigenvectors:
    \[ 8,
    \begin{pmatrix}
        1 \\
        0
    \end{pmatrix}
    \qquad -1,
    \begin{pmatrix}
        0 \\
        1
    \end{pmatrix}
    \]
    \emph{Saddle node}
    \begin{center}
        \includegraphics[width=0.6\linewidth]
        {images/fb95366674ab11ec.png}
    \end{center}
    \ul{$(4, 0)$}
    \[
    \begin{pmatrix}
        \dot{\xi} \\
        \dot{\zeta}
    \end{pmatrix}
    =
    \begin{pmatrix}
        -8 & -8 \\
        0 & 3
    \end{pmatrix}
    \begin{pmatrix}
        \xi \\
        \zeta
    \end{pmatrix}
    \]
    Eigenvalues / eigenvectors:
    \[ -8,
    \begin{pmatrix}
        1 \\
        0
    \end{pmatrix}
    \qquad 3,
    \begin{pmatrix}
        8 \\
        -11
    \end{pmatrix}
    \]
    \emph{Saddle node}
    \begin{center}
        \includegraphics[width=0.6\linewidth]
        {images/2e94bd4874ac11ec.png}
    \end{center}
    \ul{$(1, 3)$}
    \[
    \begin{pmatrix}
        \dot{\xi} \\
        \dot{\zeta}
    \end{pmatrix}
    =
    \begin{pmatrix}
        -2 & -2 \\
        3 & 0
    \end{pmatrix}
    \begin{pmatrix}
        \xi \\
        \zeta
    \end{pmatrix}
    \]
    Eigenvalues: $-1 \pm i \sqrt{5}$ (note
    $\mathrm{Re}(\lambda) < 0$). \\
    \emph{Stable spiral}
    \begin{center}
        \includegraphics[width=0.6\linewidth]
        {images/b94cee4274ac11ec.png}
    \end{center}
    At $\begin{pmatrix}1 \\ 0\end{pmatrix}$,
    \[
    \begin{pmatrix}
        \dot{\xi} \\
        \dot{\zeta}
    \end{pmatrix}
    =
    \begin{pmatrix}
        -2 & -3 \\
        3 & 0
    \end{pmatrix}
    \begin{pmatrix}
        1 \\
        0
    \end{pmatrix}
    =
    \begin{pmatrix}
        -2 \\
        3
    \end{pmatrix}
    \]
    Note that the 3 in the last vector shows that
    $\dot{\zeta} > 0$ so the spiral must be
    anticlockwise.
    \begin{center}
        \includegraphics[width=0.6\linewidth]
        {images/56c8d0a074ad11ec.png}
    \end{center}
% \end{example*}

\newpage
\section{Partial Differential Equations}

PDEs - differential equations with multiple
independent variables.

\subsection{First-order wave equation}

Consider $\psi(x, t)$ satisfying
\[ \pfrac{\psi}{t} - c \pfrac{\psi}{x} = 0
\tag{$*$} \]
($c$ is a constant with dimensions of speed, $x$
is a 1 dimensional position and $t$ is a time)
Solve with \emph{method of characteristics}. How
does $\psi$ vary along path $x(t)$, so $\psi(x(t),
t)$?
\begin{align*}
    \dfrac{\psi}{t}
    &= \pfrac{\psi}{t} + \dfrac{x}{t}
    \pfrac{\psi}{x} \\
    &= \pfrac{\psi}{x} \left( c + \dfrac{x}{t} \right)
    &&\text{using ($*$)}
\end{align*}
If choose $x(t)$ such that
\[ \dfrac{x}{t} = -c \]
\[ \implies x = x_0 - ct \]
where $x_0$ is a constant along the path (label
paths). Then
\[ \dfrac{\psi}{x} = 0 \implies \text{$\psi =
\text{constant}$ along the path} \]
Paths $x = x_0 - ct$ are the characteristics of
($*$).
\begin{center}
\begin{tsqx}
    (-0.2,0)->>(1.5,0)
    label $x$ @ (1.6,0)
    label $t$ @ (0,1.1)
    (0,-0.2)->>(0,1.0)
    ! for (int i = 0; i < 5; ++i) {
        (0.5+0.25*i,0)--(0.25*i,1.0) red
    ! }
    (0.25,0)--(0,0.5) red
    (1.5,0.5)--(1.25,1.0) red
    label $x(t) = x_0 - ct$ @ (1.7,0.3) +col:blue
    label $x_0$ @ (1.25,-0.1) +col:blue
\end{tsqx}
\end{center}
As $\psi$ is constant along characteristics,
general solution of ($*$) is
\[ \psi(x, t) = f(x_0) = f(x + ct) \]
Simply translates the $x$-dependence of $\psi$ at
$t = 0$ to left by $ct$ at time $t$.
\begin{center}
    \includegraphics[width=0.6\linewidth]
    {images/7b935b16529911ec.png}
\end{center}
Left-moving wavelike solution.

\begin{example*}[unforced wave equation]
    \[ \pfrac{\psi}{t} - c \pfrac{\psi}{x} = 0 \]
    with $\psi(x, 0) = x^2 - 3$. General solution:
    $\psi = f(x + ct)$. At $t = 0$, $\psi(x, 0) =
    f(x) = x^2 - 3$ so
    \[ \psi(x, t) = (x + ct)^2 - 3 .\]
\end{example*}

\begin{example*}[forced wave equation]
    \[ \pfrac{\psi}{t} + 5 \pfrac{\psi}{x} = e^{-t} \]
    with $\psi(x, 0) = e^{-x^2}$. Characteristics
    have
    \[ \dfrac{x}{t} = 5 \implies x = x_0 + 5t \]
    Along these
    \[ \dfrac{\psi}{t} = e^{-t} \]
    \[ \implies \psi = f(x_0) - e^{-t} \]
    where $f$ is an arbitrary function that is
    constant on each characteristic. General
    solution:
    \[ \psi(x, t) = f(x - 5t) \]
    Initial conditions:
    \[ \psi(x,0) = f(x) - \ub{e^{-0}}_{1} = e^{-x^2} \]
    \[ \implies f(x) = 1 + e^{-x^2} \]
    Final solution:
    \[ \psi(x, t) = 1 + e^{-(x - 5t)^2} - e^{-t} \]
\end{example*}

\subsection{Second-order wave equation}

Allow propagation in both directions.
\[ \pfrac[2]{\psi}{t} - c^2 \pfrac[2]{\psi}{x} = 0 \]
Can be factored as
\[ \left( \pfrac{}{t} - c\pfrac{}{x} \right)
\left( \pfrac{}{t} + c \pfrac{}{x} \right) \psi =
0 \tag{\dag} \]
(works because partial derivatives commute). The
two first order operators commute so both
\[ \ub{f(x + ct)}_{\substack{\text{nulled
by}\\\pfrac{}{t} - c \pfrac{}{x}}} \qquad
\text{and} \qquad \ub{g(x - ct)}_{\substack{\text{nulled
by}\\\pfrac{}{t} + c \pfrac{}{x}}} \]
are solutions. General solution:
\[ \psi = f(x + ct) + g(x - ct) \]
for some arbitrary functions $f$ and $g$. \\
Most general solution? Yes! \\
Let $\zeta = x + ct$ and $\eta = x - ct$.
\[ \pfrac{}{x} = \ub{\left. \pfrac{\zeta}{x}
\right|_t}_{1} + \pfrac{}{\eta} \]
\[ \pfrac{}{t} = c \pfrac{}{\zeta} - c
\pfrac{}{\eta} \]
\[ \implies \pfrac{}{t} - c \pfrac{}{x} = - 2c
\pfrac{}{\eta} \]
and
\[ \pfrac{}{t} + c \pfrac{}{x} = 2x \pfrac{}{\zeta} \]
Wave equation (\dag) reduces to
\[ -4c^2 \frac{\partial^2 \psi}{\partial \eta
\partial \zeta} = 0 \]
\[ \implies \psi = f(\zeta) + g(\eta) \]
(Integrate twice)

\begin{example*}
    \[ \pfrac[2]{\psi}{t} - c^2 \pfrac[2]{\psi}{x}
    = 0 \]
    with
    \[ \psi(x, 0) = \frac{1}{1 + x^2} \qquad
    \text{and} \qquad \pfrac{\psi}{t}(x, 0) = 0 \]
    General solution:
    \[ \psi = f(x + ct) + g(x - ct) \]
    At $t = 0$:
    \[ \psi(x, 0) = f(x) + g(x) = \frac{1}{1 +
    x^2} \tag{$**$} \]
    \[ \pfrac{\psi}{t}(x, 0) = cf'(x) - cg'(x) = 0 \]
    \[ \implies f(x) = g(x) + A \]
    Combine with ($**$):
    \[ f(x) = \frac{1}{2(1 + x^2)} + \frac{A}{2} \]
    \[ g(x) = \frac{1}{2(1 + x^2)} - \frac{A}{2} \]
    Hence
    \[ \psi(x, t) = \half \left[ \frac{1}{1 + (x +
    ct)^2} + \frac{1}{1 + (x - ct)^2} \right] \]
    \begin{center}
    \begin{tsqx}
        (-1.5,0)->>(1.5,0)
        label $x$ @ (1.6,0)
        label $\psi$ @ (0,1.6)
        (0,0)->>(0,1.5)
        ! real f(real x) {return 0.5*2^(-x**2);}
        ! draw(graph(f,-1.5,1.5,operator ..), red);
        ! real g(real x) {return 0.5*2^(-(1 - x**2)**2);}
        ! draw(graph(g,-1.5,1.5,operator ..), blue);
        label $t = 0$ @ (0.3,1) +col:red
        label $t > 0$ @ (1.3,0.3) +col:blue
        label $\rightarrow$ @ (1,0.6) +col:blue
        label $\leftarrow$ @ (-1,0.6) +col:blue
    \end{tsqx}
    \end{center}
\end{example*}

\end{document}
