Section 6.1 - Euler's Method

Contents

Why is it often necessary to use graphical and numerical techniques?

Here is an example of an ODE which MATLAB claims to solve explicitly, but the solution is very complicated. The initial value problem is

$$y'=y^2+t^2, \quad y(0)=1.$$

Here is MATLAB's "exact solution."

y = dsolve('Dy = y^2 + t^2, y(0)=1', 't'); pretty(y)

Warning: Explicit solution could not be found; implicit solution returned. 
Warning: The solutions are parameterized by the symbols:
z = solve([limit((t*gamma(3/4)*besselK(-3/4, (t^2*1i)/2)*4i -
2^(1/2)*4^(1/4)*C3*t*PI^(3/2)*besselI(-3/4,
(t^2*1i)/2)*1i)/(4*gamma(3/4)*besselK(1/4, (t^2*1i)/2) +
2^(1/2)*4^(1/4)*C3*PI^(3/2)*besselI(1/4, (t^2*1i)/2)), t, 0) = 1], C3,
NoWarning, VectorFormat)
 To include parameters and conditions in the solution, specify the
 'ReturnConditions' option. 
/        / 3 \        /   3     \               1/4       3/2
| t gamma| - | besselk| - -, #1 | 4i - sqrt(2) 4    t z pi
\        \ 4 /        \   4     /

          /   3     \    \ /        / 3 \        / 1     \            1/4
   besseli| - -, #1 | 1i |/| 4 gamma| - | besselk| -, #1 | + sqrt(2) 4    z
          \   4     /    / \        \ 4 /        \ 4     /

     3/2        / 1     \ \
   pi    besseli| -, #1 | |
                \ 4     / /

where

          2
         t  1i
   #1 == -----
           2

The solution is very complicated and involves $z$, which is the solution of a very complicated equation. Also, an attempt to plot it, even on an interval where it should be possible, gives a warning and an empty plot:

ezplot(y,[0,0.8])
title('')

Warning: Contour not rendered for non-finite ZData

There is a discussion of this example and how to solve it using MuPad on pp;. 72-74 of Differential Equations with MATLAB.

We can solve it exactly using Maple, with the command:

dsolve({diff(y(t), t) = y(t)^2+t^2, y(0) = 1}, y(t))

The solution is:

syms t
y = -t*(-(gamma(3/4)^2-pi)/gamma(3/4)^2*besselj(-3/4,1/2*t^2)+...
    bessely(-3/4,1/2*t^2))/(-(gamma(3/4)^2-pi)/gamma(3/4)^2*besselj(1/4,1/2*t^2)...
    +bessely(1/4,1/2*t^2));

To get a better idea of what the solution is like, we can view it graphically. One method is to plot the solution.

ezplot(y,[0,0.8])
xlabel 't', ylabel 'y'
title 'Solution of dy/dt = y^2 + t^2'

Another method is to look at the direction field using the quiver command.

[T,Y] = meshgrid(-3:0.2:3,-3:0.2:3);
S = T.^2 + Y.^2;
L = sqrt(1 + S.^2);
quiver(T, Y, 1./L, S./L, 0.5), axis tight
xlabel 't', ylabel 'y'
title 'Direction field for dy/dt = y^2 + t^2'

A third method is to solve numerically and plot the resulting numerical solution. Here we do that using ode45.

f = @(t,y) y^2 + t^2;
ode45(f,[0,0.8],1)

Note that this plot looks the same as what we found when we plotted the explicit solution.

Here's another example where MATLAB finds a complicated solution.

$$y' = y^2 + t, \quad y(0)=0.$$

y = dsolve('Dy = y^2 + t, y(0) = 0','t')

 
y =
 
(airy(3, -t) + 3^(1/2)*airy(1, -t))/(airy(2, -t) + 3^(1/2)*airy(0, -t))

Note that this solution has a different appearance from the one discussed on pp. 62-63. Starting with R2012a, the solution is in terms of MATLAB's notation for the Airy functions instead of MuPAD's notation.

Here's the direction field.

[T,Y] = meshgrid(-3:0.2:3,-3:0.2:3);
S = T + Y.^2;
L = sqrt(1 + S.^2);
quiver(T, Y, 1./L, S./L, 0.5), axis tight
xlabel 't', ylabel 'y'
title 'Direction field for dy/dt = y^2 + t'

Here is a plot of the solution using ezplot.

ezplot(y,[0 2])
title('Exact solution')

Here's a plot using plot.

tvals = 0:0.01:2;
plot(tvals,double(subs(y,'t',tvals)))
xlabel 't', ylabel 'y'
title 'Exact solution'

Can this be correct?

Alternatively we can plot the numerical solution obtained with ode45.

f = @(t,y) y^2 + t;
[t,ya] = ode45(f,[0,2],0);
plot(t,ya)

Warning: Failure at t=1.986314e+00.  Unable to meet integration tolerances
without reducing the step size below the smallest value allowed (3.552714e-15)
at time t.

Can this be correct?

Here is a problem MATLAB can't solve exactly.

$$y' = \mbox{--}\,t^2 y^3 + \cos(t), \quad y(0) = 0.$$

There is no choice but to use graphical and numerical techniques to analyze this equation. Here is the direction field.

[T,Y] = meshgrid(-3:0.2:3,-3:0.2:3);
S = -T.^2*Y.^3  + cos(T);
L = sqrt(1 + S.^2);
quiver(T, Y, 1./L, S./L, 0.5), axis tight
xlabel 't', ylabel 'y'
title 'Direction field for dy/dt = -t^2 y^3 + cos(t)'

Here is a plot of the numerical solution obtained with ode45.

f = @(t,y) -t^2*y^3 + cos(t);
[t,ya] = ode45(f,[0:100],0);
plot(t,ya)

From this we see that the solution appears to tend slowly towards 0 in an oscillating manner as x gets very large.


Published with MATLAB® R2015b