Relationship Between Direction Fields and Euler's Method

Example 1: y'=sin((x^2))

Differential Equation

Clear[x,y,f]
f[x_,y_]=Sin[x^2];

To solve this by hand would require that we separate and integrate, but the integral of sin((x^2)) is not one that we know how to do.

∫dy= ∫sin((x^2))dx

DSolve[y'[x]==f[x,y],y[x],x]

Mathematica can solve it, but we might not know what a Fresnel S function is. We can use a direction field to get a global undestanding of the solution, even in cases where Mathematica can't solve the differential equation.

Direction Field

I've put grid lines on for what is to come next. I've also fiddled with the layout (number of vectors and length of vectors) to make what follows about Euler's method easier to understand. The standard direction field syntax is in the second example.

fieldeqn=VectorPlot[{1,f[x,y]},{x,-2,2},{y,-2,2},Axes->True,FrameLabel->{"x","y"},VectorPoints->9,VectorScale->{Small,Small,None}]

Explanation of Euler's Method Solution

Here is Euler's method solution to the initial value problem (remember, Euler's method solves initial value problems) for the above differential equation combined with the initial condition (Subscript[x, 0])=-2, (Subscript[y, 0])=0.

Clear[x,y,h]
x[0]=-2.0;
y[0]=0.0;
h=0.5;
x[n_]:=x[n]=x[n-1]+h
y[n_]:=y[n]=h*f[x[n-1],y[n-1]]+y[n-1]

data=Table[{x[n],y[n]},{n,0,8}]
TableForm[data]

eulersol1=ListPlot[data,PlotStyle->{PointSize[0.02],Red},AxesLabel->{"x","y"}]
eulersol2=ListPlot[data,Joined->True,PlotStyle->{Purple,Thick},AxesLabel->{"x","y"}]

euler1=Show[fieldeqn,eulersol2,eulersol1]

Euler's method acts as if the direction field gives the solution a "kick" in the direction of the arrow at each step (which is where the gridlines are that I added earlier). You move in a straight line in the direction of the arrow until you have moved a distance h in the x direction, then you get another "kick" from the direction field. Between the kicks the Euler method solution gets no input from the direction field. You can make the graph bigger to really see this.

Decreasing the step size makes the solution feel more "kicks" from the direction field, and produces a more accurate solution.

Clear[x,y,h]
x[0]=-2.0;
y[0]=0.0;
h=0.1;
x[n_]:=x[n]=x[n-1]+h
y[n_]:=y[n]=h*f[x[n-1],y[n-1]]+y[n-1]

data=Table[{x[n],y[n]},{n,0,8*5}]
TableForm[data]

eulersol1=ListPlot[data,PlotStyle->{PointSize[0.02],Red},AxesLabel->{"x","y"}]
eulersol2=ListPlot[data,Joined->True,PlotStyle->{Purple,Thick},AxesLabel->{"x","y"}]

Show[euler1,eulersol2,eulersol1]

Example 2: y'=(y^3)-(x^3)

Differential Equation

Clear[x,y,f]

DSolve[y'[x]==y[x]^3+x^3,y[x],x]

This differential equation is not separable, but we can use a direction field to get a global undestanding of the solution, even in cases where Mathematica can't solve the differential equation. And we can still use Euler's method to get a numerical solution to an initial value problem (IVP).

Direction Field

The VectorScale command makes the lines all one length, and makes them easier to see.

fieldeqn=VectorPlot[{1,y^3-x^3},{x,-2,2},{y,-2,2},Axes->True,FrameLabel->{"x","y"},VectorScale->{Tiny,Small,None}]

Euler's Method Solution

Here is Euler's method solution to the initial value problem (remember, Euler's method solves initial value problems) for the above differential equation combined with the initial condition (Subscript[x, 0])=-1, (Subscript[y, 0])=0.

Clear[x,y,h]
f[x_,y_]=y^3-x^3;
x[0]=-1.0;
y[0]=0.0;
h=0.5;
x[n_]:=x[n]=x[n-1]+h
y[n_]:=y[n]=h*f[x[n-1],y[n-1]]+y[n-1]

data=Table[{x[n],y[n]},{n,0,100}];
TableForm[data]

eulersol1=ListPlot[data,PlotStyle->{PointSize[0.02],Red},AxesLabel->{"x","y"}];

euler1=Show[fieldeqn,eulersol1]

Notice the overflow that occurs in this case, which can be for two reasons:
1) The solution is becomming infinite,
2) Newton's method is suffering from round-off error.

Regardless of the overflow errors, the Newton's method solution seems to be following the direction field pretty well in the region we have plotted. But we should still reduce h and see if we get different results.

Clear[x,y,h]
f[x_,y_]=y^3-x^3;
x[0]=-1.0;
y[0]=0.0;
h=0.02;
x[n_]:=x[n]=x[n-1]+h
y[n_]:=y[n]=h*f[x[n-1],y[n-1]]+y[n-1]

data=Table[{x[n],y[n]},{n,0,100}];
TableForm[data]

eulersol1=ListPlot[data,PlotStyle->{PointSize[0.02],Red},AxesLabel->{"x","y"}];

euler1=Show[fieldeqn,eulersol1]

Very different than what we had before. I would still be wary of trusting that value for y(1) without more investigation.