%% Second Order Linear Equations and the Airy Functions:
%   Why Special Functions are Really No More Complicated 
%   than Most Elementary Functions

% We shall consider here the most important second order 
% ordinary differential equations, namely linear equations.
% The standard format for such an equation is

%           y''(t) + p(t) y'(t) + q(t) y(t) = g(t),

% where y(t) is the unknown function satisfying the 
% equation and p, q and g are given functions, all 
% continuous on some specified interval. We say that the 
% equation is homogeneous if g = 0. Thus:

%           y''(t) + p(t) y'(t) + q(t) y(t) = 0.

% We have studied methods for solving an inhomogeneous 
% equation for which you have already solved the 
% corresponding homogeneous equation. So we shall 
% concentrate on homogeneous equations here. The equation 
% is said to be a constant coefficient equation if the 
% functions p and q are constant. Again we have studied 
% methods for dealing with those. Non-constant coefficient
% equations are more problematic, but alas, they arise 
% frequently in nature. In this lesson we shall study 
% closely one of the best known examples -- namely, 
% Airy's Equation. For the record, the solutions to that 
% equation, i.e., the Airy functions, arise in 
% diffraction problems in the study of optics,
% and also in relation to the famous Schroedinger equation
% in quantium mechanics.

% Before proceeding, let's recall some basic facts about 
% the set of solutions to a linear, homogeneous second 
% order differential equation. The most basic fact is that
% the set of solutions forms a two-dimensional vector 
% space. This means that you can find two solutions, 
% y_1 and y_2, neither of which is a multiple of the 
% other, so that all solutions are given by linear 
% combinations of these two: 

%      ay_1 + by_2, where a and b are arbitrary constants.

% Just to help you get your bearings, let's mention two 
% other two-dimensional vector spaces that you know very 
% well: 

% The Euclidean plane, where every vector can be 
% expressed as
%       a i-> + b j->
% where i-> = (1,0) and j-> = (0,1);

% or all polynomials of degree at most 1
%       a + b x.

% Now returning to second order linear homogenous 
% differential equations with constant coefficients, we 
% note, by way of examples, that all solutions of 

%       y'' + y = 0 

% are given by 

%       a cos(t) + b sin(t); 

% and all solutions of 

%       y'' +(1/t)y' + (1/t^2)y = 0

% are given by 

%       a cos(ln(t)) + b sin(ln(t)). 

% (The latter is an Euler equation.) 

% In the first example, the interval of 
% definition is the whole real line, whereas in the second, 
% we must restrict to t>0. In both instances, 
% a and b are arbitrary real numbers. 

%% Airy's Equation
% This is the equation:

%       y'' - ty = 0. 

% In this presentation we shall solve it symbolically and 
% numerically. We shall also address it graphically. 
% Furthermore, if you look at DEwM, Problem 1, p. 157, you 
% will see a qualitative method of dealing with the 
% equation, which I shall briefly recall below. Finally, 
% although it is not in the course syllabus, one can also 
% use the method of series solution to solve Airy's 
% Equation. So there are lots of ways to skin this cat.

% Well, let's try dsolve first and see what happens. 
clear all
close all

%%
dsolve('D2y - t*y =0', 't')

% How about that, Matlab solves it. But what are those 
% functions it reports as the answers? In fact, they are 
% special functions, and if I may quote from DEwM, p. 52: 

% "By special functions we mean various non-elementary 
% functions that mathematicians give names to, often
% because they arise as solutions of particularly important 
% differential equations." Recall also that elementary 
% functions are "the standard functions of calculus: 
% polynomials, exponentials and logarithms,
% trigonometric functions and their inverses, and all 
% combinations of these functions through algebraic 
% operations and compositions." The simplest
% special function is the error function 

%   erf(t) = (2/sqrt(pi))int_0^t exp(-s^2) ds.

% Many special functions have integral formulas like the 
% above and/or are specified by a differential equation 
% and/or are given by a power series.

%% Graphing Airy functions
% In order to do so, we must recognize that there is still 
% some issue in the reporting of symbolic information in 
% Matlab as a consequence of Matlab having "bought" its 
% symbolic solver, i.e., MuPAD. We can see this by typing
help airy

%%
% Note that this is one of those cases wherein there is 
% some translation issue. In any event, the function that 
% is reported as airy(0,t) or more simply airy(t), is 
% just the normal Airy function Ai(t) and the other Airy 
% function reported as airy(2,t) is the Airy function
% Bi(t). More on these below. But first let's graph these 
% functions.

figure; 
ezplot('airy(t)')
set(get(gca, 'Title'),'Fontsize', 15)
set(gca,'XTickLabel',get(gca, 'XTickLabel'),'FontSize',15)
% set(gcf, 'Position', [1 1 1320 1020])
%%
% and on a bigger interval
figure; % set(gcf, 'Position', [1 1 1920 1420])
ezplot('airy(t)', [-10 10])
set(get(gca, 'Title'),'Fontsize', 15)
set(gca,'XTickLabel',get(gca, 'XTickLabel'),'FontSize',15)

