Second-order Taylor method

The following implements the second-order Taylor method.  Note that we use Mathematica's symbolic capabilities to compute the required partial derivatives of the function f(t, y).

f[t_, y_] :=  10  y (1 - y);
ft[t_, y_] = D[f[t, y], t];
fy[t_, y_] = D[f[t, y], y];
t0 = 0. ;
y0 = 0.5;
tf = 2;
n = 20;
h = (tf - t0)/n;
tylist = {{t0, y0}};
Do[
tylist = 
          Append[tylist, 
{t + 
                h, 
tylist[[-1, 2]] + 
                h*f[t, tylist[[-1, 2]]] + 
                0.5*h^2*(ft[t, tylist[[-1, 2]]] + 
                      f[t, tylist[[-1, 2]]]*
                        fy[t, tylist[[-1, 2]]])}], {t, t0, tf - h, 
          h}];
approxsolnplot = 
ListPlot[
        tylist, 
PlotJoined -> 
          True, 
PlotRange -> 
          All, 
PlotStyle -> 
          RGBColor[1, 0, 0]]

Note that I have used the option PlotJoined->True in the ListPlot command above.

We next use code from a previous day to plot a slope field so we can lay the approximate solution on top.

tmin = 0;
tmax = 2;
ymin = 0;
ymax = 2;
n = 20;
tsl = 0.03;
tanseg[t_, 
          y_] = 
Line[{
{t, y} - 
              tsl*{1, f[t, y]}/Sqrt[1 + f[t, y]^2], {t, y} + 
              tsl*{1, f[t, y]}/
                  Sqrt[1 + 
                      f[t, y]^2]
}
];
slopefield = 
Show[
        Graphics[
Table[
            tanseg[t, y], {t, tmin, 
              tmax, (tmax - tmin)/n}, 
{y, ymin, 
              ymax, (ymax - ymin)/n}]], 
Frame -> 
          True, 
AspectRatio -> Automatic]

Show[slopefield, approxsolnplot, 
      PlotRange -> {{0, 2}, {0, 2}}]

Comparisons with Runge-Kutta

With the next set of inputs, we build a list in which each element has a time, the second-order Taylor method approximation, and two Runge-Kutta approximations of order two (specifically, midpoint and modified Euler.)

f[t_, y_] :=  t - y;
ft[t_, y_] = D[f[t, y], t];
fy[t_, y_] = D[f[t, y], y];
t0 = 0. ;
y0 = 1;
tf = 2;
n = 20;
h = (tf - t0)/n;
tylist = {{t0, y0, y0, y0}};
Do[
tylist = 
Append[
tylist, 
{t + 
                h, 
tylist[[-1, 2]] + 
                h*f[t, tylist[[-1, 2]]] + 
                0.5*h^2*(ft[t, tylist[[-1, 2]]] + 
                      f[t, tylist[[-1, 2]]]*
                        fy[t, tylist[[-1, 
                              2]]]), 
tylist[[-1, 
                    3]] + 
                h*f[t + 0.5*h, 
                    tylist[[-1, 3]] + 
                      0.5*h*f[t, 
                          tylist[[-1, 
                              3]]]], 
tylist[[-1, 
                    4]] + 0.5*h*f[t, tylist[[-1, 4]]] + 
                0.5*h*f[t + h, 
                    tylist[[-1, 4]] + 
                      h*f[t, tylist[[-1, 
                              4]]]]}
], 
{t, t0, tf - h, h}];
\(PaddedForm[
      TableForm[tylist, 
        TableHeadings -> {None, {"t", "Taylor", "Midpoint", "Modified Euler"}}], {17, 16}]\)

Next, we add Nick's Crazzee Method to the mix.

