Consider the initial value problem
We will use the 'dirfield' module to plot the direction field for this differential equation. As directed in the text, we will also look at the solution curves with = 0 and = 1. The upper curve corresponds to = 1 and the lower curve corresponds to = 0.
The MATLAB module 'dirfield.m' used to generate the figures below is available here.
We can "see" the solution with = 0 on the direction field above. It starts with the two solutions we plotted, then bends up in a somewhat symmetric rounded "V" shape as t crosses 0. Our job now is to pinpoint this value of 0 using Euler's method.
The hint in the text is to restrict 0 to an interval of width 0.1. If we use dirfield, we can see that if = 0.6, the solution curve converges, and if = 0.7 the solution curve diverges. We'll start with this information to try to further pin down the value of 0.
We will start this process with the figure above open, so we can draw our numerical solution curves to it. (Run 'dirfield' with y'=y2 - t2, then draw two solution curves with (t0, y0) = (0, 0.6) and (t0, y0) = (0, 0.7). Exit dirfield.) We will implement Euler's method in MATLAB, repeating it for various values of until we can get close to the critical value 0.
Recall the formula for Euler's method:
We will use the second form of Euler's method because it is easier to code into MATLAB and we don't have to be as picky about the indices. (We can't start with a zero entry in a matrix in MATLAB.)
>> clear; % Free memory -- this does not affect the figure
>> diffeq=inline('y^2-t^2','t','y'); % Define function
>> h=0.01; % Set stepsize
>> t0=0; % Set initial t-value
>> y0=0.6; % Set initial y-value
>> t=[0:h:5]; % Make an array of t's from 0 to 5 with stepsize h
>> num=size(t,2); % How many t's do we have?
>> y(1)=y0; % Our first data pt. is (t0,y0)
>> for i=2:num % Loop to implement Euler's method
y(i)=y(i-1)+feval(diffeq,t(i-1),y(i-1))*h; % y(i)=y(1-i)+h*f'(t(i-1),y(i-1))
>> plot(t,y,'k'); % Plot results in black on same figure
The graph above seems to be identical to the one from part (a) except that part of the lower curve is black. The black part is the numerically approximated solution we found using Euler's method. It's a really good approximation to the curve MATLAB found using it's ODE solver. Now we want to hunt for the value of 0. We don't need to repeat all of the above MATLAB commands, we just need to re-define y0, reset y(1), loop to calculate the y-values, and plot the new result. For example, if we reset y0 to .65:
>> for i=2:num
Below is a graph where this process has been repeated for various values of y0. The critical value 0 is between the values of y0 = 0.675 and y0 = 0.68, corresponding to the last curve that converges and the first one that diverges (bottom to top). The reader can repeat this process to further refine our approximation of 0.