Numerical Analysis I

Problems

1. Write a Python function called CompositeSimpson that takes as input an array \(y=(y_0, y_1, \ldots,\) \(y_n)\) and a step-size \(h\), and returns the Composite Simpson approximation of the integral \(\int_a^b f(x)\, dx\) where \(y_j = f(x_j)\) and \(x_j = a + jh\) for \(j=0,1,\ldots, n\) and \(h=\frac{b-a}{n}\). Your function should first check that \(n\) is even and \(n\geq 2\) and returns an error message if not both conditions are satisfied. For example,

      print('Error: n must be even!')
      return 0

   Test your function by choosing your all-time favorite integral \(\int_a^b f(x)\,dx\) that you can compute exactly.

2. Recall from calculus that if \(\gamma(t) = (x(t), y(t))\) is a parametrized curve then the arc length of \(\gamma\) from \(\gamma(a)\) to \(\gamma(b)\) is given by the integral \[L = \int_a^b \|\gamma'(t)\|\, dx = \int_a^b \sqrt{(x'(t))^2 + (y'(t))^2}\, dt.\] In even the simplest cases, the above integral has no closed form and numerical integration is necessary. As a simple example, a parametrization of the ellipse \(\frac{x^2}{r_1^2} + \frac{y^2}{r_2^2} = 1\) is given by \[\begin{aligned} x(t) &= r_1\cos(t)\\[2ex] y(t) &= r_2\sin(t) \end{aligned}\] for \(t\in [0, 2\pi]\), and where \(r_1>0\) and \(r_2>0\). If \(r_1=3\) and \(r_2=5\), find an estimate of the arc length of the ellipse using the Composite Simpson approximation accurate to within \(\varepsilon = 1\times 10^{-6}\). Hint: Here \(f(t) = \sqrt{(x'(t))^2 + (y'(t))^2}\); use math software to compute \(f^{(4)}(t)\) to estimate \(M=\max_{t\in [0,2\pi]} |f^{(4)}(t)|\) and use your estimate for \(M\) to find \(n\) as we did in the example in class. Wolfram would be able to compute the 4th derivative and graph it.

3. Recall that Newton's method is a numerical method to obtain an approximation to a root of an equation $f(x)=0$ using fixed-point iteration with the function \[ g(x) = x - \frac{f(x)}{f'(x)}. \] In fixed-point iteration, we generate the sequence $p_{k+1} = g(p_k)$ by supplying an initial condition $p_0$ and under certain conditions the sequence $(p_k)$ converges to a fixed-point $p$ of $g$ and thus a root of the equation $f(x)=0$ provided $f'(p)\neq 0$. Let \[ f(x) = \frac{1}{\sqrt{2\pi}} \int_0^x e^{-t^2/2}\, dt - 0.45 \] and thus \[ f'(x) = \frac{1}{\sqrt{2\pi}} e^{-x^2/2}. \] To evaluate $g(p_k)$ we need to evaluate the integral \[ \frac{1}{\sqrt{2\pi}} \int_0^{p_k} e^{-t^2/2}\, dt \] which can be done using the Composite Simpson rule. Find an approximate solution to $f(x)=0$ by performing $N=200$ iterations of Newton's method with initial condition $p_0 = 0.5$. What is your approximation $\hat{p}$ and what is $f(\hat{p})$? Suggestion: Do not use your previous fixed-point iteration/Newton Python code, it is easier to write the fixed-point iteration for loop directly and using your CompositeSimpson function.