15  A1: Solutions

Jupyter Notebooks, LaTeX, Julia & Prerequisites

15.1 A: Jupyter Notebooks

Tip

You may do this part (part A only) with the help from your classmates

Exercise 1.

• Follow the steps here to install Julia and Jupyter,
  
• Upload and run this .ipynb file,
  
• Double click on the following cell to replace YOUR NAME HERE with your name

Name = "Name"

"Name"

15.2 B: \LaTeX

Exercise 2. Write down the degree n Taylor polynomial of f about some point a.

Answer.

Let the n^\text{th} degree Taylor polynomial of f about a be denoted T_nf(x). We have

\begin{align} T_nf(x) = \sum_{j=0}^n \frac{f^{(j)}(a)}{j!} (x - a)^j. \end{align}

Exercise 3. Apply the Taylor remainder theorem to f(x) = e^x about x=0 to show that 0 \leq 1+x \leq e^x for all x \geq -1

Answer.

By the Taylor remainder theorem, there exists \xi between 0 and x such that

\begin{align} f(x) = f(0) + f'(0) x + \frac{1}{2} f''(\xi) x^2. \end{align}

Since f'(x) = f''(x) = e^x, we have: for all x\geq -1,

\begin{align} 0 \leq 1 + x \leq 1 + x + \frac{1}{2}e^{\xi} x^2 = e^x. \end{align}

Here, we have used the fact that e^\xi x^2 \geq 0.

Exercise 4. Fix \alpha, \beta > 0 and suppose that x_n is a sequence such that x_{n+1} \leq (1 + \alpha) x_n + \beta for all n \in \mathbb N. Show that

\begin{align} x_{n+1} \leq (1 + \alpha)^{n+1} \left( x_0 + \frac{\beta}{\alpha} \right) - \frac{\beta}{\alpha} \end{align}

Answer.

We will show by induction that

\begin{align} x_{n+1} &\leq (1 + \alpha)^{n+1} x_0 + \beta \sum_{j=0}^{n} (1 + \alpha)^j. \end{align}

We are given that this is true for n = 0. Supposing this statement is true for n, we have

\begin{align} x_{n+2} &\leq (1 + \alpha) x_{n+1} + \beta \nonumber\\ % &\leq (1 + \alpha)\left[ (1 + \alpha)^{n+1} x_0 + \beta \sum_{j=0}^{n} (1 + \alpha)^j \right] + \beta \nonumber\\ % &= (1 + \alpha)^{n+2} x_{0} + \beta \left[ 1 + (1+\alpha)\sum_{j=0}^{n} (1 + \alpha)^j \right] \nonumber\\ % &= (1 + \alpha)^{n+2} x_{0} + \beta \sum_{j=0}^{n+1} (1 + \alpha)^j % \end{align}

as required.

We conclude by applying the geometric series formula: \sum_{j=0}^n(1 + \alpha)^j = [(1 + \alpha)^{n+1} - 1]/ \alpha.

Answer 2. Some people defined y_n := x_{n} + \frac{\beta}{\alpha} which satisfies y_{n+1} = (1 + \alpha) y_n = (1+\alpha)^{n+1} y_0. Substituting x_{n+1} and x_0 back into this equation gives the desired result.

(This is a much nicer solution than mine)

15.3 C. Julia

Suppose u : [0,T] \to \mathbb R solves the initial value problem

\begin{align} u'(t) &= f\big( t, u(t) \big) \qquad \text{on } (0,T) \nonumber \\ u(0) &= u_0. \tag{IVP} \end{align}

Recall from L1 that, in order to approximate the solution u at mesh points t_j = j h where h := \frac{T}{n} (for j=0,\dots,n), we may approximate the integral

\begin{align} u(t_{j+1}) - u( t_j ) = \int_{t_j}^{t_{j+1}} f\big( s, u(s) \big) \mathrm{d}s. \end{align}

In Math 5485, you used the Rectangular rule: u(t_j) \approx u_j where u_0 is given as the initial condition and

\begin{align} u_{j+1} := u_j + h f( t_j, u_j ). \end{align}

This is Euler’s method:

function Euler( u0, f, T, n )
    h = T/n     
    t = 0:h:T    
    u = zeros(n+1)
    u[1] = u0
    for j = 1:n
        u[j+1] = u[j] + h * f( t[j], u[j] )
    end
    return u
end 

Euler (generic function with 1 method)

Exercise 5. Use Euler’s method to approximate the solution of \text{(IVP)} for

\begin{align} &u_0 = 1, \quad T = 1, \quad \text{and} \\ &f(t, u) := \frac{2(1 - t u)}{1 + t^2}. \end{align}

In this case the exact solution is known and given by u(t) = \frac{1 + 2t}{1 + t^2}. Use the exact solution to compute the approximation errors \max_j | u(t_j) - u_j | (where u is the exact solution and u_j is the approximation at t_j). At what rate does the approximation error decay as a function of n?

Hints:

• some of the code has been written for you,
• mesh is the array ( t_j )_{j=0}^n,
• u_exact is the exact solution evaluated on the mesh,
• You may use maximum( @. abs( u - u_exact ) ) to compute the error between the exact solution u_exact and the approximation u.

using Plots, LaTeXStrings
u0, T, n = 1, 1, 100
mesh = 0:T/n:T
f(t, u) = 2*(1-t*u)/(t^2 + 1)
u_exact = t -> (2*t+1)/(t^2 + 1)

u = Euler(u0, f, T, n)

plot( t -> (2*t+1)/(t^2 + 1), 0, 1, 
            label="Exact solution" )
plot!( mesh, u, 
            label="Approximation (Euler, n = $n)" )

N = 500
errs = zeros(N)
for n ∈ 1:N
    u = Euler( u0, f, T, n )
    mesh = 0:T/n:T
    errs[n] = maximum( @. abs( u - u_exact(mesh) ) )
end

scatter( errs, xaxis=:log, yaxis=:log, xlabel=L"n", label="error", lw=3, title="Euler's Method" )
plot!( (1:N).^(-1), label=L"n^{-1}", linestyle=:dash )

We can see from the above plot that the error decays with rate O(n^{-1}).