Instructions. Solve the problems for your course track. Formulate your answers clearly and cogently. Sketch out an approach on scratch paper first. Concisely transcribe the approach to the answer you turn in, followed by appropriate calculations and conclusions, within allotted time. Use concise, complete English sentences in the description of your approach. Each question is meant to be completely answered and transcribed from proof to final copy within thirty minutes. Concentrate foremost on clear exposition of the concept underlying your approach. Sign each page thus acknowledging honor code rules. Do not turn in scratch work.

Preamble. The course objectives set out in the syllabus are reproduced below, along with a formulary. Problems probe understanding of the course concepts, rather than rote memorization.

“There are just four mathematics ideas that are considered in this course:

Comparable simplicity is encountered in computational ideas:

Formulary. Newton method to solve $f\left(x\right)=0$: $x_{n+1}=x_n-f\left(x_n\right)/f^{\text{'}}\left(x_n\right)$.

Newton form of interpolating polynomial: $p_n\left(t\right)=\left[y_0\right]+[y_1,y_0](t-x_0)+\cdots +[y_n,\cdots ,y_0](t-x_0)...(t-x_{n-1})$

Divided differences: $\left[y_k\right]=y_k$, $[y_{k+m},...,y_k]=([y_{k+m},...,y_{k+1}]-[y_{k+m-1},...,y_k])/(x_{k+m}-x_k)$

Finite difference operators: $E^{k}f\left(x\right)=f(x+kh)$, $\Delta f\left(x\right)=(E-I)f\left(x\right)$, $\nabla f\left(x\right)=(I-E^{-1})f\left(x\right)$, $\delta f\left(x\right)=(E^{1/2}-E^{-1/2})f\left(x\right)$

Binomial expansion: ${(x+y)}^n={\sum }_{k=0}^n\left(\begin{array}{l}n\\ k\end{array}\right)x^{n-k}{y}^k$, $\left(\begin{array}{l}n\\ k\end{array}\right)=$, $n\in \mathbb{N}$

Generalized binomial expansion: ${(x+y)}^{\alpha }={\sum }_{k=0}^{\infty }\left(\begin{array}{l}\alpha \\ k\end{array}\right)x^{\alpha -k}{y}^k$, $\alpha \in ℝ$

Interpolation operator: $x_k=x_0+kh$, $p_n(x_0+\alpha h)={(I+▵)}^a{y}_0$

Householder reflector $𝑯=𝑰-2$

1Track 1

1. Newton's method to solve $f\left(x\right)=0$, uses data ${𝒟}_n=\left\{(x_n,f_n,f_n^{\text{'}})\right\}$ to construct a linear approximant $g$, and the next iterate $x_{n+1}$ is the solution of $g\left(x\right)=0$. Construct an analogous method based upon data $𝒟=\{(x_{n-1},f_{n-1}),(x_n,f_n,f_n^{\text{'}})\}$, and establish its order of convergence ($f_n=f\left(x_n\right)$, $f_n^{\text{'}}=f^{\text{'}}\left(x_n\right)$).

   Solution. Since three conditions are specified, the interpolating polynomial is quadratic, of form

   ${\displaystyle p\left(t\right)=\left[f_{n-1}\right]+[f_n,f_{n-1}](t-x_{n-1})+[f_n,f_n,f_{n-1}](t-x_{n-1})(t-x_n),}$

   with $[f_n,f_n]=f_n^{\text{'}}$. Divided difference table

   ${\displaystyle \begin{array}{lllll}i& x_i& \left[f_i\right]& [f_i,f_{i-1}]& [f_i,f_{i-1},f_{i-1}]\\ n-1& x_{n-1}& f_{n-1}& -& -\\ n& x_n& f_n& & -\\ n& x_n& f_n& f_n^{\text{'}}& \end{array}}$

   Interpolating polynomial

   ${\displaystyle p\left(t\right)=f_{n-1}+\frac{f_n-f_{n-1}}{x_n-x_{n-1}}(t-x_{n-1})+\frac{f_n^{\text{'}}(x_n-x_{n-1})-(f_n-f_{n-1})}{{(x_n-x_{n-1})}^2}(t-x_{n-1})(t-x_n).}$

   Solve $p\left(t\right)=0$

   ${\displaystyle a{t}^2+bt+c=0,}$ ${\displaystyle a=\frac{f_n^{\text{'}}(x_n-x_{n-1})-(f_n-f_{n-1})}{{(x_n-x_{n-1})}^2},b=d-a(x_{n-1}+x_n),d=\frac{f_n-f_{n-1}}{x_n-x_{n-1}},}$ ${\displaystyle c=f_{n-1}-x_{n-1}d+x_{n-1}{x}_{n}a,}$ ${\displaystyle t_{1,2}=\frac{-b\pm \sqrt{b^2-4ac}}{2a}}$

   Chose $x_{n+1}$ the root closer to $x_n$. Since the interpolation is quadratic, expected order of convergence is cubic.

2. Present a pseudo-algorithm code to compute $x_n=\sqrt{p+x_{n-1}}$ for $n>1$, $n\in \mathbb{N}$, $x_1=\sqrt{p}$. Express the limit $x={lim}_{n\to \infty }{x}_n$ in terms of $p$>0.

   Solution.

   Algorithm

   Input: $n,p$

   $x=\sqrt{p}$

   $\mathrm{for}j=2$ to $n$

   $x=\sqrt{p+x}$

   return $x$

   At the limit $x$ iteration does not change the value, i.e., $x$ is a fixed point. Solve

   ${\displaystyle x=\sqrt{p+x}\Rightarrow x=\frac{1}{2}(1+\sqrt{1+4p})\mathrm{to}\mathrm{ensure}x>0.}$

3. Find a quadrature formula

   ${\displaystyle \underset{-1}{\overset{1}{\int }}f\left(x\right)dx\cong c(f\left(x_0\right)+f\left(x_1\right)+f\left(x_2\right))}$

   that is exact for all quadratic polynomials.

   Solution. Apply moment of methods to obtain system

   ${\displaystyle \underset{-1}{\overset{1}{\int }}1dx=2=3c\Rightarrow c=\frac{2}{3}}$ ${\displaystyle \underset{-1}{\overset{1}{\int }}xdx=0=c(x_0+x_1+x_2)}$ ${\displaystyle \underset{-1}{\overset{1}{\int }}{x}^{2}dx=\frac{2}{3}=c(x_0^2+x_1^2+x_2^2)}$

   Three conditions for four unknowns $c,x_0,x_1,x_2$, one arbitrary degree of freedom. Can impose, e.g.,

   ${\displaystyle \underset{-1}{\overset{1}{\int }}{x}^{3}dx=0=c(x_0^3+x_1^3+x_2^3).}$

   Symmetry suggests $x_2=-x_0$, $x_1=0$, in which formula is exact for cubics also, and

   ${\displaystyle c=\frac{2}{3},x_0=-\frac{1}{\sqrt{2}},x_1=0,x_2=\frac{1}{\sqrt{2}}.}$

4. Determine $A,B$ such that

   ${\displaystyle y_{n+1}=y_n+h(A{f}_n+B{f}_{n+1}),}$

   is a consistent scheme to solve the ordinary differential equation $y^{\text{'}}=f\left(y\right)$.

   Solution. Replace $y_{n+1}=y(x_n+h)$, $f_{n+1}=f\left(y\left(x_{n+1}\right)\right)$, and carry out Taylor series expansions

   ${\displaystyle f\left(y_{n+1}\right)=f\left(y_n\right)+f^{\text{'}}\left(y_n\right)(y_{n+1}-y_n)+\cdots =f\left(y_n\right)+f^{\text{'}}\left(y_n\right)(h{y}^{\text{'}}\left(x_n\right)+\cdots +)}$ ${\displaystyle y\left(x_n\right)+y^{\text{'}}\left(x_n\right)h+y^{\text{'}\text{'}}\left(x_n\right)\frac{h^2}{2}+\cdots =y\left(x_n\right)+h(Af\left(y_n\right)+B[f\left(y_n\right)+h{f}^{\text{'}}\left(y_n\right)y^{\text{'}}\left(x_n\right)+\cdots ])}$

   Gather terms, identify powers of $h$, coefficients are

   ${\displaystyle \begin{array}{ll}𝒪\left(h^0\right)& 0\\ 𝒪\left(h\right)& y^{\text{'}}\left(x_n\right)-(A+B)f\left(y_n\right)\\ 𝒪\left(h^2\right)& \frac{1}{2}{y}^{\text{'}\text{'}}\left(x_n\right)-B{f}^{\text{'}}\left(y_n\right)y^{\text{'}}\left(x_n\right)\end{array}}$

   Since $y^{\text{'}}\left(x_n\right)=f\left(y_n\right)$ the condition

   ${\displaystyle A+B=1}$

   is required for consistency, satisfied, e.g., by $A=1$, $B=0$ by forward Euler, $A=0$, $B=1$ by backward Euler.

5. Determine $a,b,c,d,e$ such that

   ${\displaystyle f\left(x\right)=\{\begin{array}{ll}a{(x-2)}^2+b{(x-1)}^3& x\in (-\infty ,1]\\ c{(x-2)}^2& x\in [1,3]\\ d{(x-2)}^2+e{(x-3)}^3& x\in [3,\infty )\end{array},.}$

   is a cubic spline.

   Solution. Impose continuity of function, derivatives at knots

   ${\displaystyle \mathrm{at}x=1:a=c(\mathrm{repeated}3\mathrm{times})}$ ${\displaystyle \mathrm{at}x=3:c=d(\mathrm{repeated}3\mathrm{times}),}$

   hence $a=c=d$, $b,e$ determined from arbitrary end conditions.

2Track 2

1. Formulate and analyze algorithms for finding roots of $f\left(x\right)=0$, $f:[0,2\pi ]\to ℝ$, $f\in L^2\left(ℝ\right)$, based upon approximation of $f$ within $\mathrm{span}\left(𝒯\right)$, $𝒯=\{1,\cos t,\sin t,\cos 2t,\sin 2t,...\}$.

   Solution. The approximant is a linear combination of functions in $𝒯$

   ${\displaystyle g\left(t\right)=a_0+a_1\cos t+b_1\sin t+\cdots =\sum _{k=0}^{\infty }{a}_k\cos (kt)+\sum _{k=1}^{\infty }{b}_k\sin (kt)}$

   Construct algorithms by analogy to the monomial basis $\{1,t,t^2,\}$, where secant, Newton are obtained as truncations to $\mathrm{span}\{1,t\}$, over an interval $[x_{n-1},x_n]$ with endpoints from a root approximation sequence ${\left\{x_n\right\}}_{n\in \mathbb{N}}$. Hence consider approximant of form

   ${\displaystyle g\left(t\right)=a_0+a_1\cos t+b_1\sin t,}$

   and data $𝒟=\{(x_{n-k},f_{n-k}),k=0,1,2\}$, sorted such that $x_{n-2}⩽x_{n-1}⩽x_n$, and coefficients from solution of

   ${\displaystyle \left[\begin{array}{lll}1& \cos t_{n-2}& \sin t_{n-2}\\ 1& \cos t_{n-1}& \sin t_{n-1}\\ 1& \cos t_n& \sin t_n\end{array}\right]\left[\begin{array}{l}{a}_0\\ a_1\\ b_1\end{array}\right]=\left[\begin{array}{l}{f}_{n-2}\\ f_{n-1}\\ f_n\end{array}\right].}$

   The next iterate $x_{n+1}$ is set as the root of the approximant

   ${\displaystyle g\left(x_{n+1}\right)=a_0+a_1\cos x_{n+1}+b_1\sin x_{n+1}=0,}$

   solvable at the cost of an arctangent evaluation. Algorithm analysis seeks order of convergence of error

   ${\displaystyle e_{n+1}=x_{n+1}-x,}$

   which requires (complicated) Taylor series expansions, but order of convergence should at least equal that of the secant method (superlinear).

2. Present a recursive pseudo-code algorithm to compute to within relative error $\epsilon$

   ${\displaystyle I\left(f\right)=\underset{a}{\overset{b}{\int }}f\left(t\right)dt,}$

   by applying on subintervals the Gauss-Legendre quadrature

   ${\displaystyle \underset{-1}{\overset{1}{\int }}f\left(t\right)dt\cong f(-1/\sqrt{3})+f(1/\sqrt{3}).}$

   Solution. From H09

   Algorithm

   Recursive quadrature

   RecInt(a,b,err,f,Q)

   $c=a+(b-a)/2$

   $Q_{\mathrm{ab}}=Q(f,a,b);Q_{\mathrm{ac}}=Q(f,a,c);Q_{\mathrm{cb}}=Q(f,c,b)$

   $e=|Q_{\mathrm{ac}}+Q_{\mathrm{cb}}-Q_{\mathrm{ab}}|/|Q_{\mathrm{ac}}+Q_{\mathrm{cb}}|$

   if $e<\mathrm{err}$

   return $Q_{\mathrm{ac}}+Q_{\mathrm{cb}}$

   else

   return RecInt(a,c,err,f,Q)+RecInt(c,b,err,f,Q)

   Map $[-1,1]$ to $[a,b]$ through $s\left(t\right)=(b-a)(t+1)/2+a$, $ds=(b-a)/2dt$, hence

   ${\displaystyle Q(f,a,b)=\frac{2}{b-a}[f\left(s_1\right)+f\left(s_2\right)]}$

   with

   ${\displaystyle s_{1,2}=(b-a)(\pm \frac{1}{\sqrt{3}}+1)/2+a.}$

3. Differentiation of the interpolation operator

   ${\displaystyle \frac{d{p}_n\left(x\right)}{dx}=\frac{d}{dx}{(I+\Delta )}^{\alpha }{y}_0=\frac{1}{h}\frac{d}{d\alpha }{(I+\Delta )}^{\alpha }{y}_0=\frac{1}{h}\log (I+\Delta ){(I+\Delta )}^{\alpha }{y}_0=\frac{1}{h}\log (I+\Delta )p_n\left(x\right),}$

   allows identification of the differentiation operator as a forward difference series

   ${\displaystyle \frac{d}{dx}=\frac{1}{h}\log (I+\Delta )=\frac{1}{h}(\Delta -\frac{{\Delta }^2}{2}+\frac{{\Delta }^3}{3}-\cdots +).}$

   Find an analogous expression for the integration operator

   ${\displaystyle \int dx,}$

   and analyze the resulting quadrature formulas.

   Solution. Numerical quadrature formulas are obtained by replacing the integrand with an approximation, in this case a polynomial based upon data set $𝒟=\{(x_i=x_0+ih,y_i=f\left(x_i\right)),i=0,..,n\}$,

   ${\displaystyle \underset{a}{\overset{b}{\int }}f\left(x\right)dx=\underset{x_0}{\overset{x_n}{\int }}{p}_n\left(x\right)dx.}$

   Formal operator integration gives

   ${\displaystyle \underset{x_0}{\overset{x_n}{\int }}{p}_n\left(x\right)dx=\underset{x_0}{\overset{x_n}{\int }}{(I+\Delta )}^{\alpha }{y}_{0}dx=[h\underset{0}{\overset{n}{\int }}{(I+\Delta )}^{\alpha }d\alpha ]y_0=h\left[\frac{{(I+\Delta )}^n-I}{\log (I+\Delta )}\right]y_0.}$

   Carry out series expansions

   ${\displaystyle \frac{1}{\log (1+s)}=s^{-1}+\frac{1}{2}-\frac{s}{12}+\frac{s^2}{24}-\cdots }$ ${\displaystyle {(1+s)}^n-1=ns+\frac{n(n-1)}{2}{s}^2+\cdots +s^n}$ ${\displaystyle \frac{{(1+s)}^n-1}{\log (1+s)}=(s^{-1}+\frac{1}{2}-\frac{s}{12}+\frac{s^2}{24}-\cdots )(ns+\frac{n(n-1)}{2!}{s}^2+\frac{n(n-1)(n-2)}{3!}{s}^3+\cdots +s^n)=}$ ${\displaystyle n(1+\frac{s}{2}-\frac{s^2}{12}+\cdots )(1+\frac{n-1}{2}s+\frac{(n-1)(n-2)}{6}{s}^2+\cdots )=}$ ${\displaystyle n(1+\frac{n}{2}s+\frac{n(2n-3)}{12}{s}^2+\cdots ).}$

   The above leads to

   ${\displaystyle \underset{x_0}{\overset{x_n}{\int }}{p}_n\left(x\right)dx=nh(I+\frac{n}{2}\Delta +\frac{n(2n-3)}{12}{\Delta }^2+\cdots )y_0.}$

   Analyze truncations of the above operator series. Newton-Cotes type formulas are obtained by choosing evaluation points within the integration domain of length $L=x_n-x_0=nh$

   $n=1$.

   ${\displaystyle 𝒪\left({\Delta }^0\right):\underset{x_0}{\overset{x_1}{\int }}f\left(x\right)dx\cong \underset{x_0}{\overset{x_1}{\int }}{p}_n\left(x\right)dx\cong nhI{y}_0=L{y}_0(\mathrm{rectangle}\mathrm{rule},\mathrm{error}𝒪\left(h\right))}$ ${\displaystyle 𝒪\left({\Delta }^1\right):\underset{x_0}{\overset{x_1}{\int }}f\left(x\right)dx\cong \underset{x_0}{\overset{x_1}{\int }}{p}_1\left(x\right)dx\cong h(I+\frac{1}{2}\Delta )=L\frac{y_0+y_1}{2}(\mathrm{trapezoid}\mathrm{rule},\mathrm{error}𝒪\left(h^2\right))}$ ${\displaystyle 𝒪\left({\Delta }^2\right):\underset{x_0}{\overset{x_2}{\int }}f\left(x\right)dx\cong \underset{x_0}{\overset{x_2}{\int }}{p}_2\left(x\right)dx\cong 2h(I+\frac{1}{2}\Delta +\frac{1}{6}{\Delta }^2)=\frac{L}{6}(y_0+4y_1+y_2)(\mathrm{Simpson}\mathrm{rule},\mathrm{error}𝒪\left(h^4\right))}$

4. Construct a Gauss quadrature for integrals of form

   ${\displaystyle \underset{0}{\overset{1}{\int }}{e}^{-t𝑯}𝒇\left(t\right)dt,}$

   that is exact for cubic $𝒇\left(t\right)$, where $𝑯\in {ℝ}^{m\times m}$ is a Householder reflector.

   Solution. In the series

   ${\displaystyle e^{-t𝑯}=𝑰-t𝑯+\frac{t^2}{2!}{𝑯}^2-\frac{t^3}{3!}{𝑯}^3+\cdots ,}$

   note that ${𝑯}^2=𝑰$, i.e., reflection of a reflection is the original vector, verified algebraically by

   ${\displaystyle {𝑯}^2=(𝑰-2)(𝑰-2)=𝑰-4+4=𝑰,}$

   hence

   ${\displaystyle e^{-t𝑯}=(1+\frac{t^2}{2!}+\frac{t^4}{4!}+\cdots )𝑰-(t+\frac{t^3}{3!}+\frac{t^5}{5!}+\cdots )𝑯=\cosh t𝑰-\sinh t𝑯.}$

   The integral becomes

   ${\displaystyle \underset{0}{\overset{1}{\int }}{e}^{-t𝑯}𝒇\left(t\right)dt=\underset{0}{\overset{1}{\int }}(\cosh t𝑰-\sinh t𝑯)𝒇\left(t\right)dt=\underset{0}{\overset{1}{\int }}\cosh t𝒇\left(t\right)dt-𝑯\underset{0}{\overset{1}{\int }}\sinh t𝒇\left(t\right)dt.}$

   The integrals are approximated by Gauss quadrature

   ${\displaystyle \underset{0}{\overset{1}{\int }}\sinh t{f}_i\left(t\right)dt{\cong }_{}{u}_0{f}_i\left(x_0\right)+u_1{f}_i\left(x_1\right),\underset{0}{\overset{1}{\int }}\cosh t{f}_i\left(t\right)dt{\cong }_{}{v}_0{f}_i\left(y_0\right)+v_1{f}_i\left(y_1\right),}$

   using the orthogonal polynomial family $\{{\phi }_0,{\phi }_1,{\phi }_2\}$ with respect to scalar product

   ${\displaystyle (f,g)=\underset{0}{\overset{1}{\int }}\cosh tf\left(t\right)g\left(t\right)dt.}$

   The sample points $x_0,x_1$ (nodes) are roots of ${\phi }_2$, and weights $w_0,w_1$ are found by method of moments. Analogous for the sinh integral.

5. Construct a second-order accurate scheme to solve the integro-differential equation

   ${\displaystyle m{y}^{\text{'}\text{'}}+c{y}^{\text{'}}+ky=\underset{0}{\overset{t}{\int }}R(t-\tau )y\left(\tau \right)d\tau ,}$

   that describes the motion of a point mass $m$ subject to drag $-cv=-c{y}^{\text{'}}$, elastic force $-ky$, and internal relaxation $R\ast y$, starting from initial conditions $y\left(0\right)=y_0$, $y^{\text{'}}\left(0\right)=v_0$.

   Solution. At $t_i=ih$, a centered second-order discretization of the differential operator is

   ${\displaystyle {(m{y}^{\text{'}\text{'}}+c{y}^{\text{'}}+ky)}_i\cong m\frac{y_{i+1}-2y_i+y_{i-1}}{h^2}+c\frac{y_{i+1}-y_{i-1}}{2h}+k{y}_i}$

   A second-order approximation of the integral operator is given by the composite trapezoid rule

   ${\displaystyle \underset{0}{\overset{t_i}{\int }}R(t_i-\tau )y\left(\tau \right)d\tau =\sum _{j=1}^i\underset{t_{j-1}}{\overset{t_j}{\int }}R(t_i-\tau )y\left(\tau \right)d\tau =\frac{h}{2}\sum _{j=1}^i[R_{i-j}{y}_j+R_{i-j+1}{y}_{j-1}]=S_i,}$

   with $R_k=R(kh)$. From $y^{\text{'}}\left(0\right)=v_0$, approximate by $y_1=y_0+h{v}_0$, establishing initial values to advance the iteration

   ${\displaystyle y_{i+1}=2y_i-y_{i-1}+\frac{h^2}{m}[S_i-c\frac{y_{i+1}-y_{i-1}}{2h}-k{y}_i],i=1,2,}$