%% Problem Set D, Problem 1
%% Clear variables and figures

clear, close all

%
%% Part (a)
% Here we plot a numerical solution of Airy's equation (as a dashed line) and an exact solution of the facsimile equation (solid line) near t = 0. The two solutions agree well for t between -1 and 1, but they diverge rapidly outside this range.

airyeq = @(t, y) [y(2); t*y(1)];
[tfor, yfor] = ode45(airyeq, [0, 2], [0, 1]);
plot(tfor, yfor(:,1), '--')
hold on
[tbak, ybak] = ode45(airyeq, [0, -2], [0, 1]);
plot(tbak, ybak(:,1), '--')
syms t y(t); Dy(t) = diff(y(t), t);
facA = dsolve(diff(y(t), t, 2) == 0, y(0) == 0, Dy(0) == 1)
fplot(facA, [-2, 2])
hold off
title 'Numerical solution and facsimile near t = 0'

%
%% Part (b): numerical solution
% Here is a numerical solution to Airy's equation for t_ near -16 = -4^2. This time we use deval to extract the numerical solution (but not its derivative) at the desired values of t. We also compute the numerical solution out to t = -20 so that we can use it again later in part (d).

solb = ode45(airyeq, [0, -20], [0, 1]);
tt = -18:0.05:-14;
yy = deval(solb, tt, 1);
[min(yy) max(yy)]
plot(tt, yy, '--')
hold on

%
%% Part (b): facsimile solution
% After viewing the graph, we see that it does resemble a sine curve. However, it probably does not lie close to the solution of y'' = -16y with the same initial condition at t = 0 because the two differential equations are not close to each other near t = 0. The amplitude of the oscillation is, from the graph, about 0.6, but we can make this more precise by looking at the data. We printed above the minimum and maximum values of the numerical solution yy between t = -18 and -14, and judge the amplitude to be about 0.64. So our facsimile solution should be y = 0.64\sin(4t + c_2) for some c. Since the graph crosses zero in the increasing direction near t = -16.25, we choose 4(-16.25) + c_2 = 0, or c_2 = 65. The graphs agree quite well over the entire interval.

facB = 0.64*sin(4*t + 65);
fplot(facB, [-18, -14])
hold off
axis([-18 -14 -1 1])
title 'Numerical solution and facsimile near t = -16'

%
%% Part (c): numerical solution
% Now here is a numerical solution to Airy's equation for t near 16 = 4^2.

solf = ode45(airyeq, [0, 18], [0, 1]);
tt = 14:0.05:18;
yy = deval(solf, tt, 1);
plot(tt, yy, '--')
hold on
axis([14 18 0 4*10^21])

%
%% Part (c): facsimile solution
% The solution appears to be growing exponentially, which is also what the hyperbolic sine does far from zero. In fact, \sinh(t) = (\exp(t) - \exp(-t))/2, or approximately \exp(t)/2 for large t. So, the proposed facsimile solution is approximately c_1\exp(4t + c_2)/2 = (c_1\exp(c_2)/2)\exp(4t). This suggests that we can choose c_2 arbitrarily (say c_2 = 0) and choose c_1 so that the magnitudes of the solutions match. Below we choose c_1 to make them match at t = 18. Again the graphs agree well, though since the solutions are so much larger between t = 17 and 18 than on the rest of the interval, we can't see how well they agree over the whole interval from a standard graph. In a separate figure, we show the two solutions in a ``semilog'' graph, where the vertical axis has a logarithmic scale. On such a scale, the graph of an exponential function should be straight. We see that the two solutions do agree well over the entire interval.

c1 = 2*yy(end)/exp(4*18)
facC = @(t) c1*sinh(4*t);
fplot(facC, [14, 18])
hold off
axis([14 18 0 4*10^21])
title 'Numerical solution and facsimile near t = 16'
figure
semilogy(tt, yy, '--', tt, facC(tt), 'b')
title 'Semilog plot of solution and facsimile near t = 16'

%
%% Part (d)
% Finally, we show the numerical solution over the interval (-20,2). We use the numerical solutions solb and solf that we computed in parts (b) and (c). From the graph below, it appears that the frequency increases and the amplitude decreases as t goes to minus infinity. The increasing frequency is predicted by the analysis in part (b), where the approximate frequency K increases like \sqrt{-t}. The decreasing amplitude is not explained by the facsimile approach, at least not without some more detailed analysis.

ttb = -20:0.05:0; 
yyb = deval(solb, ttb, 1);
ttf = 0:0.05:2; 
yyf = deval(solf, ttf, 1);
plot(ttb, yyb, 'b', ttf, yyf, 'b')
axis([-20 2 -2 3])
title 'Oscillations in the solution to Airy''s equation'