% Looks like a mildly damped oscillation to the left, but
% with increasing frequency; and
% a rapid decay to zero to the right.

%%
% And the second solution
figure; % set(gcf, 'Position', [1 1 1920 1420])
ezplot('airy(2,t)')
set(get(gca, 'Title'),'Fontsize', 15)
set(gca,'XTickLabel',get(gca, 'XTickLabel'),'FontSize',15)

%%
% This time I will readjust the interval in order to see 
% the oscillation
figure; % set(gcf, 'Position', [1 1 1920 1420])
ezplot('airy(2,t)', [-10 2])
set(get(gca, 'Title'),'Fontsize', 15)
set(gca,'XTickLabel',get(gca, 'XTickLabel'),'FontSize',15)

% So both solutions manifest decreasing oscillations 
% toward minus infinity, apparently with increasing 
% frequencies, but as t -> +infinity, one solution 
% grows without bound and the other decays to zero.

%% Some Qualitative Analysis
% Now I recall Problem 1 in PSD in DEwM. It leads the 
% reader through the following reasoning. For a large 
% negative value of t, say t = -K^2, Airy's Equation 
% resembles

%       y'' + K^2 y = 0, 

% whose solutions are sinusoidal oscillations with 
% frequency K. Thus the oscillatory behavior on the 
% negative axis of the Airy functions is not surprising. 
% Moreover, as t moves toward -infty, and so has
% larger absolute value (K is getting bigger), then 
% the frequency is increasing. This analysis does not
% allow us to conclude anything about the amplitude. 

% On the other hand, for t large postive, say near K^2, 
% the equation resembles 

%       y'' - K^2 y = 0, 

% whose solutions are e^(Kt) and e^(-Kt). Thus the growth 
% at positive infinity is expected; and exactly as the 
% functions

%       a e^t + b e^(-t)

% have exactly "one direction" in which the function 
% decays at +infinity, whereas in all other directions 
% the solutions grow quickly, it is again not surprising 
% that the same behavior is manifested by the Airy 
% functions.

% Indeed, the basic Airy function airy(t) = airy(0,t) is
% exactly that special choice of the Airy functions.

%% Numerical solutions to yield a graphical presentation 

% Now we imitate the code on p. 142 of DEwM. As we saw 
% above, there are two arbitrary constants to be specified
% in the choice of an Airy function. That corresponds to 
% the fact that the second order Airy equation requires 
% two pieces of initial data to determine a specific 
% solution. Thus drawing a representative set of solutions 
% does not, like in the case of first order equations, 
% yield a set of non-intersecting curves. Still, sometimes 
% the pictures are striking and reveal the general nature 
% of solutions rather dramatically -- see e.g., the graph 
% on p. 142 of DEwM.

rhs = @(t, y) [y(2); t*y(1)];
figure; % set(gcf, 'Position', [1 1 1920 1420])
hold on
for y0 = -2:2
    for yp0 = -1:0.5:1
        [tfor, yfor] = ode45(rhs, [0 2], [y0, yp0]);
        [tbak, ybak] = ode45(rhs, [0 -2], [y0, yp0]);
        plot(tfor, yfor(:,1))
        plot(tbak, ybak(:,1))
    end
end
title('Airy Functions with Varying Initial Data', 'FontSize', 15)
set(gca,'XTickLabel',get(gca, 'XTickLabel'),'FontSize',15)

%% Stability 
% We have not considered stability of second order 
% equations, but it is not hard to envision what we would 
% mean -- small perturbations in the initial data -- both 
% position and velocity -- should lead to only small 
% perturbations in the solution curve over the long term. 
% Given what we have learned about the Airy functions, do 
% you think the Airy equation is stable?

% Not likely; just the form of the equation y'' = ty and 
% the derivative test suggests not. Let's see if we can 
% justify that assertion. We know that airy(t) is the only 
% solution of Airy's equation that decays at infinity. 
% Let's solve the equation numerically and compare what 
% we get (graphically) to the curve of airy(t):

figure; % set(gcf, 'Position', [1 1 1920 1420])
hold on
ezplot('airy(t)', [0 10])
yy0 = airy(0);
yyp0 = airy(1, 0);
[tt, yy] = ode45(rhs, [0 10], [yy0 yyp0]);
plot(tt, yy(:,1), 'r')
title('Symbolic Airy vs Numerical Approximation','FontSize',15)
set(gca,'XTickLabel',get(gca, 'XTickLabel'),'FontSize',15)

% Can you explain the graph?

%%
% This might help a little.
airy(0)
airy(1,0)
axis([0 10 -0.1 0.4])

%% The Nature of Special Functions
% Finally, I would like to convince you that special 
% functions like the Airy Function are really no more 
% mysterious than many of our elementary 
% functions -- like the sine, exponential or logarithm.

