%% Numerical and Graphical Methods for 1st Order Equations

% Thus far we have studied first order ordinary 
% differential equations for which we can find a symbolic or formula 
% solution. Years ago, that is largely all that was done 
% in a sophomore-level ODE course. But this masked the 
% fact that, indeed, most ODEs cannot be solved 
% symbolically. Just to illustrate that fact, consider 
% the simple equation

%           y' + y = t.

% This equation is solved symbolically rather easily. But 
% let us change it slightly by introducing a parameter. 
% Consider now the equation

%           y' + y^b = t,

% where b is a real number. In fact it would not take you 
% long to discover that Matlab can solve this equation 
% only in case b = 0, 1, or 2. For no other values of b, 
% positive or negative, can Matlab produce a formula 
% solution of this differential equation. Nevertheless, 
% situations arise in which we need to be able to say 
% something about solutions to equations like this. In 
% this course we develop two techniques for doing so:
% graphical and numerical. (Actually, as you have seen, 
% we also develop some qualitative techniques, but they 
% are much more ad hoc than the algebraic, graphical and 
% numerical techniques that you will encounter.)

% For the purposes of illustration, we shall consider 
% three examples: 
%       b = 1, 1/2 and -2. 

% The first example Matlab can solve of course, but not 
% the other two. For the record, we'll try them; Matlab 
% reports its failure in both instances.
clear all
close all
dsolve('Dy + y = t', 't')

%%
dsolve('Dy + sqrt(y) = t', 't')

%%
dsolve('Dy + 1/y^2 = t', 't')

%% Direction Fields
% You have learned the direction field method for first 
% order equations 
%       y' = f(t, y). 

% If one draws a small vector at the point (t_0, y_0) 
% having slope f(t_0, y_0), then one knows that the 
% solution curve to the differential equation that passes
% through (t_0, y_0) must have slope = f(t_0, y_0) there.
% If we draw lots of such small slope line segments, then
% we get a global picture of the solution curves to the
% differential equation, which we have called the
% Direction Field of the equation. Naturally, Matlab can 
% do this much better than we can freehand. So here are 
% the direction fields for the three examples in question:
clear all
close all

%%
figure; 
[T, Y] = meshgrid(0:0.2:5, -3:0.2:3);
S = -Y + T;
L = sqrt(1 + S.^2);
quiver(T, Y, 1./L, S./L, 0.5), axis ([0 5, -3 3])
title 'Direction Field for dy/dt + y = t'
set(get(gca, 'Title'),'Fontsize', 25)
set(gca,'XTickLabel',get(gca, 'XTickLabel'),'FontSize',25)
set(gcf, 'Position', [1 1 1920 1220])
hold on; T=0:0.1:4; plot(T, T-1, 'r')
figureforbis1 = gcf

% The last command will allow me to recall this figure 
% later. You can see the linear solution y = t - 1. (I
% included it drawn in red.) The other solutions are 
% obtained by adding a constant times exp(-t) to it.

%% 
% Now the first of the two equations for which we cannot 
% find a formula solution. Since we are taking the square 
% root of y, it behooves us to restrict y to positive 
% values.

figure; set(gcf, 'Position', [1 1 1920 1420])
[T, Y] = meshgrid(0:0.2:5, 0.2:0.2:3);
S = -sqrt(Y) + T;
L = sqrt(1 + S.^2);
quiver(T, Y, 1./L, S./L, 0.5), axis([0 5, 0 3])
title 'Direction Field for dy/dt + sqrt(y) = t'
set(get(gca, 'Title'),'Fontsize', 25)
set(gca,'XTickLabel',get(gca, 'XTickLabel'),'FontSize',25)
hold on; T=0:0.1:4; plot(T, T-1, 'k')

% Gee it looks sort of similar to the first case. You have 
% to look closely to see that the line y = t -1 is not 
% really a solution. (Again I graphed the line so that
% you can see that in this case it is not a solution 
% curve.) Let's refine the parameters a little 
% to see that even better.

