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