Homework 7


Directions: Add work to this notebook to solve the problems below.

Be careful to name the objects you create as described in the problem. Some problems may include some tests you can run to test your code.

This assignment will not be collected.

Some problems are derived from problems in these books:

  • LL: Programming for Computations - Python by Svein Linge and Hans Petter Langtangen, 2nd edition.
  • L: A Primer on Scientific Programming with Python by Hans Petter Langtangen, 2nd edition.
  • TAK: Applied Scientific Computing With Python by Peter R. Turner, Thomas Arildsen, and Kathleen Kavanagh.

Standard imports

In [ ]:
import numpy as np

1. Solving a Bessel equation

(Before starting this problem read Example 7, which starts at the bottom of page 94 of TAK. You may want to go back and read Example 1 on page 36.)

Consider solving the Bessel equation $$x^2 y'' + x y' + (x^2 - 1) y = 0$$ subject to the boundary conditions $y(1)=1$ and $y(15)=0$ using numerical derivatives with $N = 280$ steps between $[1, 15]$. You should begin by rewriting the differential equation in the form $$y'' + \frac{1}{x} y' + (1-\frac{1}{x^2}) y = 0.$$ Next, use finite difference approximations to the derivatives and collect like terms to obtain a tridiagonal linear system for the unknown discrete values of y.

Use the function tridiagonal_solve from the previous assignment to compute the solution $y$. Graph the solution.

Some help: The interval is divided into $280$ steps by the points in x below.

In [ ]:
N = 280
x = np.linspace(1,15,N+1)

Step size is then:

In [ ]:
h = x[1] - x[0]
print(f'h = {h:.2f}')

We are looking to find a numpy array $y$ which can be obtained by multiplying by a tridiagonal matrix to get some vector we know. Then we can use tridiagonal_solve to find $y$. Here the entries of $y$, $y_i$, should be given by $y_i=y(x_i)$, where $x_i$ is the numpy array above. Thus $y(x_i+h)=y_{i+1}$ and $y(x_i-h)=y_{i-1}$ for $i=1, \ldots, N-1$.

Observe that $$\left[\begin{array}{rrrrrr} a_0 & b_0 & & & & \\ c_0 & a_1 & b_1 & & & \\ & c_1 & a_2 & b_2 & & \\ & & \ddots & \ddots & \ddots & \\ & & & c_{N-2} & a_{N-1} & b_{N-1} \\ & & & & c_{N-1} & a_{N} \end{array}\right] \left(\begin{array}{r} y_0 \\ y_1 \\ y_2 \\ \vdots \\ y_{N-1} \\ y_{N}\end{array}\right)= \left(\begin{array}{r} a_0 y_0 + b_0 y_1 \\ c_0 y_0 + a_1 y_1 + b_1 y_2 \\ c_1 y_1 + a_2 y_2 + b_2 y_3 \\ \vdots \\ c_{N-2} y_{N-2} + a_{N-1} y_{N-1} + b_{N-1} y_{N} \\ c_{N-1} y_{N-1} + a_{N} y_{N} \end{array}\right)=\left(\begin{array}{r} r_0 \\ r_1 \\ r_2 \\ \vdots \\ r_{N-1} \\ r_{N}\end{array}\right). $$ You can find the values of $a$, $b$, $c$ and $r$ by analyzing the differential equation and the boundary conditions. (Probably you want to make the boundary conditions correspond to the first and last equation.)

In [ ]:

2. Another differential equation

Solve $y'' + 2x y' + 2y = 3e^{−x} − 2xe^{−x}$ subject to $y(−1) = e + 3/e$, $y (1) = 4/e$ using the second-order finite difference method as in the above problem. Consider $N = 10, 20, 40$ in $[−1, 1]$. Plot these solutions and the true solution $y = e^{−x} + 3e^{−x^2}$. Also, on a separate set of axes, plot the three error curves. Estimate the order of the truncation error for this method. Is it what you would expect?

In [ ]:

3. A third differential equation

Consider the differential equation $$y''(x) + \sin(x) y(x) = x^2$$ on the interval $[0, \pi]$ subject to the boundary conditions $y(0)=1$ and $y(\pi)=0$.

Given an integer $n \geq 3$, consider cutting the interval into $n$ equal-sized subintervals. You are interested in approximating the solution $y$ at the points $$x_0 = 0, x_1 = \frac{\pi}{n}, \ldots, x_{n-1}=\frac{(n-1)\pi}{n}, x_n=\pi.$$ Set $y_i = y(x_i)$. Using numerical derivatives, find an $n+1 \times n+1$ matrix $M$ and a vector $v \in {\mathbb R}^{n+1}$ such that solving $M y= v$ leads to an approximate solution $y=(y_0, \ldots, y_{n-1})$ for the differential equation.

Write a function linearization(n) which returns a pair of numpy arrays (m,v). Here m should represent the matrix $M$ from the previous paragraph and v should reprensent the vector $v$.

Write a function solution(n) which returns the $y$ that solves the equation $My=v$.

In [ ]: