
%\documentstyle{siam}
\documentclass[final,psfrag,twoside,epsfrag]{article}
\input{zmacros}
\usepackage{amsmath}
\usepackage{amssymb}
\textwidth 6.4in
\textheight 9in
\oddsidemargin -0.3in
\evensidemargin -0.3in
\topmargin -.5in

%\input{zmacros.tex}
\def\be{\begin{enumerate}}
\def\ee{\end{enumerate}}
%\newcommand{\vthin}{\vskip 5pt}


\pagestyle{empty}
\begin{document}
\hfill\vbox{\hbox{MA 780}\hbox{Homework \#4}
\hbox{\bf Due   April  8}}
\vskip .2in
\baselineskip 14pt


\be

\item This problem provides a way in computing some special functions or integrals.  (a): Use the asymptotic (linear approximation,  $i.e., \, p=1$)  and matching approach to compute the singular integral $\int_0^{\pi} \frac{ \sin x}{\sqrt{x}}\, dx $ accurately to the machine precision (that is,  as  accurate as possible).  Note that the true value is  `fresnels(sqrt(2))*sqrt(2*pi) = 1.789662938968290' using Matlab, where {\em fresnels} is  the Fresnel sine  function. 

\vthin

(b): Now use the sequence $\dsp r_n = \frac{\pi}{2^n}$ to calculate the integral: $\dsp I_n = \sum_{k=0}^{n-1} \dsp \int _{r_{k+1}}^{r_k} \frac{ \sin x}{\sqrt{x}}\, dx$ to compute the integral until $\dsp \frac{|I_n - I_{n-1}|} {|I_n|} < tol$. 

\vthin

(c): Compare  the  two methods to see which method that you think is better and why.
{\bf Hint:}  You should call ( `trapez()' in Matlab) or write a composite trapezoidal quadrature formula. 

\vthin

\item  Consider the 'mid-point' and the 'trapezoidal'  quadrature  formulas on the square: $[0, h]\times [0,h]$.  (a) Find the algebraic precision using polynomials in two dimensions including constant, linear, {\em bi-linear}, quadratic and so on. (b): Find  error estimates of the quadrature formulas for approximating $\iint_{R} f(x,y) dx dy$ in terms of $h$. {\bf Hint:} Use the Taylor expansion in 2D. 

\vthin

(c): Repeat the process for the right triangle $T_h$ determined by the three vertices: $(0,0)$, $(h,0)$, and $(0,h)$. 
 

\vthin

\item Use the Monte Carlo method to evaluate the integral $\iint_{\Omega} \sin (k x) \cos(ky) dx dy$, where $\Omega=[-1,\;1]\times[-1,\; 1]$ for $k=1$, $k=10$, and $k=100$. (a): Find the true value of the integration. (b): Implement the Monte Carlo method with $500$ sample points for 10 runs and then average the results and compare with the true values. 
(c): Analyze your computed results and provide some suggestions or conclusions. 

\vthin {\bf Hint:} You need to either transform the domain or the randomly generated numbers. 

\vthin

\item Convert the IVP problem
\eqmno
y''' + \beta (1-y^2) (y')^2 + \cos y =0, \quad y(0)=1, y'(0)=-1, y''(0)=1
\enmno
to a first order system.

\vthin

\ignore{
\item Indicate whether the following methods for the IVP $y'=f(t,y)$ is
explicit or implicit, one step or multi-step methods. Proof the consistency
of each method.
\eqmno
y_{n+1} &=& y_n + \fr{h}{2} \lp\frac{\null}{\null} f_n + f(t_n +h, y_n+ h f_n)  \rp \\
y_{n+1} &=& y_n + \fr{h}{2} \lp 3 f_n- f_{n-1} \rp
\enmno
where $f_n=f(t_n, y_n)$. }

%\pagebreak }

\item % From Problem 10, Chapter 11 of the text book.
Consider one of two linear multi-step methods ( choose one)
\begin{eqnarray*}
    y_{n+1} &=& \dsp \alpha y_{n} + \frac{h}{2}\left ( \frac{\null}{\null} 2 (1-\beta)f_{n+1} + 3\beta f_n - \beta f_{n-1} \right )
\\ \eqsp
 \mbox{or} \quad    y_{n+1} &=& \dsp \alpha y_{n-1} + \frac{h}{2}\left ( \frac{\null}{\null} 2 (1-\beta)f_{n+1} + 3\beta f_n - 2 \beta f_{n-1} \right ).
\end{eqnarray*}

\be
\item Analyze consistency and order of the methods; determining the values $\alpha^*$ and $\beta^*$ for which the resulting method has maximal order.
\item Study the zero-stability of the method with $\alpha=\alpha^*$ and $\beta=\beta^*$; write its characteristic polynomial and using MATLAB, draw its region of absolute stability in the complex plane.
\ee

\vthin

\item Consider the two-stage of RK method (it is also called
the improved Euler method):
\eqmno
K_1 &=& f(x_n, y_n) \\
K_2 &=& f(x_{n+1}, y_n + h K_1) \\
y_{n+1} &=&  y_n + \fr{h}{2} \lp K_1 + K_2 \rp
\enmno
\be
\item Show the local truncation error is order of $h^2$.
\item Discuss the zero and absolute stability of the algorithm. Find or
plot the stability region.
\ee

\vthin

 \vthin


\item Referring to mymovie.m and chemistry.m on the class web-page, use Matlab
ODE-suite to solve the following IVP
\eqmno
&& y'_1 = a - (b+1)y_1 + y_1^2 y_2 \\
&& y_2' = b y_1 - y_1^2   y_2.
\enmno
Try the parameters $(a,b)=(1,3), (1,10), (1,20), (5,1)$, $T=50, 500$. Plot
the history of the solution and phase plot. Find all the equilibriums and classify them (stable, unstable, center etc.). Do you think this problem is
stiff? A stiff problem is that the eigenvalues of the linearized system
have very different magnitude.

\vthin

\item {\bf Extra Credit:}  This problem  illustrates how we can combined different methods  to solve problems. 
 Implement the Crank-Nicolson scheme for the Lorenz equations using a non-linear system of equations solver or your own code. 
\eqmno
&& x'  = -\tau x + \tau y \\
&& y' =  - x z + r x -  y  \\
&& z' = xy - \beta z.
\enmno
Try you method at least with one set of data $r=100, \tau=10, \beta=8/3$, final time $T=20$ and an initial 
data $[1,\;0,\;0]$. Plot the history of the solution $x(t), y(t), z(t)$, and phase plot. 

\vthin

Apply one step Richardson extrapolation to the final solution and compare with that of direct computation and provide a posterior error estimate.  

\vthin

Try several typical sets of data. Write a report on your algorithm and analysis,   compare your numerical results with those from literature;  and try to interpret your results. 


\ee
\end{document}
