First order differential equations: direction field, numerical solution

clearvars;

You must first download dirfield.m and put this file in the same folder as your mlx-files.

We consider a first order differential equation of the form

Example:

We will usually write this as

Defining the function :

Here the function

is given by

In Matlab we define this function by

f = @(t,y) t^2-y^2;

Note: we must use f=@(t,y)... even if t or y don't occur in the function.

Drawing the direction field

Consider a solution

: At

it has the value

, and it must have the slope

at this point.

We evaluate

on a grid of points in the

plane and indicate this slope by a short line:

For the grid we use the t-values -2,-1.8,...,1.8,2 and the y-values 0, 0.2,...,2.8,3:

dirfield(f,-2:.2:3,0:.2:3)

xlabel('t'); ylabel('y'); title('direction field for f(t,y)=t-y^2')

Note: Try different values for the step size (here 0.2) to obtain a nice looking direction field.

A solution of the ODE corresponds to a function where the graph "fits in the direction field".

Initial value problem

In order to get a unique solution we must specify an initial condition

Here we use the initial condition

with

and

t0=0; y0=1;
hold on;          % plot over the previous graph
plot(t0,y0,'bo')  % mark initial point (t0,y0) with 'o'

Finding an (approximate) solution with ode45/ode15s

Solution to the right

We want to find the solution with the inital t-value

and the final t-value

Note: we can also have

, this means we are looking for the solution to the left of the starting point.

T = 3;

We use the Matlab function ode45:

[ts,ys] = ode45(f,[t0,T],y0);

This gives a table of values for the solution function

[ts,ys]          % In Matlab you can click on the array to see all values.
ans = 45×2
                         0                         1
        0.0502377286301916         0.952203563527682
         0.100475457260383         0.909018631633333
         0.150713185890575         0.870095720468908
         0.200950914520766         0.835161527197856
         0.275950914520766         0.789997927007707
         0.350950914520766         0.752837036323012
         0.425950914520766         0.723408045488948
         0.500950914520766         0.701553297478841
         0.575950914520766         0.687225634210634
disp([ts,ys])    % this shows the whole array of values (scroll to view all rows)
                         0                         1
        0.0502377286301916         0.952203563527682
         0.100475457260383         0.909018631633333
         0.150713185890575         0.870095720468908
         0.200950914520766         0.835161527197856
         0.275950914520766         0.789997927007707
         0.350950914520766         0.752837036323012
         0.425950914520766         0.723408045488948
         0.500950914520766         0.701553297478841
         0.575950914520766         0.687225634210634
         0.650950914520766          0.68047421046911
         0.725950914520767         0.681358732690797
         0.800950914520766         0.689922946901339
         0.875950914520766         0.706192500523763
         0.950950914520767         0.730183932029128
          1.02595091452077         0.761833262306219
          1.10095091452077         0.800967057042828
          1.17595091452077         0.847311555009548
          1.25095091452077         0.900544579165089
          1.32595091452077         0.960232171878942
          1.40095091452077          1.02579929792095
          1.47595091452077          1.09657111798456
          1.55095091452077          1.17191229352202
          1.62595091452077          1.25113906579796
          1.70095091452077          1.33348478748553
          1.77595091452077          1.41818085277629
          1.85095091452077          1.50464010680906
          1.92595091452077          1.59230514721482
          2.00095091452077          1.68061755159882
          2.07595091452077          1.76910740729906
          2.15095091452077          1.85748133322935
          2.22595091452077          1.94548824252456
          2.30095091452077           2.0329347858997
          2.37595091452077           2.1197224934657
          2.45095091452077          2.20576140146816
          2.52595091452077          2.29099595268256
          2.60095091452077          2.37547956278226
          2.67595091452077          2.45932156354618
          2.75095091452077            2.542463142599
          2.82595091452077          2.62489783924626
          2.90095091452077          2.70676667304701
          2.92571318589057          2.73370293993516
          2.95047545726038           2.7605841040507
          2.97523772863019          2.78741200514159
                         3          2.81418865015293

The first column contains t-values starting with

and ending with T, the second column contains the corresponding y-values. Use disp to show the whole array (you can scroll to view all values).

We can plot the solution using

plot(ts,ys,'b')

Solution to the left

We now try to find a solution for t between

and

[ts,ys] = ode45(f,[t0,-2],y0);
Warning: Failure at t=-1.037465e+00.  Unable to meet integration tolerances without reducing the step size below the smallest value allowed (3.552714e-15) at time t.

The Matlab warning tells us that it was not able to find the solution all the way to

, it was only able to go up to

t1 = ts(end) % last entry of vector ts

t1 =          -1.03746487541147

plot(ts,ys,'b')

plot(t1*[1 1],ylim,'k--') % plot vertical line at t=t1 (dashed black)

hold off

title('Solution of IVP y''=t-y^2, y(0)=1')

Here we show the asymptote at

as a dashed vertical line.

The solution of this IVP exists on the interval

with

NOTE: If the solution ends at a point with vertical tangent is is better to use ode15s in place of ode45. In this situation the function ode45 may take a long time and may not recognize that the solution ceases to exist at the point with vertical tangent.