Numerical Analysis I

Problems

1. Write a Python function called CompositeSimpson that takes as input an array $y=(y_0,y_1,\dots ,$ $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,\dots ,n$ and $h=\frac{b-a}{n}$. Your function should first check that $n$ is even and $n\ge 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\Vert {\gamma }^{\prime }(t)\Vert dx={\int }_a^b\sqrt{(x^{\prime }(t){)}^2+(y^{\prime }(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{array}{rl}x(t)& =r_1\cos (t)\\ y(t)& =r_2\sin (t)\end{array}$$ 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 $\epsilon =1\times {10}^{-6}$. Hint: Here $f(t)=\sqrt{(x^{\prime }(t){)}^2+(y^{\prime }(t){)}^2}$; use math software to compute $f^{(4)}(t)$ to estimate $M=\underset{t\in [0,2\pi ]}{max}|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^{\prime }(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^{\prime }(p)\ne 0$. Let $$f(x)=\frac{1}{\sqrt{2\pi }}{\int }_0^x{e}^{-t^2/2}dt-0.45$$ and thus $$f^{\prime }(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.