%% 7. Variable coefficients and adjoints %% % % \setcounter{page}{74} % %% % % In this chapter we focus on second-order problems % with variable coefficients. Locally, near a particular % value of $x$, solutions look % like what we've seen in Chapters 4--6, featuring % oscillation or exponential growth and decay depending on the % coefficients. Globally, over a wider interval of\/ $x$ values, % these structures can combine in interesting ways. % %% % % To begin, let us consider the simplest variable coefficient % of all, one that jumps abruptly from one constant % to another: % $$y'' - 10 \kern 1.3pt \hbox{sign}(x)\kern .5pt y = 0. \eqno (7.1)$$ % For $x<0$ the coefficient is $10$, and for $x>0$ it is % $-10$. Thus this is a second-order ODE with piecewise % continuous coefficients, generalizing the first-order % such problems of Chapter 2, and a solution is defined as a function % $y$ with a continuous derivative that satisfies the equation except % at the points of discontinuity.\footnote{One way to justify this condition is to % think of the discontinuous coefficient ODE as the limit of % smooth ODE\kern .5pt s with sharper and sharper transitions. The second derivative % $y''$ cannot be continuous at $x=0$, since it satisfies (7.1).} % Here is the solution on $[-10,10\kern .3pt]$ with boundary conditions % $y(-10) = y(10)=1$. % ODEformats, L = chebop(-10,10); L.op = @(x,y) diff(y,2) - 10*sign(x)*y; L.lbc = 1; L.rbc = 1; y = L\0; plot(y,CO,bvp), ylim([-2 2]) title('Fig.~7.1.~~Equation with a sign change (7.1)',FS,11) %% % \vskip 1.02em %% % There is no mystery about this picture. % On the left, the ODE is $y'' + 10y = 0$, with oscillatory solutions % $A\sin(\sqrt{10}\kern .7pt x) + B\cos(\sqrt{10}\kern .7pt x)$. % On the right, it is $y'' - 10y = 0$, with exponentially % growing and decaying solutions % $C\exp(\sqrt{10}\kern .7pt x) + D\exp(-\sqrt{10}\kern .7 pt x)$. % The solution makes a continuously differentiable % transition between these behaviors at $x=0$ and matches % the boundary conditions. % Note in particular the exponential boundary layers to the % right of $x=0$ and to the left of $x=10$. % It is not hard to work out this solution exactly by solving four % equations for the four unknowns $A, B, C, D$ (Exercise 7.1). %% % % Abrupt changes of coefficients arise frequently in % applications. Equally important in practice, and mathematically % more interesting, are cases where coefficients change % continuously, and the prototypical ODE of % this kind is the {\bf Airy equation}, % $$y'' - x\kern .1pt y = 0. \eqno (7.2)$$ % Note that $10 \kern 1pt \hbox{sign}(x)$ has given way to simply $x$. % Again we have one sign for $x<0$ and another for $x>0$, % but now with a smooth transition. Here is the solution % with the same domain and boundary conditions as before. % L.op = @(x,y) diff(y,2) - x*y; y = L\0; plot(y,CO,bvp) ylim([-20 20]), title('Fig.~7.2.~~Airy equation (7.2)',FS,11) %% % \vskip 1.02em %% % % Qualitatively speaking, this image is much like the last one. % For $x>0$, the two solutions of % the ODE are exponentially growing and decaying, slowly near % $x=0$ and faster near $x=10$. % For $x<0$, they are oscillatory, with slow oscillations near $x=0$ % and faster ones near $x=-10$. % The point $x=0$ where the equation changes % from one type to another is called a {\bf turning point}. % A notable difference from the last image % is that the amplitude is ten times larger. % This is because, by chance, we are close to an eigenvalue. % If the interval had been $[-9.7,9.7]$ instead % of $[-10,10\kern .3pt]$, the roles would have % been reversed, with the solution to (7.1) having the % high amplitude (Exercise 7.3). % %% % % The Airy equation arises in applications in % optics, quantum mechanics, and other fields in the analysis of % wave effects near interesting edges, including George % Airy's original application in 1838 related to the % colors of the rainbow.\footnote{The Airy equation is also % a prototype in the theory of asymptotics, where it illustrates % fundamental ideas related to steepest descent contours, % asymptotics-beyond-all-orders, and Stokes lines. One can also % think of it as a step toward the subject of {\em multiphysics,} % which is concerned with problems and computations coupling % physical media of different kinds.} % %% % % With almost any ODE, it is generally interesting to reduce % the magnitude of the coefficient of the highest derivative % term to see what happens, and specifically, to sharpen % the separation between local and global behavior. (An ODE % problem with a small leading coefficient is called a % {\bf singular perturbation problem}; see Chapters~18 and~20.) % Here we replace $y''$ by $0.01y''$. % L.op = @(x,y) 0.01*diff(y,2) - x*y; y = L\0; plot(y,CO,bvp) title('Fig.~7.3.~~Airy equation with leading coefficient 0.01',FS,11) %% % \vskip 1.02em %% % % \noindent % On the left, instead of 3 or 4 wavelengths, there are now % about 35, as can be verified by Chebfun: % length(max(y{-10,0},'local')) %% % % \noindent % This increase is by a factor of about $10$, that is, $\sqrt{100}$. % %% % % Like any second-order linear homogeneous ODE, % the Airy equation has a two-dimensional % vector space of solutions. One unique solution % in this vector space (up to % a constant factor) has the special property of decaying to % 0 as $t\to \infty$, and a particular normalization % of this solution is called the {\em Airy function} % $\hbox{Ai}(x)$. The curves we have seen above are almost, but % not quite exactly, multiples of $\hbox{Ai}(x)$, or % $\hbox{Ai}(10 \kern 1pt x)$ for the last example. % As it happens, there is % an integral representation for $\hbox{Ai}(x)$, % $$\hbox{Ai}(x) = {1\over \pi} \int_0^\infty \cos\Big(xs+{s^3\over 3} % \kern 1pt \Big) ds. \eqno (7.3)$$ % It can be shown that at $x=0$, $\hbox{Ai}(x)$ and its derivative % take the values % $$\hbox{Ai}(0) = {3^{-2/3}\over \Gamma(2/3)}, \quad % \hbox{Ai}'(0) = {-3^{-1/3} \over\Gamma(1/3)}, \eqno (7.4)$$ % where $\Gamma$ is the gamma function, and with % these two conditions the solution to (7.2) is % specified fully. The standard choice of % a second linearly independent solution, which goes by the label % $\hbox{Bi}(x)$, is the solution of (7.2) satisfying % $$\hbox{Bi}(0) = {3^{-1/6}\over \Gamma(2/3)}, \quad % \hbox{Bi}'(0) = {3^{1/6}\over \Gamma(1/3)}, \eqno (7.5)$$ % which has the integral representation % $$\hbox{Bi}(x) = {1\over \pi} \int_0^\infty \big( % e^{xs - s^3/3} + \sin(xs + s^3/3)\big)\kern .5pt ds. \eqno (7.6)$$ % Whereas $\hbox{Ai}(x)$ approaches $0$ as $x\to \infty$, % $\hbox{Bi}(x)$ approaches $\infty$. % As $x\to -\infty$, $\hbox{Ai}$ and $\hbox{Bi}$ have asymptotically % the same amplitude, with a phase difference of $\pi/2$, so that % one is approximately zero where the other is % approximately extremal. % We can plot these two solutions by solving the ODE % as a middle-value problem'' (this is not a standard term), with % $y$ and $y'$ prescribed by (7.4) and (7.5) at $x=0$. % L = chebop(-10,10); L.op = @(x,y) diff(y,2) - x*y; g1 = gamma(1/3); g2 = gamma(2/3); Ai0 = 3^(-2/3)/g2; Aip0 = -3^(-1/3)/g1; L.bc = @(x,y) [y(0)-Ai0; feval(diff(y)-Aip0,0)]; Ai = L\0; plot(Ai,CO,green) Bi0 = 3^(-1/6)/g2; Bip0 = 3^(1/6)/g1; L.bc = @(x,y) [y(0)-Bi0; feval(diff(y)-Bip0,0)]; Bi = L\0; hold on, plot(Bi,CO,purple), hold off title('Fig.~7.4.~~Airy functions',FS,11), ylim([-2 2]) text(4.4,0.3,'$\hbox{Ai}(x)$',FS,11) text(1.5,1.25,'$\hbox{Bi}(x)$',FS,11) %% % \vskip 1.02em %% % In applications, turning points of ODE\kern .5pt s may % arise wherever the nature of the physics'' changes. % Though the variable coefficient at such a point may not be exactly % linear, it will usually % be possible to approximate it locally by a linear function. % Consequently the Airy function is relevant to many problems % with turning points. We can illustrate % these effects by changing the variable coefficient of (7.2) % to $\sin(x)$. Specifically, consider % $$0.003 y'' - \sin(x) y = 1, \quad x\in [-4\pi, 4\pi], ~ % y(-4\pi) = y(4\pi) = 0. \eqno (7.7)$$ % The solution is charming. L = chebop(-4*pi,4*pi); L.op = @(x,y) 0.003*diff(y,2)-sin(x)*y; L.lbc = 0; L.rbc = 0; y = L\1; plot(y,CO,bvp), ylim([-40 40]) title('Fig.~7.5.~~Sinusoidal variable coefficient (7.7)',FS,11) %% % \vskip 1.02em %% % % \noindent % Though one could hardly expect to find an % exact formula for such a function, the essence of the matter % is simple enough: oscillation % in regions with $\sin(x)<0$, exponential growth and decay % in regions with $\sin(x)>0$. Locally near each turning point, % $y(x)$ is approximately equal to a suitably scaled % Airy function, and the approximation would get more precise % if the small coefficient $0.003$ were reduced further. % %% % % Airy-like equations feature % transitions between exponential and oscillatory % behavior. When a first-order term is introduced in a % variable-coefficient equation, new possibilities arise. % For example, consider the BVP % $$0.001y'' + x\kern .1pt y' + x\kern .1pt y = 0, % \quad x\in [-2,2\kern .2pt ], ~ y(-2) = -4, ~ % y(2) = 2. \eqno (7.8)$$ % The solution features an {\bf interior layer} at the turning % point $x=0$, which % would grow sharper if the coefficient $0.001$ were further reduced. % L = chebop(-2,2); L.op = @(x,y) .001*diff(y,2)+x*diff(y)+x*y; L.lbc = -4; L.rbc = 2; y = L\0; plot(y,CO,bvp) title('Fig.~7.6.~~Interior layer (7.8)',FS,11) %% % \vskip 1.02em %% % % \noindent % We shall look at this and other problems with boundary % and interior layers in Chapter~20. % %% % % Having explored some ODE\kern .5pt s with variable coefficients, % we now turn to an area of ODE theory that is of particular % interest when the coefficients are variable: % the notion of an {\bf adjoint} and the associated property % of {\bf self-adjointness}.\footnote{Not much in the rest of % the book depends on this material.} Many theoretical properties % are related to adjoints, including whether eigenvalues are % real and whether eigenfunctions are orthogonal. Many practical % matters are tied to adjoints too, % notably the behavior of solutions with respect to perturbations % and the examination of how the outputs of processes governed by % ODE\kern .5pt s depend on the inputs. % %% % % Let $L$ be a second-order linear differential operator acting on % functions on $[-1,1]$. The coefficients may be variable, % so $L$ takes the form % $$L y = a(x) y'' + b(x) y' + c(x) y,$$ % which as usual we will write more simply as % $$L y = a\kern .5pt y'' + b\kern .5pt y' + c\kern .5pt y. \eqno (7.9)$$ % Now we consider the inner product of $v$ and $Ly$, where $v$ and $y$ are % two functions defined on $[-1,1]$. % For full generality, it is necessary in this % business to include the possibility of complex functions, % because even if an operator has real coefficients, its eigenfunctions % and eigenvalues may be complex. Thus the % inner product formula includes a bar over $v$ denoting the complex % conjugate: % $$\int_{-1}^1 \bar v \kern .7pt Ly = \int_{-1}^1 a\kern .5pt \bar v\kern .5pt y'' + b % \kern .5pt \bar v \kern .5pt y' + c \kern .5pt \bar v\kern .5pt y.$$ % (We omit the $dx$ factors for simplicity.) % If we integrate the first and second terms by parts, this % becomes % $$\int_{-1}^1 \bar v \kern .7pt Ly = % \left. \vrule height 10pt width 0pt depth 0pt (a\kern .5pt \bar v % \kern .5pt y' + b\kern .5pt \bar v \kern .5pt y)\right|_{-1}^1 % + \int_{-1}^1 -(a\kern .5pt % \bar v)' y' - (b\kern .5pt \bar v)' y % + c\kern .5pt \bar v\kern .5pt y,$$ % and if we integrate the first term in the integral % by parts a second time, we get % $$\int_{-1}^1 \bar v \kern .7pt Ly = % \left. \vrule height 10pt width 0pt depth 0pt % (a\kern .5pt \bar v \kern .5pt y' -(a\kern .5pt \bar v)'y % + b\kern .5pt \bar v \kern .5pt y)\right|_{-1}^1 % + \int_{-1}^1 (a\kern .5pt\bar v)''y - (b\kern .5pt \bar v)' y % + c\kern .5pt \bar v\kern .5pt y.$$ % In other words, % $$\int_{-1}^1 \bar v \kern .7pt Ly- % \overline{L^*\kern -1pt v}\kern .7pt y = J(\kern .5pt \bar v % \kern .5pt,y), \eqno (7.10)$$ % where $L^*$ is the {\bf formal adjoint} of $L$, defined by % $$L^* \kern -1pt v = (\bar a\kern .5pt v)'' - (\bar b\kern .5pt v)' + % \bar c\kern .5pt v, \eqno (7.11)$$ % and $J(\kern .5pt \bar v \kern .5pt ,y)$ is the bilinear % {\bf conjunct} or {\bf concomitant} % of $\kern .5pt\bar v \kern .5pt$ and $y$, % $$J(\kern .3pt\bar v\kern .5pt ,y) = % \left.\vrule height 10pt width 0pt depth 0pt % (a\kern .5pt \bar v\kern .5pt y' % - (a\kern .5pt \bar v)' y + b\kern .5pt \bar v\kern .5pt y)\ % \right|_{-1}^1 . \eqno (7.12)$$ % In words, {\em $\int_{-1}^1 \bar v\kern .7pt Ly$ is the same as % $\int_{-1}^1 \overline{L^*\kern -1pt v}\kern .7pt y$, apart from % boundary effects.} % Equation (7.10) is known as {\bf Green's formula,} % and it is a univariate prototype of the multivariate % Green's formulas that one finds % in the study of partial differential equations. % %% % % Comparing (7.9) and (7.11), we see that $L^*$ differs from % $L$ in three ways: the functions are conjugated, % the sign of the odd-order derivative % term is negated, and the parentheses have % moved so that $a$ and $b$ are differentiated. % We can write out (7.11) fully as % $$L^*\kern -1pt v = \bar a \kern .5pt v'' + % (2\kern .5pt \overline{a'} \kern .3pt -\kern .3pt \bar b\kern .5pt ) v' + % (\kern .5pt \overline{a''}\kern .5pt -\kern .5pt \overline{b'} % \kern .5pt +\kern .5pt \bar c\kern .5pt )v , \eqno (7.13)$$ % and similarly (7.12) can be expanded as % $$J(\kern .5pt \bar v\kern .5pt ,y) = % \left.\vrule height 10pt width 0pt depth 0pt % a(\kern .5pt \bar v\kern .5pt y' - y\kern .5pt % \overline{v'}\kern .8pt ) + (b-a'\kern .3pt % )\kern .5pt \bar v \kern .5pt y\right|_{-1}^1 . \eqno (7.14)$$ % By taking the adjoint of (7.13), it is readily verified that % the formal adjoint of the formal adjoint is the original operator, % $$(L^*)^* = L .$$ % An operator is {\bf formally self-adjoint} if $L^* = L$, and % from (7.13) it is readily verified that if the coefficient functions are % real, then this is the case if % and only if $a' = b$. % %% % % Thus, for example, if $Ly = y'' + y' + y$, we have % $L^*\kern -.5pt v = v'' - v' + v$ and $J(\kern .5pt \bar v \kern .5pt,y) = % (\kern .5pt \bar v \kern .5pt y'-y\kern .5pt \overline{v'} % +\kern .5pt \bar v \kern .5pt y)\left.\vrule height 8pt width 0pt depth 0pt % \right|_{-1}^1$. % Clearly $L$ is not formally self-adjoint since % $L$ and $L^*$ are different. % Taking the arbitrary choices $v(x) = \exp(x)$ and $y(x) = \hbox{Ai}(x)$, % here is the left side of (7.10): % L = chebop(-1,1); L.op = @(x,y) diff(y,2) + diff(y) + y; x = chebfun('x'); v = exp(x); y = airy(x); v'*(L*y) - y'*(L'*v) %% % % \noindent % A calculation of $J(\kern .5pt \bar v \kern .5pt ,y)$ % confirms that the right side of (7.10) is the same: % Jfun = @(vb,y) vb.*diff(y) - y.*diff(vb) + vb.*y; J = @(vb,y) feval(Jfun(vb,y),1) - feval(Jfun(vb,y),-1); J(conj(v),y) %% % As a similar example with variable coefficients, for the entirely arbitrary % operator $Ly = x\kern .1pt y'' + J_0^{}(x) y' + \sec(x) y$, we have % $L^*\kern -1pt v = x v'' +(2-J_0^{}(x)) v' + (\sec(x) - % J_0'(x) )v$ and $J(\kern .5pt \bar v \kern .5pt ,y) = % x(\kern .5pt \bar v \kern .5pt y'-y % \kern .5pt \overline{v'}\kern .5pt)+(J_0^{}(x)-1) \kern .5pt \bar v \kern .5pt y\left.\vrule height 8pt width 0pt depth 0pt % \right|_{-1}^1$. % Here again is the left side of (7.10), L.op = @(x,y) x*diff(y,2) + besselj(0,x)*diff(y) + sec(x)*y; v'*(L*y) - (L'*v)'*y %% % % \noindent % and here is the right side, % Jfun = @(vb,y) x.*(vb.*diff(y)-y.*diff(vb)) + ... (besselj(0,x)-1).*vb.*y; J = @(vb,y) feval(Jfun(vb,y),1) - feval(Jfun(vb,y),-1); J(conj(v),y) %% % % The word formally'' in the expressions formal adjoint'' and % formally self-adjoint'' means: apart from boundary conditions.'' % To speak of a true adjoint without this qualification, we % need to consider not just an operator but a BVP.\ \ Let $L$ % together with certain homogeneous boundary conditions define a BVP, which % we will write with a script letter as $\cal L$.\ \ The % {\bf adjoint BVP} ${\cal L}^*$ is the BVP in which $L$ is changed to % $L^*$ and the boundary conditions are changed % to {\bf adjoint boundary conditions,} which % are homogeneous boundary conditions on $v$ such that the conjunct % (7.14) is always zero. There exists a unique set of adjoint boundary % conditions with this property, though we shall not prove this. % %% % For example, suppose a BVP for the operator $L$ of % (7.9) has boundary conditions $y(-1)=y(1) = 0$. % Then (7.14) reduces to $J(\kern .5pt \bar v \kern .5pt ,y) = % a(1)\kern .5pt \overline{ v(1)} \kern .5pty'(1) - a(-1)\kern .5pt % \overline{v(-1)} \kern .5pt y'(-1)$. % The values of $y'(1)$ and $y'(-1)$ could be anything, so to ensure % $J(\kern .5pt \bar v \kern .5pt ,y) = 0$, we must assume two conditions: % $$\overline{v(-1)}\kern .5pt =0, \quad \overline{v(1)} % \kern .5pt = 0,$$ % or equivalently, % $$v(-1) =0, \quad v(1) = 0.$$ % These are the adjoint boundary conditions. %% % % For another example, % if the boundary conditions are $y(-1) = y'(1) = 0$, (7.14) % reduces to $J(\kern .5pt \bar v\kern .5pt,y) = % -a(-1)\kern .5pt \overline{v(-1)}\kern .5pt y'(-1) -a(1)y(1) % \kern .5pt \overline{v'(1)}\kern .5pt % + [\kern .5pt b(1)-a'(1)]\kern .7pt \overline{v(1)}\kern .5pt y(1) = 0,$ and this will be zero for all % $y$ if the following two conditions hold: % $$\overline{v(-1)}=0, \quad a(1) \overline{v'(1)}\kern .5pt + % [b(1)-a'(1)]\kern .5pt \overline{v(1)} \kern .5pt = 0,$$ % or equivalently % $$v(-1)=0, \quad \overline{a(1)}\kern .5pt v'(1) + % [\kern .5pt \overline{b(1)}\kern .5pt - \kern .5pt % \overline{a'(1)}\kern .5pt ]\kern .3pt v(1) \kern .5pt = 0.$$ % So these are the adjoint boundary conditions. % In Chebfun, where each chebop may have boundary conditions attached, % the computed adjoint of a chebop is not just the formal adjoint % but the adjoint BVP, with adjoint boundary conditions. % %% % % Suppose $\cal L$ is a BVP and ${\cal L}^*$ is its % adjoint. Then for any functions % $y$ and $v$ satisfying the BCs of $\cal L$ and % ${\cal L}^*$, respectively, the right-hand side of (7.10) will vanish, giving % $$\int_{-1}^1 \kern .5pt \bar v \kern .7pt Ly = % \int_{-1}^1 \overline{L^*v} \kern .7pt y . \eqno (7.15)$$ % A BVP is {\bf self-adjoint} if it is the same as its % adjoint BVP, i.e., % $$\int_{-1}^1 \bar v\kern .7pt Ly= \int_{-1}^1 % \overline{L v}\kern .7pt y, \eqno (7.16)$$ % again for any $v$ and $y$ satisfying the appropriate boundary conditions (which % are now the same for both $v$ and $y$). From this identity we % can derive the following conclusions. % %% % % \smallskip % {\em % {\bf Theorem 7.1. Self-adjoint eigenproblems % (\textsf{\textbf{fLaSHi}}).} % Let $\cal L$ be a self-adjoint % linear homogeneous BVP.\ \ Then % the eigenvalues of\/ ${\cal L}$ are real, and % the eigenfunctions corresponding to distinct % eigenvalues are orthogonal, i.e., if $y$ and $v$ % are eigenfunctions with eigenvalues $\lambda\ne \mu$, % then $\int_{-1}^1 \bar v\kern .5pt y^{} = 0$. % } % \smallskip % %% % % \def\qed{\vrule height4pt width3pt depth4pt} % {\em Proof.} % Let $y$ be an eigenfunction of $\cal L$ with % eigenvalue $\lambda$. Then (7.16) gives % $$\lambda \int_{-1}^1 \bar y \kern .5pt y % = \int_{-1}^1 \kern .5pt \bar y \kern .7pt Ly = % \int_{-1}^1 \overline{Ly} \kern .7pt y = % \bar\lambda \int_{-1}^1 \bar y \kern .5pt y .$$ % Since $y$ is nonzero (part of the definition of % an eigenfunction), the integral of $\bar y \kern .5pt y$ is nonzero, % and thus this identity implies $\lambda = \bar\lambda$, % which establishes the first claim. Now let $v$ be another % eigenfunction of $\cal L$ with eigenvalue $\mu\ne \lambda$. % Then another calculation of the same pattern gives % $$\lambda \int_{-1}^1 \bar v \kern .5pt y % = \int_{-1}^1 \kern .5pt \bar v \kern .7pt Ly = % \int_{-1}^1 \overline{Lv} \kern .7pt y = % \mu \int_{-1}^1 \bar v \kern .5pt y .$$ % Since $\mu \ne \lambda$, this establishes the claim of orthogonality. % \qed % \smallskip % %% % If a self-adjoint operator $\cal L$ has real coefficients, then % the eigenfunctions can be taken to be real, % so the complex conjugates in the calculations above can be dropped. %% % % \begin{center} % \hrulefill\\[1pt] % {\sc Application: adjoints and optimization} \\[-3pt] % \hrulefill % \end{center} % %% % % Whenever something has to be perturbed, optimized, or designed % and there is a differential equation involved, there's a good chance % there may be a role for the adjoint. Here we'll give just a hint % of how this may work. An interesting paper from which to learn % more is M. B. Giles and N. A. Pierce, An introduction to the adjoint % approach to design, {\em Flow, Turbulence and Combustion} 65 % (2000), pp.~393--415. % %% % We show the idea with an % artificial example. Suppose $y$ is the solution of % $$y'' + \sin(x/10) y' + y = f(x) = \exp(-(x-k)^2), \quad % y(0 ) = y(100) = 0 \eqno (7.17)$$ % for some integer $k$ between $1$ and $100$, and we want to know, % which choice of $k$ will % maximize the integral $I(f) = \int_0^{100} y(x)\kern .7pt dx \kern % .7pt ?$ % You can think of this as a prototype of a problem of finding the maximal % lift on an airfoil in various flow conditions. %% % The obvious approach is to simply solve all 100 of these BVPs. % However, this will take a while, as we can see by solving % just the first 10 of them. L = chebop(@(x,y) diff(y,2) + sin(x/10)*diff(y) + y, [0 100]); L.lbc = 0; L.rbc = 0; x = chebfun('x',[0 100]); format compact, format short, tic for k = 1:10 f = exp(-(x-k)^2); y = L\f; disp([k sum(y)]) end toc %% % % \noindent % Use of the adjoint enables us to solve one BVP rather % than 100 of them. The key observation is that since $I(f)$ % is just a single number that depends linearly on $f$ (i.e., % $I$ is a linear functional), it % must be possible to write $I$ in the form % $$I(f) = \int_0^{100} v \kern .2pt f$$ % for some function $v$ defined on $[\kern .5pt 0,100\kern .5pt]$. % (This property of linear functionals % is called the {\em Riesz representation theorem.} % In this application all the functions are real, so we dispense % with complex conjugates.) % In fact, if $g$ is the constant function $1$ on % $[\kern .5pt 0,100\kern .5pt]$, then % $v$ is the solution of the adjoint BVP ${\cal L}^*v = g$. % To derive this characterization, we note that % $$I(f) = \int_0^{100} g\kern .5pt y .$$ % Therefore the two equations % $${\cal L} y = f, \quad {\cal L}^* v = g \eqno (7.18)$$ % imply % $$I(f) = \int_0^{100} g\kern .5pt % y = \int_0^{100} [{\cal L}^* v] y = \int_0^{100} v \kern .5pt % {\cal L} \kern .7pt y = \int_0^{100} v \kern .5pt f . % \eqno (7.19)$$ % By taking advantage of this formula, instead of solving % 100 ODE BVPs, we can now evaluate 100 inner products. % First we confirm that the first 10 of these give % the same numbers as before. % tic, v = L'\1; for k = 1:10 f = exp(-(x-k)^2); disp([k v'*f]) end %% % % \noindent % Next we evaluate all 100 inner products and % plot the maximal solution, which turns out to correspond to $k=32$. % tic, v = L'\1; d = []; for k = 1:100 f = exp(-(x-k)^2); d(k) = v'*f; end [maxint, kmax] = max(d), toc f = exp(-(x-kmax)^2); y = L\f; plot(y,LW,1,CO,bvp), xlabel('$x$',FS,10), ylabel('$y(x)$',FS,10) title(['Fig.~7.7.~~The solution of (7.17) with $k=32$, ',... 'maximizing $I(f)$'],FS,11) %% % \vskip 1.02 em %% % % The function $v$ has an interpretation: it quantifies % the effect on $I(f)$ of the value of $f$ % at each $x\in [\kern .5pt 0,100\kern.5pt]$. % Consequently $v$ is also % the {\em most efficient shape} that $f$ may take, % measured in the 2-norm, % if our goal is to maximize $I(f)$. That is, suppose we % want to find a function $f$ with $\|f\|_2^{} = 1$ such that $I(f)$ % is maximized. % Applying the Cauchy--Schwarz inequality to (7.19) gives % $$I(f) \le \|v\|_2^{} \kern .5pt \|f\|_2^{} ,$$ % with equality if $f$ is a multiple of $v$. So the maximal % value of $I(f)$ for $\|f\|_2^{} = 1$ is achieved % by the choice $f = v/\|v\|_2^{}$, and the value % is equal to $\|v\|_2^{}$, % maxval = norm(v) %% % % \noindent % A plot of $v$ shows that $I(f)$ is % most sensitive to the values of $f(x)$ for $x\approx 32$, as is % consistent with Figure 7.7. % plot(v,LW,1,'k') xlabel('$x$',FS,10), ylabel('$v(x)$',FS,10) title(['Fig.~7.8.~~The adjoint solution $v$ ',... 'governs sensitivity of $I(f)$'],FS,11) %% % % \medskip % \smallskip % {\sc History.} % George Biddell Airy, who used Airy functions in analysis of the % colors of the rainbow in 1838, was one of the leading astronomers % and mathematical physicists of the 19th century. As Astronomer % Royal in Britain during 1835--1881, nearly half a % century, he was largely responsible for the adoption of the longitude % line through Greenwich as the earth's prime meridian. % \smallskip % %% % % {\sc Our favorite reference.} Though there are a hundred % functional analysis books that discuss adjoints in the abstract, % it is hard to find more down-to-earth discussions of adjoint % differential equations in the style of this book. % One source we like is Max D. Gunzburger, {\em Perspectives in % Flow Control and Optimization,} SIAM, 2003. As promised in % the title, this book gives fascinating perspectives on the % meaning of adjoints --- even in the opening 8 pages. % \smallskip % %% % % \begin{displaymath} % \framebox[4.7in][c]{\parbox{4.5in}{\vspace{2pt}\sl % {\sc Summary of Chapter 7.} % Algebraically, second-order % linear ODE\kern .5pt s with variable coefficients have the same structure % as those with constant coefficients: there is a two-dimensional % space of solutions, and an IVP can always be solved but % a BVP will have a unique solution if and only if it does not % correspond to an eigenvalue of the associated homogeneous problem. % Locally, solutions approximate solutions of ODE\kern .5pt s % with constant coefficients, especially when the coefficient of % the second derivative term is small. % Interesting transitions occur at turning points, where a coefficient % changes sign; the prototype is $x=0$ for the Airy equation % $y'' - x\kern .1pt y = 0$. A linear ODE BVP $\cal L$ has an adjoint % ${\cal L}^*$ consisting of the formal adjoint operator % together with adjoint boundary conditions. If $\cal L$ is % self-adjoint, its eigenvalues are real and its eigenfunctions % corresponding to distinct eigenvalues are orthogonal. % \vspace{2pt}}} % \end{displaymath} % %% % % \small\parindent=0pt\parskip=1pt % {\em Exercise $7.1$. Exact solution of the jump problem $(7.1)$.} % {\em (a)} Determine the exact solution of (7.1) with $y(-10)=y(10) = 1$. % {\em (\kern .7pt b)} Determine the eigenvalues of the coresponding % linear operator with boundary conditions $y(-10) = y(10) = 0$. % \par % {\em Exercise $7.2$. An ODE of Euler.} % Consider the ODE $y'' + (a/t)y' + (b/t^2) y = 0$. Show % that it can be reduced to a simpler problem by the substitution % $\tau = \log(t)$ and use this method to find the general solution. % \par % {\em \underline{Exercise $7.3$}. Dependence on interval length.} % {\em (a)} What is the value of the ratio of $\max_x^{} |y(x)|$ for (7.2) divided % by the same quantity for (7.1)? {\em (\kern .7pt b)} What happens % to this ratio if the interval is changed % from $[-10,10\kern .5pt]$ to $[-9.7,9.7\kern .3pt]\kern 1pt$? % \par % {\em \underline{Exercise $7.4$}. Intersection of $\hbox{\rm Ai}(x)$ and % $\hbox{\rm Bi}(x)$.} Use Chebfun {\tt roots} to calculate the % value $x_0^{}$ closest to $0$ for which % $\hbox{Ai}(x_0^{}) = \hbox{Bi}(x_0^{})$. In Matlab you can calculate % $\hbox{Ai}(x)$ and $\hbox{Bi}(x)$ with \verb|airy(x)| % and \verb|airy(2,x)|, respectively. % \par % {\em \underline{Exercise $7.5$}. Airy functions with % many oscillations.} Solve the Airy equation % on the interval $[-D,0\kern .3pt]$ with boundary % conditions $y(-D) = 1$, $y(0)=0$ for % $D = 10, 20, 40, \dots, 320$. Plot the solution in % each case and measure how long it takes Chebfun to % solve the differential equation. For comparison, measure how long % it takes to construct a chebfun for % $\hbox{Ai}(x)$ on $[-32,0\kern .3pt ]$ directly from Matlab's {\tt airy} % command. % \par % {\em \underline{Exercise $7.6$}. Sine oscillations.} % The solution of (7.7) shows seven maxima % in $[-\pi,0\kern .5pt]$ in Fig.~7.5. What does the % count become if the coefficient $0.003$ % is reduced by a factor of $16\kern .7pt$? % \par % {\em Exercise $7.7$. Taylor series solution of Airy equation.} % Exercise 2.10 mentioned a classic analytical solution method for % linear ODE\kern .5pt s, the use of Taylor series. Suppose % $y(t) = \sum_{k=0}^\infty a_k^{} t^k$ is a convergent Taylor series % solution for (7.2) near $t=0$. % {\em (a)} Show that $a_2^{} = a_5^{} = a_8^{} = \cdots = 0$. % {\em (\kern .7pt b)} Give a formula for $a_3^{}, a_6^{}, a_9^{},\dots$ % in terms of $a_0^{}$. % {\em (c)} Give a formula for $a_4^{}, a_7^{}, a_{10}^{},\dots$ % in terms of $a_1^{}$. % {\em (d)} Show that the radius of convergence of the series is % $\infty$, implying that Airy functions are analytic throughout % the complex $t\kern .5pt$-plane. % \par % {\em \underline{Exercise $7.8$}. Eigenvalues of the Airy operator.} % {\em (a)} Calculate the roots of ${\rm Ai}(x)$ in $[-10,10\kern .3pt]$. % {\em (\kern .7pt b)} Calculate the first six eigenvalues % of the Airy operator $L: y \mapsto y'' - x\kern .1pt y$ with % homogeneous boundary conditions on $[\kern .3pt 0,\infty)$. % (It suffices to replace $\infty$ by 15.) % {\em (c)} Explain this coincidence with an % analytical verification. To get the idea it % may help to plot the eigenfunctions. % \par % {\em \underline{Exercise $7.9$}. Legendre equation and polynomials.} % If $n$ is a nonnegative integer, the solution to the BVP % $[(1-x^2)y']' + n(n+1) y = 0,$ $y(-1) = (-1)^n,$ % $y(1) = 1$ has a special property: it is a polynomial, known % as the {\em Legendre polynomial\/} of degree~$n,$ for which % the standard notation is $P_n^{}$. % {\em (a)} Compute $P_{10}^{}$ and $P_{40}^{}$ by solving the % BVP and plot them. You can confirm that your calculation % is successful by comparing these functions with % {\tt legpoly(10)} and {\tt legpoly(40)}. % {\em (\kern .7pt b)} The Legendre polynomials are orthogonal: % $\int_{-1}^1 P_m^{}(x)P_n^{}(x) \kern .5pt dx$ is equal % to $2/(2n+1)$ if $m=n$ and $0$ otherwise. % Confirm this numerically by computing the % inner products of $P_{10}^{}$ with itself, % $P_{10}^{}$ with $P_{40}^{}$, % and $P_{40}^{}$ with itself. (In Chebfun, the % inner product of {\tt u} and {\tt v} is \verb|u'*v|.) % \par % {\em \underline{Exercise $7.10$}. An oscillatory pulse.} % {\em (a)} Plot the solution of $0.01y''-0.01y'-x\kern .1pt y=1$ for $x\in [-10,10\kern .5pt]$ % with $y(-10) = y(10)=0$. Determine (i) the maximum value, (ii) the % width of the highest peak measured at half its height, and (iii) the % distance over which the oscillations decay by a factor of 10. % {\em (\kern .7pt b)} Produce the same plots and numbers in the three alternative cases where % the $y''$ coefficient is reduced to $0.005$, or the $y'$ coefficient, % or both. % \par % {\em Exercise $7.11$. Condition for self-adjointness.} % {\em (a)} % Use (7.13) to verify that if the coefficient functions are real, % then the operator $L$ of (7.9) is formally % self-adjoint if and only if $a' = b$. % {\em (\kern .7pt b)} What is the analogous condition if % the coefficients are complex? % \par % {\em \underline{Exercise $7.12$}.~~Self-adjoint and nonself-adjoint.} % Let $\cal L$ be the BVP $e^xy''+e^x y' + x\kern .2pt y=0$ % with boundary conditions % $y(-1)=y'(1)=0$. % {\em (a)} Show analytically that $\cal L$ is self-adjoint. % {\em (\kern .7pt b)} Compute % the first six eigenvalues and eigenfunctions of $\cal L\kern .3pt$, % showing that the % eigenvalues are real and the eigenfunctions are numerically orthogonal. % {\em (c)} % Show that if the ODE is changed to $e^x y''+ y' + x\kern .2pt y=0$, the eigenvalues % remain real but the eigenfunctions are no longer orthogonal, and verify % that the eigenvalues of $\cal L$ are the same as those of ${\cal L}^*$. % {\em (d)} % On the other hand, it can be shown by an adaptation of the argument of % Theorem~7.1 that the $j\kern .3pt$th eigenfunction of ${\cal L}^*$ is orthogonal % to the $k\kern .3pt$th eigenfunction of $\cal L$ for $j \ne k$. Verify % this numerically. % \par % {\em \underline{Exercise $7.13$}.~~Optimizations with the adjoint.} % {\em (a)} Repeat the optimization illustrated in the Application, but maximizing % $I(f)=\int_0^{100} x\kern .4pt y(x)\kern .7pt dx$. % {\em (\kern .7pt b)} Repeat again, but maximizing % $I(f)=\int_{30}^{50} y(x)\kern .7pt dx$. %