Numerical Solutions of Differential Equations

Enter the following command to load the necessary package.

<< Graphics`PlotField`

Part 1.

No responses needed in this section.

Part 2. Numerical Methods in a Helper Application

Enter the following command to define the function f that describes the right-hand side of the differential equation (dy/dt)= f [t, y].

f[t_, y_] :=  2*Cos[t] - t*y;
DE1  =  y'[t] == f[t, y[t]];

Construct the direction field for the differential equation.

fieldPlot  =  
        PlotVectorField[{1, 2*Cos[t] - t*y}, {t, -1, 8}, {y, -1, 3}, 
          ScaleFunction -> (1 &), PlotPoints -> 20, 
          Axes -> True, HeadLength -> 0, 
          HeadWidth -> 0];

Euler's Method

First, we use Euler's Method to solve the initial-value problem (IVP)

	(dy/dt) = 2 Cos[t] - t*y, with y[0] = 2

on the interval [0,8].  We set the number n of steps and the step size Δ, and we calculate the t-coordinates t[k] = k*Δ.

n  =  20;
Delta  =  8/n;
Clear[t];
t[k_] :=  k*Delta;

Next we define the recursive procedure for getting each y from the preceding point.  Then, we enter the initial condition, calculate the approximate solution points, and collect the points into a list.

Clear[y];
y[0]  =  2;
y[k_] :=  y[k]  =  
          y[k - 1]  +  N[f[t[k - 1], y[k - 1]]]*Delta;
Table[y[k], {k, n}];
Data[n]  =  Table[{t[k], y[k]}, {k, 0, n}];

Now we plot the approximation solution, giving it a name that depends on the number of points.  If the points look too small, just enlarge the plot.

Egraf[n]  =  
        ListPlot[Data[n], 
          PlotStyle -> {RGBColor[1, 0, 0], 
              PointSize[0.015]}];

After you have displayed the approximations for n = 20, 40, and 80, you may display them all together.  Scroll back to where n was set to 20, and change it, first to 40, then to 80.  Each time you plot, change the symbol and the color so you can tell the plots apart when you display them together.

Show[Egraf[20], Egraf[40], Egraf[80]];

Use the text cell below to describe the changes you see as n goes from 20 to 40 to 80.

The Improved Euler's Method

The following commands repeat what we did with Euler's Method, but this time the recursive procedure uses the Improved Euler (or Heun) Method.

n  = 20;
Delta  = 8/n;
Clear[t, y];
t[k_] :=  k*Delta;
y[0] = 2;
\(\(y[k_] := \(y[k]\  = \ 
          y[k - 1]\  + \ 
            N[0.5*\((f[t[k - 1], y[k - 1]] + 
                      f[t[k], y[k - 1] + f[t[k - 1], y[k - 1]]*Delta])\)]*
              Delta\);\)
\)

IEdata[n]  =  Table[{t[k], y[k]}, {k, n}];
IEgraf[n]  =  
        ListPlot[IEdata[n], 
          PlotStyle -> {RGBColor[1, 0, 0], 
              PointSize[0.015]}];

As you did before, scroll back to where n was set to 20, and change it, first to 40, then to 80.  Also, change the symbol and the color in the plot so you can tell the plots apart when you display them together.

After you have displayed the approximations for n = 20, 40, and 80, you may display them all together.

Show[IEgraf[20], IEgraf[40], IEgraf[80]];

Use the text cell below to describe the changes you see as n goes from 20 to 40 to 80.

A Mathematica Numerical ODE Solver

Mathematica's default numerical ODE solver uses methods by Adams and Gear.  Consult an advanced numerical analysis text for details.  Here we specify for Mathematica to use the popular Runge-Kutta-Fehlberg 4-5 order method.

Clear[y]
NDSolve[{DE1, y[0] == 2}, y, {t, 0, 8}, 
        Method -> RungeKutta];
Plot[Evaluate[y[t] /. %], {t, 0, 8}, 
        DisplayFunction -> Identity];
RKgraf  =  
        Show[%, fieldPlot, 
          DisplayFunction -> $DisplayFunction];

Now compare RKgraf with the results from the other numerical schemes.Vary the second entry to display each of the preceding plots along with RKgraf.

Show[RKgraf, Egraf[20]];

Answers to questions 5 and 6:

Part 3. The Logistic Equation

The following lines clear Mathematica's memory.  Please do not make alterations to these lines.

Unprotect[In, Out];
Clear[In, Out];
Protect[In, Out];

We define the logistic differential equation as DE2.

f[t_, y_] :=  y*(1 - y);
DE2  =  D[y[t], t] == f[t, y[t]];

Plot the direction field and Runge-Kutta solution starting at [0,0.1].

NDSolve[{DE2, y[0] == 0.1}, y, {t, 0, 8}];
solutionPlot  =  
        Plot[Evaluate[y[t] /. %], {t, 0, 8}, 
          PlotStyle -> RGBColor[1, 0, 0], 
          DisplayFunction -> Identity];

fieldPlot  =  
        PlotVectorField[{1, f[t, y]}, {t, 0, 8}, {y, 0, 1.5}, 
          ScaleFunction -> (1 &), PlotPoints -> 20, 
          Axes -> True, PlotPoints -> 20, AspectRatio -> 1, 
          DisplayFunction -> Identity];

RKplot  =  
        Show[solutionPlot, fieldPlot, 
          DisplayFunction -> $DisplayFunction];

DSolve[{DE2, y[0] == 0.1}, y[t], t]

Enter your exact solution formula as g[t].  You may use Mathematica for calculations if necessary.  Incidentally, e is predefined in Mathematica as the reserved variable E.

g[t_] :=  ?;

The following commands construct a plot of your function g[t] (in blue) and overlay the function plot on the numerical solution (red) and direction field plot.

gplot  =  
        Plot[g[t], {t, 0, 8}, PlotStyle -> RGBColor[0, 0, 1], 
          DisplayFunction -> Identity];
Show[RKplot, gplot, 
        DisplayFunction -> $DisplayFunction];

Scroll back and change the initial point to [0, 1.5], and re-plot the numerical solution.  (Note that the direction field need not be re-plotted.)  Recalculate your solution function g, and answer question 3 in the text cell below.

Answer:

Part 4. Another Initial Value Problem

Copy, paste, and edit Mathematica commands from Part 3 into the input cell below to carry out the calculations for Part 4.