Problems

Optimization

Newton's Method

Write a function called newton which takes input parameters $f$, $x_0$, $h$ (with default value 0.001), tolerance (with default value 0.001) and max_iter (with default value 100). The function implements Newton's method to approximate a solution of $f(x) = 0$. In other words, compute the values of the recursive sequence starting at $x_0$ and defined by

$$ x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)} $$

Use the central difference formula with step size $h$ to approximate the derivative $f'(x_n)$. The desired result is that the method converges to an approximate root of $f(x)$ however there are several possibilities:

The sequence reaches the desired tolerance $|f(x_n)| \leq \mathtt{tolerance}$ and newton returns the value $x_n$.
The number of iterations exceeds the maximum number of iterations max_iter, the function prints the statement "Maximum iterations exceeded" and returns None.
A zero derivative is computed $f'(x_n) = 0$, the function prints the statement "Zero derivative" and returns None.

Numerical Differentiation

Central Difference Formula

Write a function called derivatives which takes input parameters $f$, $a$, $n$ and $h$ (with default value h = 0.001) and returns approximations of the derivatives $f'(a)$, $f''(a)$, $\dots$, $f^{(n)}(a)$ (as a NumPy array) using the formula

$$ f^{(n)}(a) \approx \frac{1}{2^n h^n} \sum_{k=0}^n (-1)^k {n \choose k} f \left( a + ( n - 2k ) h \right) $$

Use either scipy.misc.factorial or scipy.misc.comb to compute $n$ choose $k$:

$$ {n \choose k} = \frac{n!}{k!(n-k)!} $$

Taylor Polynomials

Write a function called taylor which takes input parameters $f$, $a$, $n$ and $L$ and plots both $f(x)$ and the Taylor polynomial $T_n(x)$ of $f(x)$ at $x=a$ of degree $n$

$$ T_n(x) = \sum_{k=0}^n \frac{f^{(k)}(a)}{k!}(x - a)^k $$

on the interval $[a-L,a+L]$ (in the same figure).

Numerical Integration

Trapezoid Rule

Find $f''(x)$ for $f(x) = \ln( \ln x)$.
Prove $|f''(x)| \leq \displaystyle \frac{2}{e^2}$ for $x \geq e$.
Write a function called log_log which takes input parameters u and abs_tolerance such that $u \geq e$ and abs_tolerance is a positive number (with default value 0.0001). The function uses the trapezoid rule to compute and return an approximation of the integral $$ \int_e^u \ln( \ln x) dx $$ The number $N$ of subintervals used in the trapezoid rule must be large enough to guarantee that the approximation is within abs_tolerance of the true value. You may use the function scipy.integrate.trapz.

Simpson's Rule

Find $f''''(x)$ for $f(x) = e^{-x^2}$.
Plot $f''''(x)$ for $x \in [0,5]$. Determine a bound $M$ such that $|f''''(x)| \leq M$ for $x \geq 0$.
Write a function called erf which takes input parameters u and abs_tolerance (with default value 0.0001) such that both are positive numbers. The function uses Simpson's rule to compute and return an approximation of the integral $$ \int_0^u e^{-x^2} dx $$ The number $N$ of subintervals used in Simpson's rule must be large enough to guarantee that the approximation is within abs_tolerance of the true value. You may use the function scipy.integrate.simps.

Differential Equations

Lorenz Equations

The Lorenz equations are the system of nonlinear differential equations

\begin{align} \frac{dx}{dt} &= \sigma(y - x) \\ \frac{dy}{dt} &= x(\rho - z) - y \\ \frac{dz}{dt} &= xy - \beta z \end{align}

where $\sigma$, $\rho$ and $\beta$ are positive numbers. Write a function called lorenz which takes input parameters sigma, rho, beta, u0, t0, tf, N and plot_vars (with default value [0,1]). The function computes and plots a numerical approximation of the corresponding solution of the Lorenz equations using the function scipy.integrate.odeint. The input parameters are:

sigma, rho and beta define the parameters $\sigma$, $\rho$ and $\beta$
u0 is a list of numbers of length 3 defining the initial conditions $[x(t_0),y(t_0),z(t_0)]$
t0 is the start of the interval of integration $[t_0,t_f]$
tf is the end of the interval of integration $[t_0,t_f]$
N is an integer specifying the number of evenly spaced points from $t_0$ to $t_f$ (inclusively) over which to compute the solution of the system
plot_vars is a list of length 2 specifying which 2 components to plot where $x=0$, $y=1$ and $z=2$. For example, if plot_vars is $[0,1]$ then plot the solution $x$ versus $y$. If plot_vars is $[1,2]$ then plot the solution $y$ versus $z$. Note $x$ versus $y$ means $x$ is the horizontal axis and $y$ is the vertical. Default value is [0,1] which plots $x$ versus $y$.

The function lorenz returns a 2D NumPy array with 4 columns:

column at index 0 is the array of $N$ evenly spaced $t$ values from $t_0$ to $t_f$ (inclusively)
column at index 1 is the array of $x$ values of the solution
column at index 2 is the array of $y$ values of the solution
column at index 3 is the array of $z$ values of the solution

Damped Oscillator

Write a function called damping which takes input parameters m, b, k, F, u0, t0, tf and N. The function uses scipy.integrate.odeint to compute a numerical approximation of the corresponding solution of the nonlinear damping equation:

$$ m y'' + b |y'| y' + ky = F(t) $$

The input parameters are:

m, b and k are positive numbers in the nonlinear damping equation
F is a function of one variable $F(t)$ in the nonlinear damping equation
u0 is a list of numbers of length 2 defining the initial conditions $[y(t_0),y'(t_0)]$
t0 is the start of the interval of integration $[t_0,t_f]$
tf is the end of the interval of integration $[t_0,t_f]$
N is an integer specifying the number of evenly spaced points from $t_0$ to $t_f$ (inclusively) over which to compute the solution

The function damping plots the approximation of the solution $y(t)$ and returns a 2D Numpy array with 2 columns:

column at index 0 is the array of $N$ evenly spaced $t$ values from $t_0$ to $t_f$ (inclusively)
column at index 1 is the array of $y$ values of the solution