%% 21. Into the complex plane %% % % \setcounter{page}{265} % %% % Most of the familiar mathematical functions are defined for % complex arguments as well as real ones. For example, here are the % Taylor series of two functions about the point $t_0=0$: % $$e^{\kern .3pt t} = 1 + t + {t^2\over 2} + {t^3\over 3!} + \cdots, \quad % {1\over 1+t^2} = 1 - t^2 + t^4 - t^6 + \cdots.$$ % To define the value of a function when $t$ is complex, % we can use the same series, so long as $t$ is close % enough to $t_0$ for convergence. %% % % We say that a function $y(t)$ is {\bf analytic} at a point % $t_0$ if it has a Taylor series at~$t_0$ that converges to $y$ in some % neighborhood of $t_0$. If $y$ is initially defined just for real values, % then the neighborhood will be an interval $(a,b)$ with % $a < t_0 < b$. Once we've got a convergent power series, it can % be used to define values in the complex plane too, % and this is the process called % {\bf analytic continuation}. It is a basic fact % of complex variables that every power series has a % {\bf radius of convergence} $r\in [\kern .3pt 0,\infty]$ with the property % that the series converges for all real or complex % numbers $t$ with $|\kern .5pt t-t_0|r$. Thus a power % series always converges inside a disk, the % {\bf disk of convergence,} which is the largest disk % centered at $t_0$ inside which the function is analytic. % %% % % For $y(t) = e^{\kern .3pt t}$, the disk of convergence has radius % $r=\infty$, regardless % of the choice of $t_0$, because the exponential function is analytic % for all values of $t$. A function $y(t)$ that % is analytic for all real and complex values % of $t$ is said to be {\bf entire}. For $1/(1+t^2)$ and % $t_0=0$, we have $r=1$, because the function has % singularities at $t=\pm i$ (poles), where the value blows % up to $\infty$. For $t_0=1$ with the same function, we % would have $r=\sqrt 2$. % %% % % It is a general idea of mathematics that one can learn % things about a function by examining its singularities in the % complex plane. For example, how would you explain to % a student in a calculus class {\em why\/} the function % $y(t) = 1/(1+t^2)$ has a Taylor series % that converges just for $t\in (-1,1)$? The answer % is the presence of those singularities at $\pm i$, and % if your student doesn't know about the complex plane, it % is not clear how you could really give a satisfactory explanation. % %% % % Here is another example. The function $y(t) = \tanh(10\kern .7pt t)$ % makes a rapid transition from $-1$ to $+1$ near $t=0$. % ODEformats, y = chebfun('tanh(10*t)'); plot(y), ylim([-1.5 1.5]) title('Fig.~21.1.~~A function $y(t)$ with a rapid transition',FS,11) %% % \vskip 1.02em %% % % \noindent % How does this transition come about? % One way to understand it is to note that~$y$ has poles % in the complex plane quite near to $t=0$, at $\pm \pi i/20,$ $\pm 3\pi i/20,$ % $\pm 5\pi i/20,\dots.$ % As $t$ increases along the real axis, $y$ % changes quickly as it passes these special points. % A contour plot of $|y(t)|$ for $t$ ranging over a rectangular region % in the complex plane reveals the four poles of $y$ % closest to the origin. % x = linspace(-1.2,1.2,201); y = linspace(-0.5,0.5,101); [xx,yy] = meshgrid(x,y); zz = xx + 1i*yy; ff = tanh(10*zz); levels = [1 2 4]; contour(x,y,abs(ff),levels,'k',LW,.7), axis equal hold on, plot(pi*.05i*(-3:2:3),'.r',MS,8), hold off title(['Fig.~21.2.~~In the complex $t$-plane, ' ... 'there are poles nearby'],FS,11) %% % \vskip 1.02em %% % % \noindent % If the poles were ten times closer, $y$ would make its transition ten times faster. % %% % For another example, what can we make of the function % $$y(t) = -\log (1.1 + \cos(\pi t))? \eqno (21.1)$$ % A plot reveals spikes near odd integer values of $t$. f = chebfun('-log(1.1+cos(pi*t))',[-4 4]); plot(f), ylim([-3 3]) title(['Fig.~21.3.~~Another function $y(t)$ ' ... 'with regions of rapid transition'],FS,11) %% % \vskip 1.02em %% % % \noindent % As before, these points of rapid change are % associated with nearby singularities in the complex plane, which % lie in pairs just above and below each odd integer. % x = linspace(-4,4,201); y = linspace(-1.5,1.5,101); [xx,yy] = meshgrid(x,y); zz = xx + 1i*yy; ff = -log(1.1+cos(pi*zz)); levels = 2.5:.5:7; contour(x,y,abs(ff),levels,'k',LW,.7), axis equal hold on, plot(acos(1.1)/pi+(-3:2:3),'.r',MS,8) plot(-acos(1.1)/pi+(-3:2:3),'.r',MS,8), hold off title(['Fig.~21.4.~~Now the nearby complex ' ... ' singularities are branch points'],FS,11) %% % \vskip 1.02em %% % % \noindent % This contour plot requires more interpretation than the previous % one. Before, each singularity $t_s$ was a % {\bf pole,} with $y\sim C/(t-t_s)$ for some constant $C$. % In this second % contour plot, however, the two rows of singularities % are not poles but {\bf branch points}. A function with a branch % point is not analytic and single-valued in any punctured % neighborhood of the point, because if you continue it % analytically all the way around, you reach a different value from % the one you started with. % To make a function $y$ with branch points single-valued and analytic, % it is customary to restrict the domain by introducing {\bf branch cuts} connecting the % branch points, along which $y$ is not defined. Given a function % with certain branch points, there is no unique % set of branch cuts that the mathematics forces you to choose. % Computer programming languages, however, generally % take $\log(t)$ and $t^{1/2}$ to have % branch lines extending along the negative real axis from $0$ to % $-\infty$, and this determines where the branch cuts lie for % more complicated functions like (21.1). % In the figure, the branch lines have ended up % as vertical rays extending upward to $\infty$ from the % branch points in the upper half-plane % and downward to $\infty$ from the branch points in the lower % half-plane. % %% % % Everything we have said in the last three % pages about functions in general applies to % the particular case of solutions of ODE\kern .5pt s. Many of them make sense % in the complex plane, and their behavior in the plane can % both be interesting in its own right and also shed light on % behavior, such as rapid transitions, on the real axis. % %% % As the simplest possible example, consider the IVP % $$y' = y, \quad y(0) = 1 . \eqno (21.2)$$ % The unique solution is $e^{\kern .3pt t}$, an entire function in the complex % $t$-plane, and thus for example at $t=i$ we have % $$y(i) = \exp(i) = \cos(1) + i\kern .3pt \sin(1) % \approx 0.5403 + 0.8415i .$$ % We can interpret (21.2) in the complex plane % in two ways. One is to % regard it as a differential equation that applies not % just for real values of $\kern .3pt t$ but also % complex ones. The value $y(i)$, for example, might % be determined by integrating the ODE along the imaginary line segment % extending from $t=0$ to $t=i$. % Chebfun doesn't work with complex intervals directly, but % we can achieve the necessary effect by parametrizing the % interval $[\kern .3pt 0, i\kern .3pt ]$ as $i s$ for % $s\in [\kern .3pt 0,1]$. The equation becomes % $${dy\over ds} = {dt\over ds} {dy \over dt} = i\kern .3pt y, % \eqno (21.3)$$ % or in Chebfun, L = chebop(0,1); L.lbc = 1; L.op = @(s,y) diff(y) - 1i*y; y = L\0; y(1) %% % \vskip 1.02em %% % The other interpretation of (21.2) for complex $t$ is that % we may start from the solution $y(t)$ for real $t$, and then % analytically continue it into the complex plane. This idea defines the % same function as before, provided one uses the same branch % cuts for the analytic continuation as for the ODE solution. %% % For example, here we solve (21.2) in the usual manner % on the real interval % $[\kern .3pt 0,1]$, and then evaluate the result % at the complex point $t = i$. L = chebop(0,1); L.op = @(t,y) diff(y) - y; L.lbc = 1; y = L\0; y(1i) %% % \vskip 1.02em %% % % \noindent % Since Chebfun's solutions of ODE\kern .5pt s are numerical rather % than symbolic, it is entirely appropriate to wonder % why this worked! % The reason is that Chebfun's solution function $y$ % is actually a polynomial representation of $\exp(t)$ % on $[\kern .3pt 0,1]$, which also approximates the % function in a larger region of the complex plane. % So Chebfun has performed a numerical version of % analytic continuation for us. % As a rule, this trick works reasonably well for complex % values of $\kern .3pt t$ close to the interval of definition, % but not for values further out, and certainly not for % values of $\kern .3pt t$ that lie further from the interval % than some singularities of $y$. % %% % Let us illustrate this by considering an example with a singularity. % The nonlinear IVP % $$y' = -2\kern .3pt t\kern .3pt y^2, \quad t\in [-4,4], % ~~ y(-4) = {1\over 17} \eqno (21.4)$$ % has the unique solution we examined at the start of the chapter, % $$y(t) = {1\over 1+t^2},$$ % with poles at $t = \pm i$. The solution looks as it should, N = chebop(-4,4); N.lbc = 1/17; N.op = @(t,y) diff(y) + 2*t*y^2; y = N\0; plot(y,CO,ivpnl), ylim([0 1.2]) title(['Fig.~21.5.~~Solution of an ODE ' ... 'passing near a complex singularity'],FS,11) %% % \vskip 1.02em %% % % \noindent % At $t=i/2$, the correct value of $y$ is $4/3$, % as Chebfun confirms, % y(1i/2) %% % % \noindent % At $t=3i/2$, however, on the far side of the pole at $i$, % the correct value of $y$ is $-4/5$, whereas % Chebfun gets a huge incorrect value. % y(3i/2) %% % % \noindent % One can never expect a simple evaluation like this to work, % if the point of evaluation is further from the interval of % computation than some of the singularities of the % function.\footnote{This statement can be made mathematically % precise with the use of {\em Bernstein ellipses,} ellipses % in the complex $t$-plane whose foci are at the two endpoints % of the interval of definition of the chebfun. A chebfun % will normally give a degree of good approximation in the largest % region of analyticity bounded by such an ellipse, % as discussed in Chapter~8 of % Trefethen, {\em Approximation Theory and Approximation Practice,} % SIAM, 2013.} % %% % % There is an alternative approach in Chebfun, however, that can often % be used for numerical analytic continuation even beyond singularities, % especially when the singularities are poles. The Chebfun command % {\tt aaa} approximates a function not by the usual % Chebfun method involving polynomials, but by a rational function % $r(t) = p(t)/q(t)$ of some adaptively determined type $(n,n)$, which means % that the numerator and denominator degrees are $\le n$.\footnote{See % Nakatsukasa, S\ete, and Trefethen, The AAA algorithm % for rational approximation,'' submitted manuscript, 2016.} For example, the command % [r,pol] = aaa(y,'tol',1e-8); %% % % \noindent % constructs a rational approximation to $y$ and evaluates its poles, which % match those of $y$, % pol %% % % \noindent % This time, the value at $t = 3i/2$ comes out correct. % r(3i/2) %% % % Although rational approximations will usually not give % precise information, especially % since singularities of solutions to ODE\kern .5pt s in the complex % plane are usually more complicated than just poles, they % are often good at giving a rough idea of the nature of such % singularities near the domain of approximation. For example, in Chapter 13 % we plotted one component of a solution to the Lorenz equations (13.1), % N = chebop(0,5); N.op = @(t,u,v,w) [diff(u)-10*(v-u); ... diff(v)-u*(28-w)+v; diff(w)-u*v+(8/3)*w]; N.lbc = [-15; -15; 20]; [u,v,w] = N\0; plot(u,CO,ivpnl); title('Fig.~21.6.~~A trajectory of the Lorenz equations',FS,11) %% % \vskip 1.02em %% % % \noindent % You might guess that these spikes are associated with nearby % singularities in the complex plane, and here are the poles of the % {\tt aaa} rational approximation. % Note that the sharpest spikes correspond to % singularities closest to the real axis. % [r,pol] = aaa(u,'tol',1e-8); plot(pol,'.r',MS,10), xlim([0 5]), axis equal title(['Fig.~21.7.~~Nearby complex singularities ' ... 'found by rational approximation'],FS,11) %% % \vskip 1.02em %% % % One should not imagine that the function plotted in Fig.~21.6 truly % has poles at all the points shown in Fig.~21.7. Rather, it is likely that % it has a singularity --- probably not just a pole --- near the % point in each group that lies closest to the real axis. % The rest of the dots are probably lining up along branch cuts. % %% % % Solutions to linear ODE\kern .5pt s have no singularities unless the leading-order % coefficient passes through $0$ or % the coefficients are singular, but solutions to % nonlinear ODE\kern .5pt s may have all kinds of singularities, which % may be {\em movable} in the sense that they % appear at locations dependent on initial data. % A natural question for a mathematician is, are there % nonlinear ODE\kern .5pt s whose movable singularities are only % poles, never branch points? It turns out that a % full answer to this question is known in the case of second-order ODE\kern .5pt s % of the form $y'' = F(t,y,y')$, where $F$ is a rational function. % Paul Painlev\'e showed that all such equations possessing % this property that the movable singularities are poles % can be organized into 50 classes, 44 of which % can be reduced to other known functions.\footnote{Painlev\'e % was not your average mathematician. In 1917 and again % in 1925, he served as Prime Minister of France, and he is % buried in the Panth\'eon in Paris.} The remaining % six classes are called the {\bf Painlev\'e equations,} and % their solutions are known as {\bf Painlev\'e transcendents.} % The {\bf Painlev\'e I} equation is % $$y'' = 6y^2 + t, \eqno (21.5)$$ % and the {\bf Painlev\'e II} equation is % $$y'' = 2y^3 + ty + \alpha, \eqno (21.6)$$ % where $\alpha$ is a parameter. % %% % % Let us take a look at (21.5), the Painlev\'e I equation. % In Chapter 3, we saw that % the first-order ODE $y'=y^2$ and various generalizations % have solutions which blow up to $\infty$. % The blowup points are simple % poles. For (21.5), a second-order equation of a similar form, % one gets double poles instead of simple poles. For example, % here is a solution to (21.5) on the interval $[-15,2.27\kern .5pt]$ % with middle conditions'' $y(0) = y'(0) = 0$. For % $t<0$, the solution is oscillatory and nonsingular, but % for $t>0$ the curve is approaching a double pole at $t\approx 2.6$. % N = chebop(0,2.27); N.op = @(t,y) diff(y,2)-6*y^2-t; N.lbc = [0;0]; y = N\0; plot(y,CO,ivpnl) N = chebop(-15,0); N.op = @(t,y) diff(y,2)-6*y^2-t; N.rbc = [0;0]; y = N\0; hold on, plot(y,CO,ivpnl), hold off axis([-15 3 -4 10]) title('Fig.~21.8.~~Solution to Painlev\''e I eq.~(21.5)',FS,11) %% % \vskip 1.02em %% % % \noindent % The reader of this chapter will suspect that if $y(t)$ % oscillates like this on the negative $t$-axis, there must % be singularities nearby in the left half of the complex % $t$-plane. This is true, and in fact, there is an infinite % array of double poles extending to $\infty$ in all directions. % %% % % The Painlev\'e I equation does have smooth solutions for certain % boundary data, which we can isolate most easily by % solving a BVP.\ \ Here is an example: % $$y'' = 6y^2 + t, \quad t \in [-24,0\kern .3pt], ~ y(-24) = 2, ~y(0) = 0. % \eqno (21.7)$$ % N = chebop(-24,0); N.op = @(t,y) diff(y,2) - 6*y^2 -t; N.lbc = 2; N.rbc = 0; y = N\0; plot(y,CO,bvpnl) title('Fig.~21.9.~~Solution to Painlev\''e I eq.~(21.7)',FS,11) %% % \vskip 1.02em %% % % \noindent % Because of the exponential sensitivities introduced by the % nonlinear term, one is hardly likely to find such smooth % solutions by solving an IVP.\ \ To illustrate this % sensitivity, here we superimpose % on the figure the solution to % an IVP (or if you prefer a {\em final-value problem}) % corresponding to the same solution just plotted, except that the boundary % conditions are both specified at the right-hand boundary, with % a small perturbation introduced in the derivative: % $$y'' = 6y^2 + t, \quad t \in [-24,0\kern .3pt], ~ y(0) = 0, ~y'(0)= % 0.999999\alpha. % \eqno (21.8)$$ % Here $\alpha\approx -0.451427$ is the value of $y'(0)$ corresponding to the solution % of (21.7). % yp = diff(y); alpha = yp(0); N.lbc = []; N.rbc = [0; 0.999999*alpha]; y = N\0; hold on, plot(y,CO,ivpnl), hold off title('Fig.~21.10.~~Another nearby solution superimposed',FS,11) %% % \vskip 1.02em %% % % \noindent % The new solution matches the smooth one approximately for % $t\in [-5,0\kern .3pt]$, but then diverges to an oscillatory form. % %% % % It is impossible to understand Painlev\'e equations very fully % by looking just on the real axis, however, as is shown beautifully % by the complex plane explorations of the paper % by Fornberg and Weideman cited as our favorite reference below. % %% % % \begin{center} % \hrulefill\\[1pt] % {\sc Application: Jacobi sine function}\\[-3pt] % \hrulefill % \end{center} % %% % % Let $m\in [\kern .3pt 0,1)$ be a parameter and consider the nonlinear % second-order initial-value problem % $$y'' = -(1+m)y + 2m y^3, \quad y(0) = 0, ~ y'(0) = 1 . % \eqno (21.9)$$ % If $m$ is zero, the equation is just $y'' = -y$ and the % solution is $y(t) = \sin(t)$. As~$m$ increases, the % nonlinearity becomes stronger. Here is what we find for % $m = 0.998$ over the interval $t\in [\kern .3pt 0,100\kern .3pt ]$. % N = chebop(0,100); N.lbc = [0;1]; m = 0.998; N.op = @(y) diff(y,2) + (1+m)*y - 2*m*y^3; y = N\0; plot(y,LW,1,CO,ivpnl), ylim([-2 2]), grid on title(['Fig.~21.11.~~Solution of (21.9): ' ... 'Jacobi sine function $\hbox{sn}(t,0.998)$'],FS,11) %% % \vskip 1.02em %% % % \noindent % Like $\sin(t)$, this function oscillates % between $-1$ and $1$, but its shape is squarer. % The period $T$ is about three times longer than for $\sin(t)$, % reflecting how much time this function spends lingering % near $\pm 1$: % [~,maxima] = max(y,'local'); T = maxima(3) - maxima(2) %% % % Now $\sin(t)$ is completely smooth, analytic throughout % the complex $t$-plane. The function $y(t)$, on the other hand, % makes rapid transitions between $-1$ and~$1$, which will become % more abrupt as $m$ increases toward~1.\ \ We may accordingly guess that if % this function is analytically continued into the complex plane, there will be % singularities above and below the $t$-axis near $t=0, T/2, T, \dots .$ % A call to {\tt aaa} confirms this prediction. In the next figure, % the dots show the poles of the {\tt aaa} rational approximant, and % the vertical lines mark ${\hbox{Re}\kern .3pt}(t) = 0, T, \dots , 6T$. % [r,pol] = aaa(y,'tol',1e-8); plot([-10 110],[0 0],LW,.7,CO,.75*[1 1 1]), grid off, xlim([-10 110]), axis equal, hold on for k = 0:6 plot(k*T*[1 1],[-30 30],LW,.7,CO,.75*[1 1 1]) end plot(pol,'.r',MS,10) title(['Fig.~21.12.~~Poles near the real axis ' ... 'found by {\tt aaa}'],FS,11) %% % \vskip 1.02em %% % % This is a striking configuration. In fact, the function $y$ defined % by (21.9) is {\em doubly periodic} in the complex $t$-plane, % periodic not only in the real but also the imaginary direction. % (The poles of the mathematically exact function % keep going forever in a perfectly regular array, % but {\tt aaa} has just estimated some of % the poles near $[\kern .3pt 0 , 100\kern .3pt]$.) % It is the function known as the {\em Jacobi sine function} with % parameter $m = 0.998$, written $\hbox{sn}(t,0.998)$. % Although it is analytic for real values of $\kern .3pt t$, it has poles for % complex~$t$, and they lie in an infinite doubly periodic array. % Here is the period in the imaginary direction (close to $\pi$, but % different): % imagpol = sort(abs(imag(pol))); T2 = 2*min(abs(imag(pol))) %% % % \noindent % If we include horizontal lines in the plot too, we % get an image marking fundamental domains of periodicity. % for k = -10:10 plot([-10 110],k*T2*[1 1],LW,.7,CO,.75*[1 1 1]) end plot(pol,'.r',MS,10), hold off title(['Fig.~21.13.~~Lines mark both real and ' ... 'imaginary periods'],FS,11) %% % \vskip 1.02em %% % % \noindent % Here, arbitrarily, we verify the double periodicity % by comparing $r(t)$ for $t = 5-i$ and the same value plus % $T$, $iT_2$, and $T+iT_2$. (See Exercise 21.8.) % t = 8-1i; r(t+[0; T; 1i*T2; T+1i*T2]) %% % % Doubly periodic functions, also known as {\em elliptic functions}, are % among the core topics of complex analysis. % %% % % \smallskip % {\sc History.} The idea of giving physical meaning to imaginary % space or time has proved fruitful for quite a few % problems of engineering and physics, from the design of transonic airplane wings % to the explanation of the Big Bang. Stephen Hawking presented the latter idea % in his 1988 blockbuster book {\em A Brief History of Time}. The elegant, % singularity-free frame of reference in which to understand the % universe, Hawking suggests, % is one in which $\kern .3pt t\kern .3pt$ runs in the % imaginary direction. In that direction, he % proposes, the universe's boundary conditions are periodic, % so there are no boundaries and % no singularities. It's only if you turn % a right angle and do analytic continuation into the real-$t$ direction --- % analytically continuing % back to around $t = -\hbox{14,000,000,000}$ years, to be precise --- % that you encounter the mother of all singularities. % \smallskip % %% % % {\sc Our favorite reference.} There is an excellent classic % book by Einar Hille on {\em Ordinary Differential Equations in the Complex Domain,} % but our favorite reference for this topic, with spectacular figures % illustrating the functions in question, is Fornberg and % Weideman, A numerical methodology for the Painlev\'e equations,'' % {\em Journal of Computational Physics,} 2011. % %% % % \smallskip % \begin{displaymath} % \framebox[4.7in][c]{\parbox{4.5in}{\vspace{2pt}\sl % {\sc Summary of Chapter 21.} Many ODE\kern .5pt s % make sense for complex as well as real arguments. % Solutions\/ $y(t)$ for complex\/ $t$ can be defined by applying the % ODE in the complex plane, or by analytic continuation from % solutions for real\/ $t$. If $y(t)$ changes rapidly as\/ $t$ % varies through certain real values, there are usually one % or more singularities of\/ $y$ nearby in the complex $t$-plane. % Typically such singularities will be branch points, requiring % associated branch cuts, but in the case of Painlev\'e equations, % the only movable singularities (i.e., lying at data-dependent % locations) are poles. % \vspace{2pt}}} % \end{displaymath} % %% % \smallskip\small\parskip=1pt\parindent=0pt % {\em Exercise $21.1$. Painlev\'e I change of variables.} % Suppose $y(t)$ is a solution to (21.5). Show that another % solution is $\omega^3 \kern .3pt % y(\omega\kern .3pt t)$, where $\omega = \exp(2\pi i/5)$. % \par % {\em \underline{Exercise $21.2$}. Multiple solutions to Painlev\'e BVP.}\ \ Compute % and plot another solution to (21.7) by starting from the % initial guess $y(t) = -\exp(-10(t+1)^2)$. What is its minimum value? % \par % {\em \underline{Exercise $21.3$}. Some singularities are just pseudo-singularities.} % {\em (a)} Repeat the calculation of Figure 20.4 for the ODE of (20.6), % $\varepsilon y'' + x\kern .1pt y' + x\kern .1pt y = 0$ with $\varepsilon = 0.001$. % Then plot the poles of the AAA approximants of the solution $y$, as in Figures % 21.7 and 21.12, for ${\tt tol} = 10^{-4}$, $10^{-6}$, and $10^{-8}$. % You will see that the poles of these approximations vary with the % tolerance. They are not approaching any actual singularities of the % true solution $y$, for $y$, despite % its steep interior layer near $x=0$, cannot have any singularities since the % equation is linear with analytic coefficients and a nonzero % leading-order coefficient. % {\em (\kern .7pt b)} In such a case % $y$, though analytic as a function of $x$, must increase % rapidly in magnitude as $x$ leaves the real axis. (One can prove % this via {\em Cauchy's estimate.}) To confirm % this, examine values of $x$ along the % imaginary axis. First, make the substitution $x = i s$ % as in (21.3) and write down the form that the ODE takes when % transformed to the $s$ variable. % {\em (c)} Now solve (20.6) numerically % for $x \in [\kern .3pt 0, 0.2 \kern .3pt i \kern .2pt ]$, that is, % $s \in [\kern .3pt 0, 0.2 \kern .2pt ]$, % taking as initial data the appropriately % transformed values $y(0)$ and $y'(0)$ from % the solution of {\em (a)}. Make a semilogy plot of % $|y(x)|$ against $|x|$. How big is $|y(0.2\kern .3pt i)|\kern .5pt ?$ % \par % {\em \underline{Exercise $21.4$}. Bypassing blowup.} As we have discussed % at several points in the book, positive % solutions to $y' = y^2$ blow up to $\infty$ in finite time. The situation % changes, however, if $y$ is even very slightly complex. % {\em (a)} Draw a quiver plot of $y' = y^2$ in the % complex $y$-plane and describe the shape of the complex % trajectories in this plot. % {\em (\kern .7pt b)} Solve $y' = y^2 + f$ with $y(0)=1$ % for $t\in [\kern .3pt 0,1]$, where % $f$ is the complex random function generated by {\tt rng(1)}, % \verb| f = 0.01*randnfun([0 1],'norm','complex')|. % Plot the solution with {\tt axis equal}. (For % analysis of such problems see % Herzog and Mattingly, `Noise-induced stabilization of % planar flows I'', {\em Electronic Journal of Probability,} 2015.) % \par % {\em \underline{Exercise $21.5$}. Jacobi sine function.} % A chebfun {\tt y2} representing the Jacobi sine function % as plotted in Fig.~21.11 could be constructed directly % from the Matlab function \verb|ellipj(t,0.998)|. % Do this and plot {\tt abs(y-y2)} to give an indication of accuracy of % this Chebfun ODE solution. %