%%
figure; set(gcf, 'Position', [1 1 1920 1420])
[T, Y] = meshgrid(0:0.2:6, 0.2:0.2:5);
S = -sqrt(Y) + T;
L = sqrt(1 + S.^2);
quiver(T, Y, 1./L, S./L, 0.5), axis([0 6, 0 5])
title 'Dir Field for dy/dt + sqrt(y) = t; a closer look'
set(get(gca, 'Title'),'Fontsize', 25)
set(gca,'XTickLabel',get(gca, 'XTickLabel'),'FontSize',25)
figureforbisonehalf = gcf

% The non-linearity is somewhat clearer now.
% Still, all the solutions seem to be squeezing together 
% as time increases. We'll see that more clearly when we 
% generate a numerical solution below.

%%
% Let's put the line in again just to see more clearly 
% that there is no linear solution.

figure; set(gcf, 'Position', [1 1 1920 1420])
[T, Y] = meshgrid(0:0.2:6, 0.2:0.2:5);
S = -sqrt(Y) + T;
L = sqrt(1 + S.^2);
quiver(T, Y, 1./L, S./L, 0.5), axis([0 6, 0 5])
title 'Dir Field for dy/dt + sqrt(y) = t; a closer look'
set(get(gca, 'Title'),'Fontsize', 25)
set(gca,'XTickLabel',get(gca, 'XTickLabel'),'FontSize',25)
hold on
T=0:0.1:6; plot(T, T-1, 'k')

%%
% Here's an exercise for you:
% Show that 
%       y = (1/4)t^2

% solves the equation. That is the limiting behavior of 
% all the solution curves. But Matlab cannot figure that 
% out on its own. However,
plot(T, 0.25.*T.^2)

%%
% And now let's look at the third case on which Matlab's 
% symbolic solver also chokes. Since y appears in 
% the denominator of the equation, we had better avoid 
% y = 0.

figure; set(gcf, 'Position', [1 1 1920 1420])
[T, Y] = meshgrid(0:0.2:5, -3:0.4:3);
S = -1./Y.^2 + T;
L = sqrt(1 + S.^2);
quiver(T, Y, 1./L, S./L, 0.5), axis([0 5, -3 3])
title 'Direction Field for dy/dt + 1/y^2 = t'
set(get(gca, 'Title'),'Fontsize', 25)
set(gca,'XTickLabel',get(gca, 'XTickLabel'),'FontSize',25)

% Well we can clearly see the problem across the line 
% y = 0. Perhaps I have been too clever by choosing the 
% step size = 0.4 so that I would miss y = 0 in the 
% computations. But there clearly is a problem. Above the
% t-axis, the solution curves come down then seem to go 
% up off to infinity. Below the t-axis, the solution 
% curves also come down initially, but then they seem to 
% be crashing into the t-axis -- but perhaps they are 
% just tending toward it asymptotically; it's hard to 
% tell. We'll come back to this when we compute numerical
% solutions.

%% Numerical Solutions

% Okay, it's time to get going on numerical methods for 
% these equations. The first can be handled symbolically, 
% as we know, but we'll solve it numerically just to 
% establish the paradigm. Of course to work numerically
% with ode45, we must solve an initial value problem and 
% so we are confronted with a choice of initial value of 
% the time paramter t. I'll choose t = 0.

f = @(t, y) -y + t

% The following commands produce solution curves on the 
% interval [0, 5] for y(0) = -3, -2.8, ... 2.8, 3.
[t, y] = ode45(f, [0 5], -3:0.2:3);
figure
plot(t, y)

% Just what we expect, since the general solution was 
%       t - 1 + c*exp(-t). 
% Let's overlay the direction field.

%%
figure(figureforbis1)
hold on 
plot(t, y)
title 'dir field and soln curves for dy/dt + y = t'
hold off