f[t_, y_] :=  t - y;
ft[t_, y_] = D[f[t, y], t];
fy[t_, y_] = D[f[t, y], y];
t0 = 0. ;
y0 = 1;
tf = 2;
n = 4;
h = (tf - t0)/n;
tylist = {{t0, y0, y0, y0, y0}};
Do[
tylist = 
Append[
tylist, 
{t + 
                h, 
tylist[[-1, 2]] + 
                h*f[t, tylist[[-1, 2]]] + 
                0.5*h^2*(ft[t, tylist[[-1, 2]]] + 
                      f[t, tylist[[-1, 2]]]*
                        fy[t, tylist[[-1, 
                              2]]]), 
tylist[[-1, 
                    3]] + 
                h*f[t + 0.5*h, 
                    tylist[[-1, 3]] + 
                      0.5*h*f[t, 
                          tylist[[-1, 
                              3]]]], 
tylist[[-1, 
                    4]] + 0.5*h*f[t, tylist[[-1, 4]]] + 
                0.5*h*f[t + h, 
                    tylist[[-1, 4]] + 
                      h*f[t, tylist[[-1, 
                              4]]]], 
tylist[[-1, 
                    5]] - 1*h*f[t, tylist[[-1, 5]]] + 
                2*h*f[t + 0.25*h, 
                    tylist[[-1, 5]] + 
                      0.25*h*f[t, 
                          tylist[[-1, 
                              5]]]]}
], 
{t, t0, tf - h, h}];
\(PaddedForm[
      TableForm[tylist, 
        TableHeadings -> {None, {"t", "Taylor", "Midpoint", "Modified Euler", "Nick's Crazzee"}}], {17, 16}]\)

The function we use above doesn't provide much of a challenge since the partial derivatives of f(t, y) are constants.  Below are examples with more interesting functions.  Note that I've substituted Euler's method for the second-order Taylor method in the following.

f[t_, y_] :=  1 - y^2;
t0 = 0. ;
y0 = 0;
tf = 1;
n = 10;
h = (tf - t0)/n;
tylist = {{t0, y0, y0, y0, y0}};
Do[
tylist = 
Append[
tylist, 
{t + 
                h, 
tylist[[-1, 2]] + 
                h*f[t, tylist[[-1, 
                        2]]], 
tylist[[-1, 3]] + 
                h*f[t + 0.5*h, 
                    tylist[[-1, 3]] + 
                      0.5*h*f[t, 
                          tylist[[-1, 
                              3]]]], 
tylist[[-1, 
                    4]] + 0.5*h*f[t, tylist[[-1, 4]]] + 
                0.5*h*f[t + h, 
                    tylist[[-1, 4]] + 
                      h*f[t, tylist[[-1, 
                              4]]]], 
tylist[[-1, 
                    5]] - 1*h*f[t, tylist[[-1, 5]]] + 
                2*h*f[t + 0.25*h, 
                    tylist[[-1, 5]] + 
                      0.25*h*f[t, 
                          tylist[[-1, 
                              5]]]]}
], 
{t, t0, tf - h, h}];
PaddedForm[
      TableForm[tylist, 
        TableHeadings -> {None, {"t", "Euler", "Midpoint", "\
Modified Euler",  "Nick's Wacko"}}], {17, 16}]

f[t_, y_] :=  
        Exp[t] (y^3 + t y^3 + 1)/3/
              y^2/(t Exp[t] - 6);
t0 = 0. ;
y0 = 1;
tf = 1;
n = 10;
h = (tf - t0)/n;
tylist = {{t0, y0, y0, y0, y0}};
Do[
tylist = 
Append[
tylist, 
{t + 
                h, 
tylist[[-1, 2]] + 
                h*f[t, tylist[[-1, 
                        2]]], 
tylist[[-1, 3]] + 
                h*f[t + 0.5*h, 
                    tylist[[-1, 3]] + 
                      0.5*h*f[t, 
                          tylist[[-1, 
                              3]]]], 
tylist[[-1, 
                    4]] + 0.5*h*f[t, tylist[[-1, 4]]] + 
                0.5*h*f[t + h, 
                    tylist[[-1, 4]] + 
                      h*f[t, tylist[[-1, 
                              4]]]], 
tylist[[-1, 
                    5]] - 1*h*f[t, tylist[[-1, 5]]] + 
                2*h*f[t + 0.25*h, 
                    tylist[[-1, 5]] + 
                      0.25*h*f[t, 
                          tylist[[-1, 
                              5]]]]}
], 
{t, t0, tf - h, h}];
PaddedForm[
      TableForm[tylist, 
        TableHeadings -> {None, {"t", "Euler", "Midpoint", "\
Modified Euler",  "Nick's Texas Two-Step"}}], {17, 16}]