%% 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.