\documentclass[12pt]{article}
\usepackage{amssymb}
\usepackage{amsmath}
\usepackage{equations}
\usepackage{own}
%\setlength{\textwidth}{6.5in}
%\renewcommand{\baselinestretch}{1.1}
%\renewcommand{\theequation}{\thesection.\arabic{equation}}
% see page 19 of Advanced LaTeX course by Michel Goossens for nifty picture
%\voffset -1in
%\hoffset -1in
\oddsidemargin 0.55625in % Left margin on odd-numbered pages.
\evensidemargin 0in % Left margin on even-numbered pages.
\reversemarginpar % put notes in opposite margin than defaults
\textheight 8.70in
\addtolength{\textheight}{-\baselineskip}
\textwidth 5.98in
\topmargin -.28in % distance from top of page to header - 1in
\headheight 0.5in % height of box containing the header
\headsep \baselineskip % distance between head and body
%\footheight 0.5in % height of foot
\footskip 0.5in % bottom of last line to bottom of foot
\topsep=0in % get rid of extra vertical space in list modes
\partopsep=0in %
\parsep=0.25in %
\parskip 0.1 in
\title{Perturbation Methods\\
MATH 6620--1 -- Spring 2010 \\
Homework 1}
\date{Due Thursday, February 25 at 2:00 PM}
\renewcommand{\theenumi}{\alph{enumi}}
\newcommand{\mathi}{\mathrm{i}}
\newcommand{\expe}{\mathrm{e}}
\newcommand{\difd}{\mathrm{d}}
\newcommand{\grad}{\boldsymbol{\nabla}}
\newcommand{\cross}{\boldsymbol{\times}}
\newcommand{\dive}{\grad \cdot}
\newcommand{\bx}{\mathbf{x}}
\newcommand{\bv}{\mathbf{v}}
\newcommand{\bL}{\mathbf{L}}
\newcommand{\bu}{\mathbf{u}}
\newcommand{\bxinit}{\mathbf{x}_{\text{in}}}
\newcommand{\bvinit}{\mathbf{v}_{\text{in}}}
\newcommand{\rinit}{r_{\text{in}}}
\newcommand{\thinit}{\theta_{\text{in}}}
\newcommand{\rp}{r^{\prime}}
\newcommand{\ssp}{s^{\prime}}
\newcommand{\xip}{\xi^{\prime}}
\newcommand{\tp}{t^{\prime}}
\newcommand{\hax}{\hat{\bx}}
\newcommand{\tilp}{\tilde{p}}
\begin{document}
\maketitle
This homework has 155 points plus 60 bonus points available but, as always, homeworks are graded out of 100 points. Full credit will generally be awarded for a solution only if it is both correctly and efficiently presented using the techniques covered in the lecture and readings, and if the reasoning is properly explained. If you used software or simulations in solving a problem, be sure to include your code, simulation results, and/or worksheets documenting your work. If you score more
than 100 points, the extra points do count toward your homework total.
While you are welcome to use computational tools to automate routine calculations (again, providing a printout of your work), on this homework please do not use software to automate the development of perturbation methods.
Part of the credit is showing that you understand the mechanics of how to develop these perturbation methods.
\section{Perturbing Kepler (55 points)}
This problem concerns the development of an approximate explicit solution for the orbit of one body about another under purely gravitational interaction. For simplicity, we assume one of the bodies is so massive that it may be treated as fixed, and therefore we only follow the dynamics of the lighter body (think of a planet orbiting about a star). Actually this approximation can easily be dispensed with by the concept of reduced mass, but we won't pursue this here.
We'll begin with some background for where the implicit equations of motion arise, and the problem is to convert them to explicit form under the assumption of small eccentricity. We start with Newton's equations for motion for the position $ \bx $ and velocity $ \bv $ of the orbiting lighter body as a function of time $ t $:
\begin{align}
m \frac{\difd \bv}{\difd t} &= - \grad U (|\bx|) \label{eq:newton} \\
\frac{\difd \bx}{\difd t} &= \bv. \nonumber
\end{align}
Here we have the gravitational potential
\begin{equation}
U(r) = -G M m/r, \label{eq:gravity}
\end{equation}
with $ G = 6.67 \times 10^{-11} \mathrm{m}^3/\mathrm{s}^2 \mathrm{kg}$ the universal gravitational constant, $ M$ the mass of the heavier fixed body and $ m $ the mass of the lighter orbiting body. We supplement these differential equations with some initial data
\begin{align}
\bx (t=0) &= \bxinit, & \bv (t=0) &= \bvinit.
\end{align}
One checks that the angular momentum $ \bL = m \bx \cross \bv $ is conserved by these dynamics ($ \difd \bL/\difd t = 0 $), which implies in particular that the motion takes place forever in the plane perpendicular to this conserved angular momentum vector $ \bL $. Taking polar coordinates $ (r,\theta) $ in this plane, the equations of motion become:
\begin{subequations}
\begin{align}
m \frac{\difd^2 r}{\difd t^2} &= - U^{\prime} (r) + \frac{L^2}{m r^3}, \\
m r^2 \frac{\difd \theta}{\difd t} &= L \label{eq:thetat}
\end{align}
\end{subequations}
where $ L = |\bL| $ and $ U^{\prime} (r) = \frac{\difd U (r)}{\difd r} $. These equations can be integrated by noting that they conserve energy $ E = \frac{m}{2} \left(\frac{\difd r}{\difd t}\right)^2 + \frac{L^2}{2 m r^2} + U(r) $. Note that this energy $ E < 0 $ for a bound orbit. Solving this equation for $ \frac{\difd r}{\difd t } $, we obtain a nonlinear first order equation for $ r(t) $:
\begin{equation}
\frac{\difd r}{\difd t} = \sqrt{\frac{2}{m} (E - U(r)) - \frac{L^2}{m r^2}}, \label{eq:rt}
\end{equation}
which can be solved implicitly by separation of variables:
\begin{equation}
t = \int_{\rinit}^{r} \frac{\difd \rp}{\sqrt{\frac{2}{m} (E - U(\rp)) - \frac{L^2}{m \rp{}^2}}} \label{eq:tr}
\end{equation}
where $ \rinit = |\bxinit| $.
Passing to phase plane description implied by the first order equations (\ref{eq:thetat}) and (\ref{eq:rt}) for $ r $ and $ \theta $, we obtain
\begin{equation*}
\frac{\difd \theta}{\difd r} = \frac{L}{mr^2 \sqrt{{\frac{2}{m} (E - U(r)) - \frac{L^2}{mr^2}}}},
\end{equation*}
which can also be solved by separation of variables:
\begin{equation*}
\theta = \thinit + \int_{\rinit}^{r} \frac{L \difd \rp}{{m\rp{}^2 \sqrt{{\frac{2}{m} (E - U(\rp)) - \frac{L^2}{m\rp{}^2}}}}}
\end{equation*}
where $ \thinit $ is the initial angular position. Using now the explicit expression (\ref{eq:gravity}) for the gravitational potential, this integral can be computed through a change of variables to $ \ssp = 1/\rp $, completion of the square in the radical, and a trigonometric substitution:
\begin{subequations}
\begin{equation}
\frac{p}{r} = 1 + e \cos (\theta - \theta_c) \label{eq:rtheta}
\end{equation}
where we have introduced the constants
\begin{equation*}
p = \frac{L^2}{G M m^2},
\end{equation*}
which is half of what is called the \emph{latus rectum} (one way of describing the size of the orbit)
and
\begin{equation*}
e = \sqrt{1 + \frac{2 E L^2}{G^2 M^2 m^3}}
\end{equation*}
is the \emph{eccentricity} (how far the orbit deviates from a circular shape ($e = 0$)).
$\theta_c $ is some constant depending on the initial angular position. One can show that if
$ E < 0 $, the orbit is generally an ellipse (or a circle as a special case when $ e =0 $) with major axis $ a = p/(1 -e^2) $ and minor axis $ b= p/\sqrt{1 - e^2} $.
Likewise, the using the explicit form of the gravitational potential (\ref{eq:gravity}), the integral in Eq.~(\ref{eq:tr}) may be computed through the substitution $ \rp = a (1 - e \cos \xip) $ to give the following parametric solution for how the distance $ r $ of the orbiting body from the fixed body depends on time $ t $:
\begin{align}
r &= a (1 - e \cos \xi), \label{eq:rxi}\\
t &= \sqrt{\frac{a^3}{GM}} (\xi - e \sin \xi) + t_c \label{eq:txi}
\end{align}
\label{eq:par}
\end{subequations}
where $ t_c $ is a constant depending on the initial position of the orbiting body.
So by varying the parameter $ \xi $, we can trace out parametric plots of $ r = r(\xi) $, $ t = t(\xi) $, and $ \theta = \theta (r(\xi)) $ from Eqs.~(\ref{eq:par}). These parametric plots would describe exact solutions to the equations (\ref{eq:newton}).
Suppose though that we wanted to derive explicit expressions for the time-trajectories $ r (t) $ and $ \theta (t) $
of the orbiting body. In fact, this can be approximately done when the eccentricity $ e $ is suitably small.
\begin{enumerate}
\item (\textbf{5 points}) Set up an exactly solvable base problem from which you will develop a perturbation approximation for $ r(t) $ and $ \theta (t) $. Note you do not need to go back to the original differential equation (\ref{eq:newton}); you may freely use any of the above results (which made no assumptions about $ e $) without rederiving them. Derive the formulas for $ r(t) $ and $ \theta (t) $ in your base problem.
\item (\textbf{5 points}) From the exact solution of your base problem, read off scales to normalize (intelligently nondimensionalize) the equations you will use to determine perturbation approximations for $ r(t) $ and $ \theta (t) $ for small eccentricity. Identify the small nondimensional parameter that will be used to generate your perturbation approximation.
\item (\textbf{15 points}) \label{part:pertalg} Obtain an approximate expression for $ r(t) $ and $ \theta (t) $ which includes a leading order term plus a nontrivial correction and an error estimate.
\item (\textbf{10 points}) Plot your approximate solution for $ r(t) $ and $ \theta (t) $ as well as the exact solutions described parametrically in Eq.~(\ref{eq:par}). Show graphically that your approximation improves as your nondimensional small parameter becomes smaller.
\item (\textbf{10 points}) Actually a somewhat awkward but explicit solution for $ r(t) $ and $ \theta (t) $ can be obtained by observing that Eq.~(\ref{eq:txi}) can be solved to give
\begin{equation}
\xi (t) = \sqrt{\frac{GM}{a^3}} (t - t_c) + 2 \sum_{n=1}^{\infty} \frac{1}{n} J_n (n e)
\sin\left(n \sqrt{\frac{GM}{a^3}} (t-t_c)\right)
\end{equation}
where
\begin{equation}
J_n (z) = \sum_{k=0}^{\infty} \frac{(-1)^k}{k!} \frac{(z/2)^{n+2k}}{(k+n)!}
\end{equation}
are Bessel functions. By substituting this exact expression for $ \xi(t) $ into (\ref{eq:rxi}), we obtain an explicit solution $ r(t) $, and substituting this into Eq.~(\ref{eq:rtheta}), we can extract an explicit expression for $ \theta (t) $. Show that your approximation is consistent with this exact solution.
\item (\textbf{10 points}) Derive an approximation for $ r(t) $ and $ \theta (t) $ by a conceptually distinct method than you used in part~\ref{part:pertalg}.
\end{enumerate}
\section{Nonlinear Beam Deflection (55 points)}
Supposedly, an equation for the steady displacement $ u(x) $ of a beam under a sinusoidal loading is given by the equation
\begin{equation}
A \frac{\difd^4 u}{\difd x^4} - B \left[\int_{0}^L \left(\frac{\difd u}{\difd x}\right)^2 \, \difd x\right] \frac{\difd^2 u}{\difd x^2}
= f \sin (\pi x/L )
\label{eq:beam}
\end{equation}
where $ A $ is the product of the elastic modulus and second moment of area, $ B $ is somehow related to these quantities too, $ f $ describes the strength of the loading, and $L $ is the length of the beam. This is supposed to be a simplification of the full beam equation which is OK if the loading is not too large, but honestly I don't see how to derive it. Taking it for granted, though, the appropriate boundary conditions for pinned ends at $ x = 0, L $ are:
\begin{equation}
u(0) = u (L) = \frac{\difd^2 u}{\difd x^2} (x=0) = \frac{\difd^2 u}{\difd x^2} (x=L) = 0.
\end{equation}
Consider now the problem in a situation where the loading magnitude $ f $ is somehow small.
\begin{enumerate}
\item (\textbf{5 points}) Set up a base problem that is exactly solvable but has a nontrivial solution. Note that in this case you can't just take the $ f \rightarrow 0$ limit because that leaves a base problem which just has a trivial zero solution. So you can't neglect the forcing in the base problem, but how does the equation simplify if the forcing is very small? That is, what simpler approximate equation should be valid to leading order when the forcing is small enough?
\item (\textbf{5 points}) Solve your simplified base problem exactly.
\item (\textbf{10 points}) Use the solution of your base problem to choose scales for the variables in equation (\ref{eq:beam}), and perform a corresponding normalization (or intelligent nondimensionalization) of the problem.
Identify a nondimensional parameter which should be small for ``small loading." \label{part:ndode}
\item (\textbf{15 points}) Obtain an approximate solution to your nondimensional equation which contains
a leading order term, one nontrivial correction term, and an error estimate.\label{part:pertode}
\item (\textbf{10 points}) Derive an approximate solution to the original (dimensional) equation (\ref{eq:beam}), and explain the conditions under which the approximation should be expected to be valid.
\item (\textbf{10 points}) Use a different approximation method from the one you used in part~\ref{part:pertode} to generate an approximation with a leading order term plus a correction. (Here it's fine to start with the nondimensional equation from part~\ref{part:ndode} and stop with a nondimensionalized approximation.)
\end{enumerate}
\section{Flow Through Corrugated Channel (45 points plus 20 bonus points)}
A corrugated channel is described by the curves $ y = \pm h(x) $ with $ h(x) = W/2 + C \cos 2 \pi x/\ell $, where
$ W $ is the average width of the channel, $ C $ is the size of the corrugations, and $ \ell $ describes the wavelength of the corrugations.
The steady flow of a fluid through this channel can be described (more or less; see part~\ref{part:fixpde}) by the following two-dimensional partial differential equation for the two-dimensional fluid velocity vector field $ \bu (x,y) $:
\begin{equation}
\mu \Delta \bu (x,y) = - p^{\prime} \hax \text{ for } |y| < h(x),
\end{equation}
where $ \mu $ is the constant dynamic viscosity, $ p^{\prime} $ is a constant prescribed pressure gradient, and $ \hax $ is a unit vector along the $ x $ direction. This equation expresses force balance between viscosity and pressure.
The fluid velocity has no-slip boundary conditions along the corrugated channel:
\begin{equation}
\bu (x,y) = 0 \text{ for } y = \pm h(x). \label{eq:noslip}
\end{equation}
The channel extends infinitely far along positive and negative $ x $, and the appropriate boundary conditions are that the fluid velocity remains bounded as $ x \rightarrow \pm \infty $. The objective of this problem is to compute an approximation, when the corrugation is relatively small, for the fluid flux through a cross section
\begin{equation}
q(x) = \int_{-h (x)}^{h(x)} \bu (x,y) \cdot \hax \, \difd y
\end{equation}
\begin{enumerate}
\item (\textbf{10 points}) Set up and solve a base problem appropriate for studying the effects of small corrugations in the channel.
\item (\textbf{10 points}) Use the solution of your base problem to choose scales to normalize (intelligently nondimensionalize) the given problem (with corrugations). What nondimensional small parameter is appropriate for studying the effects of small corrugations?
\item (\textbf{25 points}) Develop a perturbation expansion to obtain an expression for $ q(x) $ that includes a leading order term, one nontrivial correction term, and an error estimate.
\item (\textbf{20 bonus points}) \label{part:fixpde}You might be concerned that the flux $ q(x) $ varies with $ x $, which implies that the fluid is compressible. For a compressible fluid, a rather more complicated set of equations would need to be solved (including a variable for the fluid density). Rather, fix up your perturbation calculation by rather using the real equation for an incompressible fluid flow (useful for most applications):
\begin{align}
\left. \begin{aligned}
\mu \Delta \bu (x,y) &= - p^{\prime} \hax + \grad \tilp (x,y) \\
\dive \bu (x,y) &= 0
\end{aligned} \right\} \text{ for } |y| < h(x),
\end{align}
with no-slip boundary conditions (\ref{eq:noslip}). Here $ \tilp (x,y) $ is a periodic perturbation to the pressure which is needed to enforce the incompressibility of the velocity field.
\end{enumerate}
\section{Operation Counts (40 bonus points)}
This problem arose from a comment from a student regarding a more rigorous comparison of the computational complexity of the series method and method of iteration in solving regular perturbation problems.
\begin{enumerate}
\item (\textbf{20 bonus points}) Consider some class of algebraic or transcendental regular perturbation problems, and count the number of operations needed to obtain an approximation to a certain order of accuracy with the series method and with the method of iteration. Since they have different computational expense, keep track separately of addition/subtraction, multiplication/division, and root-taking (as well as other transcendental operations). Can you make a general inference about which method is more computationally efficient? For maximum bonus credit, your answer should be in theoretical terms, i.e., a general formula for the operational counts for the approximation methods in terms of the type of equation from the class of equations you consider. You can receive a lesser amount of bonus credit by studying particular examples, either theoretically or by computation.
\item (\textbf{20 bonus points}) Apply the same considerations to some class of perturbed differential equations to obtain some conclusions about the relative computational efficiency of the method of iteration or series method.
\end{enumerate}
\end{document}