% In fact, how do you define sin(x)? In most caclulus 
% books, the function sin(x) is defined by saying: let x 
% measure an angle in radians, then draw a right triangle 
% with that angle, whereupon sin(x) is the ratio of the 
% length of the opposite side over the hypotenuse. If I 
% ask you to tell me what sin(sqrt(3)) is, do you really 
% have a good feeling for that value? So, many advanced 
% calculus books attempt to put the definition of the 
% sine function on a firmer analytic footing. They define 
% the sine function as follows. First define the function 
% ArcSin(t) by the integral formula:

%   ArcSin(t) = int_0^t 1/sqrt(1-s^2)) ds, -1 < t < 1.

% Then they engage in some calculus to establish that this 
% is a differentiable function, monotone increasing, 
%   InvSin(-1) = -pi/2,     InvSin(1) = pi/2, 
% and the tangent line is vertical at the two endpoints.
% Finally, they define sin(x) to be the inverse function of 
% ArcSin(t), which is then defined on [-pi/2, pi/2] and 
% finally they extend it to the whole real line by invoking 
% periodicity. 

% Now is that any simpler than the method we have used to 
% define airy(t)? I won't go through the derivation
% but I will tell you that the basic Airy function can also 
% be obtained by an integral. Here is the formula:

%   airy(t) = (1/pi) int_0^infty cos((1/3)s^3 + ts) ds.

% Not an easy integral, but we can deal with it numerically 
% when necessary. This process is perhaps a little more 
% complicated technically than:

%   ln(t) = int_0^t 1/s ds, 
%   exp(t) = inverse function of log(t);

% but aren't we talking about more or less the same kind 
% of object.

% The moral of the story: we can deal with airy(t) and most 
% special functions in the same way that we deal with 
% elementary functions: 

% graphically, numerically, analytically and even 
% symbolically on occasion. 
% You should not be intimidated by special functions.

% There are more special functions than you can imagine: 
% Bessel functions, Legendre functions, hypergeometric 
% functions, the Riemann Zeta function and a score more. 
% Many of these, but not all, are solutions of 
% second order homogeneous non-constant coefficient 
% equations. You may encounter some of them in your 
% physics or engineering courses. Let's just see if Matlab
% can deal with the Bessel functions. Those are the 
% solutions of the equation:

%       t^2 y''(t) + t y'(t) + (t^2 - n^2) y(t) = 0

% for different choices of an integer n.

dsolve('t^2*D2y + t*Dy +(t^2-n^2)*y=0')

%%
% I leave it to you to explore other special functions in 
% Matlab, but let us just draw the fundamental solutions 
% of the Bessel equation for n = 0.
figure; % set(gcf, 'Position', [1 1 1920 1420])
ezplot('besselj(0,t)')
title('Bessel function of first kind','FontSize', 15)
set(gca,'XTickLabel',get(gca, 'XTickLabel'),'FontSize',15)
%%
% Looks like a sinusoidal; let's redraw on a bigger 
% interval.

ezplot('besselj(0,t)', [-10 10])

%%
% and even bigger interval

ezplot('besselj(0,t)', [-100 100 -0.5 1])

%%
% Damped oscillation clearly. Actually, the oscillatory 
% behavior is not surprising. For t large, the middle 
% term in the equation 

%   t y'' + y' + t y =0

% is negligible and so the sinusoidal behavior is evident. 
% Also, it looks like the solution curve has the same 
% behavior in both directions. In fact bessel functions 
% of order zero are even. That's a good exercise: show 
% that if y(t) satisfies the Bessel equation (with n  = 0),
% then so does y(-t).

% Now for the other solution:
figure; % set(gcf, 'Position', [1 1 1920 1420])
ezplot('bessely(0,t)')
title('Bessel function of second kind','FontSize', 15)
set(gca,'XTickLabel',get(gca, 'XTickLabel'),'FontSize',15)
%%
% that's interesting; why no negative values?

ezplot('bessely(0,t)', [-5 5])

%%
% There are no negative values because bessely(0,t) has a 
% logarithmic singularity at t = 0. So Matlab discards the 
% negative values, which are just a mirror reflection of 
% the positive values. Let's see the singularity more 
% clearly.
figure; ezplot('bessely(0,t)', [0 0.01])
figure; ezplot('bessely(0,t)', [0 0.00001])

%%
% In fact lim bessely(0,t) as t-> 0+ is -Inf;
% although the divergence is very slow (logarithmic).

%%
% Let's look at one more -- the hypergeometric equation:

dsolve('t*(1 - t)*D2y + (c-(1 + a + b)*t)*Dy - a*b*y=0')
pretty(ans)
% solves the hypergeometric equation. I leave further 
% experimentation with special functions as solutions of 
% second order homogeneous equations to you; there are 
% lots of examples in Boyce & DiPrima. In fact, most of 
% the examples appear in the chapter on series solutions. 
% Since that topic is not covered in our syllabus, I will 
% have to leave it to you to look at on your own or to 
% encounter in other math or engineering courses.