# Worksheet 4


```
%matplotlib inline
```

## Question 1

Convert the ODE

\begin{equation}
  y''' + x y'' + 3 y' + y = e^{−x}
\end{equation}

into a first order system of ODEs.

### Answer Question 1

Step by step we introduce

\begin{align}
  u &= y' \\
  v &= u' \\
    &= y''.
\end{align}

We can therefore write the ODE into a system of ODEs. The first order ODEs for $y$ and $u$ are given by the definitions above. The ODE for $v$ is given from the original equation, substituting in the definition of $u$ where appropriate, to get

\begin{align}
  \begin{pmatrix} y \\ u \\ v \end{pmatrix}' & = \begin{pmatrix} u \\ v \\ e^{-x} - x y'' - 3 y' - y \end{pmatrix} \\
  & = \begin{pmatrix} u \\ v \\ e^{-x} - x v - 3 u - y \end{pmatrix}.
\end{align}

## Question 2

Show by Taylor expansion that the backwards differencing estimate of $f(x)$,

\begin{equation}
  f(x) = \frac{f(x) − f(x − h)}{h}
\end{equation}

is first order accurate.

### Answer Question 2

We have the Taylor series expansion of $f(x − h)$ about $x$ is

\begin{equation}
  f(x − h) = f(x) − h f'(x) + \frac{h^2}{2!} f''(x) + {\mathcal O} (h^3).
\end{equation}

Substituting this in to the backwards difference formula we find

\begin{align}
  \frac{f(x) - f(x - h)}{h} & = \frac{f(x) - f(x) + h f'(x) + \frac{h^2}{2!} f''(x) + {\mathcal O} (h^3)}{h} \\
  & = f'(x) + {\mathcal O} (h)
\end{align}

Therefore the difference between the exact derivative $f'$ and the backwards difference estimate is $\propto h$ and hence the finite difference estimate is first order accurate.

## Question 3

Use Taylor expansion to derive a symmetric or central difference estimate of $f^{(4)}(x)$ on a grid with spacing $h$.

### Answer Question 3

For this we need the Taylor expansions

\begin{align*}
    f(x + h) & = f(x) + h f^{(1)}(x) + \frac{h^2}{2!} f^{(2)}(x) +
    \frac{h^3}{3!} f^{(3)}(x) + \frac{h^4}{4!} f^{(4)}(x) +
    \frac{h^5}{5!} f^{(5)}(x) + \dots \\
    f(x - h) & = f(x) - h f^{(1)}(x) + \frac{h^2}{2!} f^{(2)}(x) -
    \frac{h^3}{3!} f^{(3)}(x) + \frac{h^4}{4!} f^{(4)}(x) -
    \frac{h^5}{5!} f^{(5)}(x) + \dots \\
    f(x + 2 h) & = f(x) + 2 h f^{(1)}(x) + \frac{4 h^2}{2!} f^{(2)}(x) +
    \frac{8 h^3}{3!} f^{(3)}(x) + \frac{16 h^4}{4!} f^{(4)}(x) +
    \frac{32 h^5}{5!} f^{(5)}(x) + \dots \\
    f(x - 2 h) & = f(x) - 2 h f^{(1)}(x) + \frac{4 h^2}{2!} f^{(2)}(x) -
    \frac{8 h^3}{3!} f^{(3)}(x) + \frac{16 h^4}{4!} f^{(4)}(x) -
    \frac{32 h^5}{5!} f^{(5)}(x) + \dots 
\end{align*}

By a central or symmetric difference estimate we mean that the coefficient of $f(x \pm n h)$ should have the same magnitude. By comparison with central difference estimates for first and second derivatives we see that for odd order derivatives the coefficients should have opposite signs and for even order the same sign.

So we write our estimate as

\begin{equation*}
    f^{(4)}(x) \simeq A f(x) + B \left( f(x + h) + f(x - h) \right)
    + C \left( f(x + 2 h) +  f(x - 2 h) \right)
\end{equation*}

and we then need to constrain the coefficients $A, B, C$. By looking at terms proportional to $h^s$ we see

\begin{align*}
    h^0: && 0 & = A + 2 B + 2 C \\
    h^1: && 0 & = 0 \\
    h^2: && 0 & = B + 4 C \\
    h^3: && 0 & = 0 \\
    h^4: && \frac{1}{h^4} & = \frac{B}{12} + \frac{16 C}{12}. 
\end{align*}

This gives three constraints on our three unknowns so we cannot go to higher order. Solving the equations gives

\begin{equation*}
    A = \frac{6}{h^4}, \qquad B = -\frac{4}{h^4}, \qquad C =
    \frac{1}{h^4}.
\end{equation*}

Writing it out in obvious notation we have

\begin{equation*}
    f_1^{(4)} = \frac{1}{h^4} \left( 6 f_i - 4 (f_{i+1} + f_{i-1}) +
      (f_{i+2} + f_{i-2}) \right).
\end{equation*}

## Question 4

State the convergence rate of Euler's method and the Euler predictor-corrector method.

### Answer Question 4

Euler's method converges as $h$ and the predictor-corrector method as $h^2$.

## Question 5

Explain when multistage methods such as Runge-Kutta methods are useful.

### Answer Question 5

Multistage methods require only one vector of initial data, which must be provided to completely specify the IVP; that is, the method is self-starting. It is also easy to adapt a multistage method to use variable step sizes; that is, to make the algorithm adaptive depending on local error estimates in order to keep the global error within some tolerance. Finally, it is relatively easy to theoretically show convergence. Combining this we see that multistage methods are useful as generic workhorse algorithms and in cases where the function defining the IVP may vary widely in behaviour, so that adaptive algorithms are required.

## Question 6

Explain the power method for finding the largest eigenvalue of a matrix. In particular, explain why it is simpler to find the absolute value, and how to find the phase information.

### Answer Question 6

The idea behind the power method is that most easily seen by writing out a generic vector ${\bf x}$ in terms of the eigenvectors of the matrix $A$ whose eigenvalues we wish to find,

\begin{equation}
  {\bf x} = \sum_{i=1}^N a_i {\bf e}_i,
\end{equation}

where we assume that the eigenvectors are ordered such that the associated eigenvalues have the order $|\lambda_1 | > |\lambda_2 | \ge |\lambda_3 | \ge \dots \ge |\lambda_N |$. Note that we always assume that there is a unique eigenvalue $\lambda_1$ with largest magnitude.

We then note that multiplying this generic vector by the matrix $A$ a number of times gives

\begin{equation}
  A^k {\bf x} = \lambda_1^k \sum_{i=1}^N a_i \left( \frac{\lambda_i}{\lambda_1} \right)^k {\bf e}_i.
\end{equation}

We then note that, for $i = 1$, the ratio of the eigenvalues $(\lambda_i / \lambda_1)^k$ must tend to zero as $k \to \infty$. Therefore in the limit we will "pick out" $\lambda_1$.

Of course, to actually get the eigenvalue itself we have to essentially divide two vectors. That is, we define a sequence $x^{(k)}$ where the initial value $x^{(0)}$ is arbitrary and at each step we multiply by $A$, so that

\begin{equation}
  x^{(k)} = A^k x^{(0)}.
\end{equation}

It follows that we can straightforwardly get $\lambda_1$ by looking at "the ratio of successive iterations". E.g.,

\begin{equation}
  \lim_{k \to \infty} \frac{ \| {\bf x}_{k+1} \| }{ \| {\bf x}_k \| } = | \lambda_1 |.
\end{equation}

This only gives information about the magnitude as we have used the simplest way of getting from a vector to a real number, the absolute value. To retain information about the phase we need to replace the absolute value of the vectors with some linear functional such as the sum of the coefficients.

## Coding Question 1

Apply Euler's method to the ODE

\begin{equation}
  y' + 2y = 2 − e^{−4 x}, \qquad y(0) = 1.
\end{equation}

Find the value of $y(1)$ (analytic answer is $1 − (e^{−2} − e^{−4})/2$) and see how your method converges with resolution.

### Answer Coding Question 1


```
def Euler(f, y0, interval, N = 100):
    """Solve the IVP y' = f(x, y) on the given interval using N+1 points (counting the initial point) with initial data y0."""
    import numpy as np
    
    h = (interval[1] - interval[0]) / N
    x = np.linspace(interval[0], interval[1], N+1)
    y = np.zeros((len(y0), N+1))
    y[:, 0] = y0
    
    for i in range(N):
        y[:, i+1] = y[:, i] + h * f(x[i], y[:, i])
    
    return x, y

def fn_q1(x, y):
    """Function defining the IVP in question 1."""
    import numpy as np
    
    return 2.0 - np.exp(-4.0*x) - 2.0*y

# Now do the test
import numpy as np

exact_y_end = 1.0 - (np.exp(-2.0) - np.exp(-4.0)) / 2.0

# Test at default resolution
x, y = Euler(fn_q1, np.array([1.0]), [0.0, 1.0])
print "Error at the end point is ", y[:, -1] - exact_y_end

import matplotlib.pyplot as plt
fig = plt.figure(figsize = (12, 8), dpi = 50)
plt.plot(x, y[0, :], 'b-+')
plt.xlabel('$x$', size = 16)
plt.ylabel('$y$', size = 16)

# Now do the convergence test

levels = np.array(range(4, 10))
Npoints = 2**levels
abs_err = np.zeros(len(Npoints))
for i in range(len(Npoints)):
    x, y = Euler(fn_q1, np.array([1.0]), [0.0, 1.0], Npoints[i])
    abs_err[i] = abs(y[0, -1] - exact_y_end)

# Best fit to the errors
h = 1.0 / Npoints
p = np.polyfit(np.log(h), np.log(abs_err), 1)
    
fig = plt.figure(figsize = (12, 8), dpi = 50)
plt.loglog(h, abs_err, 'kx')
plt.loglog(h, np.exp(p[1]) * h**(p[0]), 'b-')
plt.xlabel('$h$', size = 16)
plt.ylabel('$|$Error$|$', size = 16)
plt.legend(('Euler Errors', "Best fit line slope {:.3}".format(p[0])), loc = "upper left")

plt.show()
```

## Coding Question 2

Apply the standard RK4 method to the above system, again checking that it converges with resolution.

### Answer Coding Question 2


```
def RK4(f, y0, interval, N = 100):
    """Solve the IVP y' = f(x, y) on the given interval using N+1 points (counting the initial point) with initial data y0."""
    import numpy as np
    
    h = (interval[1] - interval[0]) / N
    x = np.linspace(interval[0], interval[1], N+1)
    y = np.zeros((len(y0), N+1))
    y[:, 0] = y0
    
    for i in range(N):
        k1 = h * f(x[i]          , y[:, i])
        k2 = h * f(x[i] + h / 2.0, y[:, i] + k1 / 2.0)
        k3 = h * f(x[i] + h / 2.0, y[:, i] + k2 / 2.0)
        k4 = h * f(x[i] + h      , y[:, i] + k3)
        y[:, i+1] = y[:, i] + (k1 + k4 + 2.0 * (k2 + k3)) / 6.0
    
    return x, y

def fn_q2(x, y):
    """Function defining the IVP in question 2."""
    import numpy as np
    
    return 2.0 - np.exp(-4.0*x) - 2.0*y

# Now do the test
import numpy as np

exact_y_end = 1.0 - (np.exp(-2.0) - np.exp(-4.0)) / 2.0

# Test at default resolution
x, y = RK4(fn_q1, np.array([1.0]), [0.0, 1.0])
print "Error at the end point is ", y[:, -1] - exact_y_end

import matplotlib.pyplot as plt
fig = plt.figure(figsize = (12, 8), dpi = 50)
plt.plot(x, y[0, :], 'b-+')
plt.xlabel('$x$', size = 16)
plt.ylabel('$y$', size = 16)

# Now do the convergence test

levels = np.array(range(4, 10))
Npoints = 2**levels
abs_err = np.zeros(len(Npoints))
for i in range(len(Npoints)):
    x, y = RK4(fn_q1, np.array([1.0]), [0.0, 1.0], Npoints[i])
    abs_err[i] = abs(y[0, -1] - exact_y_end)

# Best fit to the errors
h = 1.0 / Npoints
p = np.polyfit(np.log(h), np.log(abs_err), 1)

fig = plt.figure(figsize = (12, 8), dpi = 50)
plt.loglog(h, abs_err, 'kx')
plt.loglog(h, np.exp(p[1]) * h**(p[0]), 'b-')
plt.xlabel('$h$', size = 16)
plt.ylabel('$|$Error$|$', size = 16)
plt.legend(('RK4 Errors', "Best fit line slope {0:.3}".format(p[0])), loc = "upper left")

plt.show()
```

## Coding Question 3

Write a code using the power method and inverse power method to compute the largest and smallest eigenvalues of an arbitrary matrix. Apply it to a random $n = 3$ matrix, checking that the correct answer is found. How does the number of iterations required for convergence to a given level vary with the size of the matrix?

### Answer Coding Question 3


```
def PowerMethod(A, tolerance = 1e-10, MaxSteps = 100):
    """Apply the power method to the matrix A to find the largest eigenvalue in magnitude."""
    import numpy as np
    import numpy.linalg as la
    
    n = np.size(A, 0)
    # Simple initial value
    x = np.ones(n)
    x /= la.norm(x)
    ratio = 1.0
    for k in range(MaxSteps):
        ratio_old = ratio
        x_old = x.copy()
        x = np.dot(A, x)
        ratio = np.sum(x) / np.sum(x_old)
        x /= la.norm(x)
        if (abs(ratio - ratio_old) < tolerance):
            break
    
    return [ratio, k]

def InversePowerMethod(A, tolerance = 1e-10, MaxSteps = 100):
    """Apply the inverse power method to the matrix A to find the smallest eigenvalue in magnitude."""
    import numpy as np
    import scipy.linalg as la
    
    n = np.size(A, 0)
    # Simple initial value
    x = np.ones(n)
    x /= la.norm(x)
    ratio = 1.0
    for k in range(MaxSteps):
        ratio_old = ratio
        x_old = x.copy()
        x = la.solve(A, x)
        ratio = np.sum(x) / np.sum(x_old)
        x /= la.norm(x)
        if (abs(ratio - ratio_old) < tolerance):
            break
    
    return [1.0/ratio, k]

# Test on a random 3x3 matrix
import numpy as np
import scipy.linalg as la

A = np.random.rand(3,3)
max_lambda, iterations_max = PowerMethod(A)
min_lambda, iterations_min = InversePowerMethod(A)
eigenvalues, eigenvectors = la.eig(A)
print "Computed maximum and minimum eigenvalues are", max_lambda, min_lambda
print "True eigenvalues are", eigenvalues

# Now check how the number of iterations depends on the matrix size.
# As we are computing random matrices, do average of 10 attempts
MinMatrixSize = 3
MaxMatrixSize = 50
Attempts = 10
iterations = np.zeros((MaxMatrixSize - MinMatrixSize + 1, Attempts))
for n in range(MinMatrixSize, MaxMatrixSize):
    for a in range(Attempts):
        A = np.random.rand(n, n)
        ratio, iterations[n - MinMatrixSize, a] = PowerMethod(A)

import matplotlib.pyplot as plt
ii = np.mean(iterations, 1)
nn = np.array(range(MinMatrixSize, MaxMatrixSize))
fig = plt.figure(figsize = (12, 8), dpi = 50)
plt.plot(range(MinMatrixSize, MaxMatrixSize+1), np.mean(iterations, 1), 'kx')
plt.xlabel('Matrix Size')
plt.ylabel('Mean number of iterations')

plt.show()
```

We see that the number of iterations is practically unchanged with the size of the matrix.


```
from IPython.core.display import HTML
def css_styling():
    styles = open("../../IPythonNotebookStyles/custom.css", "r").read()
    return HTML(styles)
css_styling()
```




<style>
    @font-face {
        font-family: "Computer Modern";
        src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');
    }
    div.cell{
        width:800px;
        margin-left:16% !important;
        margin-right:auto;
    }
    h1 {
        font-family: Verdana, Arial, Helvetica, sans-serif;
    }
    h2 {
        font-family: Verdana, Arial, Helvetica, sans-serif;
    }
    h3 {
        font-family: Verdana, Arial, Helvetica, sans-serif;
    }
    div.text_cell_render{
        font-family: Gill, Verdana, Arial, Helvetica, sans-serif;
        line-height: 110%;
        font-size: 120%;
        width:700px;
        margin-left:auto;
        margin-right:auto;
    }
    .CodeMirror{
            font-family: "Source Code Pro", source-code-pro,Consolas, monospace;
    }
/*    .prompt{
        display: None;
    }*/
    .text_cell_render h5 {
        font-weight: 300;
        font-size: 12pt;
        color: #4057A1;
        font-style: italic;
        margin-bottom: .5em;
        margin-top: 0.5em;
        display: block;
    }

    .warning{
        color: rgb( 240, 20, 20 )
        }  
</style>




> (The cell above executes the style for this notebook. It closely follows the style used in the [12 Steps to Navier Stokes](http://lorenabarba.com/blog/cfd-python-12-steps-to-navier-stokes/) course.)

0497f055451dfd2ed8228bda697a6b9a00dac896
ipynb
Jupyter Notebook
Worksheets/Worksheet4_Notebook.ipynb
srinivasvl81/NumericalMethods
7bc25d882f287a21e9646c676c95672414c980ff
CC-BY-3.0
2021-12-01T09:15:04.000Z
2021-12-01T09:15:04.000Z
Worksheets/Worksheet4_Notebook.ipynb
indranilsinharoy/NumericalMethods
989e0205565131057c9807ed9d55b6c1a5a38d42
MIT
Worksheets/Worksheet4_Notebook.ipynb
indranilsinharoy/NumericalMethods
989e0205565131057c9807ed9d55b6c1a5a38d42
MIT
2019-11-08T06:52:24.000Z
2019-11-08T06:52:24.000Z
Qwen/Qwen-72B
1. YES
2. YES