% Let's also note that if we solve the equation for y' 
% and take the partial wrt to y, we get -1, which is < 0.
% So according to Theorem 5.2 in Section 5.3 of DEwM, 
% the equation is stable everywhere as the graphic 
% suggests.

%%
% Now let's move on to the first of the two equations 
% without symbolic solutions.  We'll choose the parameters 
% to match those that gave us the best direction field 
% earlier in the session. Specifically, we'll restrict
% the t-axis to [0, 6] and the initial y-values to vary 
% between 0.2 and 5, and we'll overlay the direction field. 

% Here's the code:
g = @(t, y) -sqrt(y) + t;
figure(figureforbisonehalf)
[t, y] = ode45(g, [0 6], 0.2:0.2:5);
hold on
plot(t, y)
title 'dir field and soln curves for dy/dt + sqrt(y) = t'

% Again, if we solve for y' and take the partial wrt y, 
% we obtain -(1/2)(1/sqrt(y)), which for positive y, 
% is again < 0, so we have stablity as the graphic 
% suggests.

%%
% By the way, I suggested that the asymptotic curve for 
% all the solutions was y = (1/4)t^2. Let's superimpose 
% that curve on the graph.
Z = 0:0.2:5;
plot(Z, (1./4).*Z.^2, 'k','LineWidth', 2)

%%
% Finally, let's look at the most intractable case and see 
% whether we can clear up the mystery that we encountered 
% when we drew the direction field. Because of the term 
% 1/y^2 in the differential equation, we know that there 
% are problems on the t-axis. Let's break the problem 
% into two domains, initial data above the t-axis and 
% then below the t-axis. In each case, I'll first redraw 
% the direction field, then plot some numerical solution 
% curves and finally combine them. Before we do so, let's
% apply the partial wrt y test (i.e., Theorem 5.2 from 
% Section 5.3 of DEwM. The partial is 2y^(-3), which 
% suggests stability below the t-axis and instability 
% above. Let's handle the stable case first:

%%
figure; set(gcf, 'Position', [1 1 1920 1420])
[T, Y] = meshgrid(0:0.2:5, -3:0.2:-0.2);
S = -1./Y.^2 + T;
L = sqrt(1 + S.^2);
quiver(T, Y, 1./L, S./L, 0.5), axis([0 5, -3 0])
title 'Dir Field for dy/dt + 1/y^2 = t below the t-axis'
set(get(gca, 'Title'),'Fontsize', 25)
set(gca,'XTickLabel',get(gca, 'XTickLabel'),'FontSize',25)
figureforbisminustwobottom = gcf

% Not so clear what is going on. Best guess is that 
% solution curves that start far down on the y-axis tend 
% toward the t-axis asymptotically; and those that start 
% on the y-axis up near zero first decline and then move 
% back up toward the t-axis asymptotically also. There 
% also seem to be solution curves that come off the 
% t-axis vertically down and then tend back toward it 
% asymptotically. We'll verify some of this by first 
% looking at a single curve, then a family and then by 
% superimposing the direction field.

%%
h = @(t, y) -1./y.^2 + t
figure; set(gcf, 'Position', [1 1 1920 1420])
ode45(h, [0 25], -0.02)

% The sharp turn is consistent with the 
% direction field. But it's not clear that the asymptote 
% is the t-axis. Can you explain from the differential 
% equation why that must be so? 
%           y' + 1/y^2 = t

% Now for more
% solution curves superimposed on the direction field.

%%
figure(figureforbisminustwobottom)
[t, y] = ode45(h, [0 6], -3:0.2:-0.2);
hold on
plot(t, y)
title 'dir field and soln curves for dy/dt + 1/y^2 = t'
hold off

%%  Observations
% I want to examine this picture a little more closely 
% and make three observations:

% 1. The solution curves appear to be of two distinct 
% types: those that start on the y-axis, perhaps decrease
% slightly for a bit and then increase slowly toward the 
% t-axis; and those that seem to start on the t-axis, 
% come off it vertically down, and then also increase 
% slowly back to the t-axis. Moreover, there seems to be 
% a solution curve that separates the two behaviors. This
% is very reminiscent of the idea of a "separatrix," 
% which we shall study later in the course.

% 2. For the former of the two types, the domain of 
% definition of the solution curve appears to be the 
% positive t-axis. But for the latter type, the solution 
% curve is defined on an interval of the form 
% (alpha, infty), for some positive number alpha. 
% Those curves come off the t-axis with infinite 
% negative slope.

% Question: Why for large t must the separatrix in fact 
% be asymptotic to
%       y = -1/sqrt(t)?

% Answer: Stare at the diff eqn:
%        y' + (1/y^2) = t.
% y' is tending to zero, so 1/y^2 is asymptotic to t for 
% large t.

% 3. All the solution curves are coming together as 
% t->infty; which is what we expect since below the 
% t-axis, we know from the partial test, that we have 
% stability.
%%
% Our last chore is to examine the situation above the 
% t-axis where the equation is unstable.

figure; set(gcf, 'Position', [1 1 1920 1420])
[T, Y] = meshgrid(0:0.2:5, 0.2:0.2:3);
S = -1./Y.^2 + T;
L = sqrt(1 + S.^2);
quiver(T, Y, 1./L, S./L, 0.5), axis([0 5, 0 3])
title 'Dir Field for dy/dt + 1/y^2 = t above the t-axis'
set(get(gca, 'Title'),'Fontsize', 25)
set(gca,'XTickLabel',get(gca, 'XTickLabel'),'FontSize',25)
figureforbisminustwotop = gcf

% Looks like the reverse of the previous situation where
% instead of the arrows coming together they are flying 
% apart. Again, let's look at a typical problematic curve 
% and then several curves.

%%
figure; set(gcf, 'Position', [1 1 1920 1420])
ode45(h, [0 20], 0.02)

%%
% You really see the instability there; the numerical 
% solution jumped across the t-axis and started traveling 
% on a solution curve below it. Here's an example where 
% you can see it even better
figure; set(gcf, 'Position', [1 1 1920 1420])
ode45(h, [0 20], 0.75)

%% 
% If we can stay far enough above the t-axis at our 
% initial value, we can get a more reasonable graph.
figure(figureforbisminustwotop)
set(gcf, 'Position', [1 1 1920 1420])
[t, y] = ode45(h, [0 6], 1.5:0.1:3);
hold on
plot(t, y)
title 'dir field and soln curves for dy/dt + 1/y^2 = t'

%%
% Observations that are very similar to those above hold 
% here:

% 1. Above a certain point the solution curves decrease 
% slightly, then go flying off to infinity; but below 
% that point, they crash into the t-axis. It's a little 
% harder to see, but once again, there is a separatrix 
% curve that separates the two behaviors. 
% Reasoning as above, you should be able to see that for 
% large t, that curve is asymptotic to 
%       y = 1/sqrt(t).

%% In fact here is graphical evidence
W=0.2:0.1:5;
plot(W,1./sqrt(W), 'k','LineWidth', 2)

%%
% And here is the real separatrix:
[t, y] = ode45(h, [4 -2], 0.5);
plot(t,y, 'gr', 'Linewidth', 4)
%Can you explain what I did?

%% 2. The former solution curves are defined over the 
% positive t-axis, but the latter curves are defined on 
% an interval of the form (0, beta), for some positive 
% number beta.

% 3. The diagram illustrates the instability of the 
% equation above the t-axis.

%% Final Thoughts

% Try to visualize the two portraits (above and below 
% the t-axis) together. They fit together with vertical 
% tangents to the solution curves at the t-axis. Try 
% printing the two portraits and hold them up together.

% Exercise:
% See if you can figure out what the solution curves look
% like for negative t and try to do some analysis of the 
% curves analogous to what I did here for positve t.