ticks[min_,max_]:=Join[
Table[{i,If[i==0, ,Style[i,12]],{0.01,0}},{i,Ceiling[min],Floor[max], 2}],(* Numbered ticks *)
Table[{j, ,{0.01,0}},{j,Round[min],Round[max-1],1}]]         (* Un-numbered ticks *)


(* 2D Options *)
SetOptions[
Plot,
ImageSize->350,
AxesOrigin->{0,0},
AspectRatio->1,
AxesStyle->Directive[Medium, Bold],
Ticks->ticks
];

SetOptions[
ParametricPlot,
ImageSize->350,
AxesOrigin->{0,0},
AspectRatio->1,
AxesStyle->Directive[Medium, Bold],
Ticks->ticks
];

SetOptions[
VectorPlot,
ImageSize->350,
AxesOrigin->{0,0},
AspectRatio->1,
Frame->False,
Axes->True,
AxesStyle->Directive[Medium, Bold],
Ticks->ticks
];

(* 3D Options *)

SetOptions[
Plot3D,
ImageSize->600,
BoxRatios->1,Boxed->False,AxesOrigin->{0,0,0},
AxesStyle->Directive[Medium, Bold],
Ticks->ticks, 
Mesh->None,
ViewPoint->{2,0.9,1},
ViewVertical->{0,0,1}
];

SetOptions[
ContourPlot3D,
ImageSize->600,
BoxRatios->1,Boxed->False,AxesOrigin->{0,0,0},
AxesStyle->Directive[Medium, Bold],
Ticks->ticks, 
Mesh->None,
ViewPoint->{2,0.9,1},
ViewVertical->{0,0,1}
];

SetOptions[
ParametricPlot3D,
ImageSize->600,
BoxRatios->1,Boxed->False,AxesOrigin->{0,0,0},
AxesStyle->Directive[Medium, Bold],
Ticks->ticks, 
Mesh->None,
ViewPoint->{2,0.9,1},
ViewVertical->{0,0,1}
];

SetOptions[
RegionPlot3D,
ImageSize->600,
BoxRatios->1,Boxed->False,AxesOrigin->{0,0,0},
AxesStyle->Directive[Medium, Bold],
Ticks->ticks, 
Mesh->None,
ViewPoint->{2,0.9,1},
ViewVertical->{0,0,1}
];

SetOptions[
VectorPlot3D,
ImageSize->350,
AxesOrigin->{0,0,0},
AspectRatio->1,
BoxRatios->1,
Axes->True,
AxesStyle->Directive[Medium, Bold],
Ticks->ticks,
ViewPoint->{2,0.9,1},
ViewVertical->{0,0,1}
];

SetOptions[
RevolutionPlot3D,
ImageSize->350,
AxesOrigin->{0,0,0},
AspectRatio->1,
BoxRatios->1,
Boxed->False,
Axes->True,
AxesStyle->Directive[Medium, Bold],
Ticks->ticks,
ViewPoint->{2,0.9,1},
ViewVertical->{0,0,1}
];

Euler’s Method

The Idea

Suppose you want to solve the differential equation (dy/dx)=(e^(-(x^2))) where y(0)=1.  Well, first off, good luck!  You’ll quickly realize that it has a unique solution, namely

({{({{y(x)=1+(∫_0^x)(e^(-(u^2)))du}}),}})

but unfortunately this is not one of our elementary integrals.  Mathematica simply defines a new function,

Integrate[Exp[-u^2], {u, 0, x}]

So, we need to approximate.  This is where Euler’s Method can help.  The idea:  

Consider the initial value problem

({{({{(dy/dx)=f(x,y),       y((Subscript[x, 1]))=(Subscript[y, 1]).}})}})

Pick an initial point, say ((Subscript[x, 1]), (Subscript[y, 1]))

The slope through this point ((Subscript[x, 1]), (Subscript[y, 1])) is given by f((Subscript[x, 1]), (Subscript[y, 1])).

Move a small distance along the line segment through this point to the next point, call it creatively ((Subscript[x, 2]), (Subscript[y, 2])).  Note that (Subscript[y, 2]) is simply an approximation to y((Subscript[x, 2])), that is (Subscript[y, 2])≈y((Subscript[x, 2])).

At the new point ((Subscript[x, 2]), (Subscript[y, 2])), we have a new slope given by f((Subscript[x, 2]), (Subscript[y, 2])).

Move a small distance along the line segment through this point to the next point, call it creatively ((Subscript[x, 3]), (Subscript[y, 3])).  Note that (Subscript[y, 3]) is simply an approximation to y((Subscript[x, 3])), that is (Subscript[y, 3])≈y((Subscript[x, 3])).

Etc.

A more Mathematica version of this process is:

Pick a step size, say h>0.  Usually this is “small.”

Pick an initial point, say ((Subscript[x, 1]), (Subscript[y, 1])).

Your next point is ((Subscript[x, 2]), (Subscript[y, 2])), where (Subscript[x, 2])=(Subscript[x, 1])+h  (just step over a bit) and (Subscript[y, 2])=(Subscript[y, 1])+ h · f((Subscript[x, 1]), (Subscript[y, 1])).

Etc

Your n+(1^rst) point is ((Subscript[x, n+1]), (Subscript[y, n+1])), where (Subscript[x, n+1])=(Subscript[x, n])+h  (just step over a bit) and (Subscript[y, n+1])=(Subscript[y, n])+ h · f((Subscript[x, n]), (Subscript[y, n])).

Hence, the numbers (Subscript[y, 1]), (Subscript[y, 2]), (Subscript[y, 3]), ..., (Subscript[y, n]), ... form an approximate solution to y(x).  If we pair them with appropriate x values, we can even graph them as points.

Example

Let us return to the example.  Suppose you want to approximate the solution to the differential equation (dy/dx)=(e^(-(x^2))) where y(0)=1.

Pick a step size, say h>0.  For us, let us choose h=1.

Pick an initial point, say ((Subscript[x, 1]), (Subscript[y, 1]))=(0,1) since y(0)=1.

Your next point is ((Subscript[x, 2]), (Subscript[y, 2])), where (Subscript[x, 2])=(Subscript[x, 1])+h=0+1 = 1.  (just step over a bit) and

({{({{ (Subscript[y, 2])=(Subscript[y, 1])+ h · f((Subscript[x, 1]), (Subscript[y, 1]))=1+1·f(0,1)=1+1·1=2.}})}})

Your next point is ((Subscript[x, 3]), (Subscript[y, 3])), where (Subscript[x, 3])=(Subscript[x, 2])+h=1+1 = 2.  (just step over a bit) and

({{({{ (Subscript[y, 3])=(Subscript[y, 2])+ h · f((Subscript[x, 2]), (Subscript[y, 2]))=2+1·f(1,2)=2+1·(e^(-(1^2)))≈2.3678794411714423`.}})}})

Etc

Your n+(1^rst) point is ((Subscript[x, n+1]), (Subscript[y, n+1])), where (Subscript[x, n+1])=(Subscript[x, n])+h  (just step over a bit) and (Subscript[y, n+1])=(Subscript[y, n])+ h · f((Subscript[x, n]), (Subscript[y, n])).

Note that we can use Mathematica to compute some values for us!

N[Pi]
N[Pi, 100]

2+Exp[-1]  (* This will give us an exact value!*)
N[2+Exp[-1]]  (* Remeber that the function N[] gives you an approximation! *)

A better way would be to define a function for f(x,y):

f[x_, y_]=N[Exp[-x^2]]  (* I've added the N[] so that I don't have to keep typing it *)

f[0,1]
f[1,2]

Here is the Code to Perform Euler’s Method

f[x_, y_]=Exp[-x^2];  (* Define the right hand side of our derivative, (dy/dx)=f(x,y) *)

Note that the following Mathematica Module is a new function that can be called.  It accepts the following (local) variables:

funct:  Euler’s Method will be used on this function.

initX: The left endpoint of the interval.

initY: This is (Subscript[y, 0]), the initial value of the dependent variable.

finalX: The right endpoint of the interval.

numsteps: The number of steps or subintervals.  Note that the step size is calculated in the module.

basicEuler[funct_,initX_, initY_,finalX_,numsteps_]:=Module[{f=funct,beginPoint =initX,w0=initY, endPoint = finalX,n=numsteps},
x=beginPoint;
w=w0;
h=(endPoint-beginPoint)/n;  (*The is the step size. *)
xList={x};  (*Initialize our list of t-values *)
wList ={w}; (*Initialize our list of w-vaules: the approximations to the solution *)

(* Here we loop from 1 to n performing Euler's method and appending the results to the lists *)
For[j=1,j<n+1,j++,
m1=h*f[x,w];
x=x+h;
w=w+m1;
xList=Append[xList,x];
wList=Append[wList,w];
];
Return[{xList,wList}]
]

Here is the Code to Perform Euler’s Method:  Case h = 1 (number of steps is 4)

my1=basicEuler[f, 0,1, 4,4]

This generated a list with two elements: the list of t-values and the list of corresponding approximation values.  In order to graph the approximate solution we will need to create a list of ordered pairs.  Do this with the Table command, remembering that we had 4 steps (so 5 points).

xyList1=Table[{my1[[1,k]],my1[[2,k]]},{k,1,5}]

Now we can plot the approximation to the solution of the IVP:

eulerHisOne  =(* Note that you can save plots for later! *)
ListPlot[xyList1,
Joined->True,
Mesh->All,PlotMarkers->{Automatic,Medium},
PlotStyle->{Blue,Thick},
AxesOrigin->{0,0},
PlotLegends->LineLegend[{Blue}, {"H = 1"}]
];

(* Plot our list of ordered pairs *)
Show[eulerHisOne,
PlotRange->{-0.5, 2.6},
AxesLabel->{x, y},
AxesStyle->Bold,
PlotLabel->Style[Text["Euler's Method"], 16 ]
]

N[xyList1]

Here is the Code to Perform Euler’s Method:  Case h = 0.5 (number of steps is 8)

my2=basicEuler[f, 0, 1,4,8];
xyList2 = Table[{my2[[1,k]],my2[[2,k]]},{k,1,9}];

eulerHisOneHalf=(* Note that you can save plots for later! *)
ListPlot[xyList2,
Joined->True,
Mesh->All,PlotMarkers->{"▲",14}, (* Note that "▲" is found in the SpecialCharacters menu *)
PlotStyle->{Red,Thick},
AxesOrigin->{0,0},
PlotLegends->LineLegend[{Red}, {"H = 0.5"}]
];

(* Plot our list of ordered pairs *)
Show[eulerHisOneHalf,
PlotRange->{-0.5, 2.6},
AxesLabel->{x, y},
AxesStyle->Bold,
PlotLabel->Style[Text["Euler's Method"], 16]
]

Several Euler’s Method Plots Compared

Show[
eulerHisOne, eulerHisOneHalf,
PlotRange->{-0.5,2.6},
AxesLabel->{x, y},
PlotLabel->Style[Text["Euler's Method Comparison"], 16]
]

Direction Field for this Problem

dField=VectorPlot[
{1, f[x,y]}, {x, -2, 4}, {y, -2, 2.6},
Axes->True,
AxesLabel->{x,y},
VectorScale->{Small, Small, None},
VectorStyle->Gray,
StreamPoints->Coarse,
StreamStyle->Directive[Darker[Green],Thick,"Line"],
StreamScale->None
];  (* Note that if I don't use a semi-colon here, it will save the plot and still show the plot. *)

Show[dField,
PlotLabel->Style["Direction Field, (dy/dx) = (e^(-(x^2)))", Bold, 16]
]

Comparison of the Direction Field with Euler’s Method.

Show[eulerHisOne,eulerHisOneHalf,dField,
PlotLabel->Style["Direction Field  and Euler's Method", Bold, 16],
PlotRange->{{-2.5, 4.5}, {-2.5, 2.5}},
AspectRatio->1]

Assignment

Use Mathematica to complete the following problems.  Your notebook should include written comments using appropriate styles, graphs that are clearly labeled and work that is well organized.

Problem

Consider the differential equation

({{({{(dy/dx)=(e^x) cos(y).}})}})

Find all of the equilibria (places where the slope is zero and the solution is a constant function)...note that there are a LOT of them.  Please give a formula for all equilibria.

Generate a direction field which shows where at least 3 perhaps 4 of these equilibria are located. Use -4≤x≤6 and something appropriate for y.

Use Euler’s Method to solve estimate the solutions with the following step sizes: h=1, 0.5, 0.25, 0.125 over the interval 0≤x≤4. Choose your initial condition to be near, but just above the largest negative equilibrium (y(0)=-(3/2) would be a good one to use).

Create plots using Grid to illustrate the 4 cases generated by the various h.  Make sure that your horizontal axis ends at x=6.

Explain your results as h changes.  Use your own words, don’t plagiarize others!

Overlay a plot using Euler’s Method when h=0.125 with the direction field.

Finally, use DSolve or NDSolve to find a solution for this differential equation using your same initial guess as above.  Plot this and compare to Euler’s Method when h=0.125.  Use your own words, don’t plagiarize others!!

Existence and Uniqueness of Solutions

Here is some of the theory that says when it is “ok” to look for an approximation to a solution using a numerical method.  The “problems” listed below may be discussed in class, but do not need to be completed as part of this lab.

Existence Theorem

Suppose f(t,y) is a continuous function in a rectangle of the form {(t,y)|a<t<b, c<y<d} in the t y -plane.  If ((Subscript[t, 0]),(Subscript[y, 0]))  is a point in this rectangle, then there exists an ϵ>0 and a function y(t) defined for (Subscript[t, 0])-ϵ<t<(Subscript[t, 0])+ϵ that solves the initial value problem

({{({{(dy/dt)=f(t,y),       y((Subscript[t, 0]))=(Subscript[y, 0]).}})}})

Problem

Consider the differential equation

({{({{(dy/dt)=3 (y^(2/3)),    y(0)=0.}})}})

Find all solutions.  What is the domain of existence?  How large can ϵ be in the existence theorem?

Problem

Consider the differential equation ({{(dy/dx)=(y/(1-x)),    y(0)=1.}})

Find all solutions.   What is the domain of existence?  How large can ϵ be in the existence theorem?

Generate a direction field.  Ensure that the direction field uses -2≤x≤2.

Now consider the initial value problem

({{({{(dy/dx)=(y/(1-x)),    y(2)=1.}})}})

What is the domain of existence?  How large can ϵ be in the existence theorem?

Uniqueness Theorem

Suppose f(t,y) and ((∂f)/(∂y)) are continuous functions in a rectangle of the form {(t,y)|a<t<b, c<y<d} in the t y -plane.  If ((Subscript[t, 0]),(Subscript[y, 0]))  is a point in this rectangle and if (Subscript[y, 1])(t) and (Subscript[y, 2])(t) are two functions that solve the initial-value problem,

({{({{(dy/dt)=f(t,y),       y((Subscript[t, 0]))=(Subscript[y, 0]).}})}})

for all t in the interval (Subscript[t, 0])-ϵ<t<(Subscript[t, 0])+ϵ (where ϵ is a positive number), then

({{({{(Subscript[y, 1])(t)=(Subscript[y, 2])(t)}})}})

for  (Subscript[t, 0])-ϵ<t<(Subscript[t, 0])+ϵ.  That is, the solution to the initial-value problem is unique.

Problem

Consider the differential equation

({{({{(dy/dt)=(y+1)  (y-1).}})}})

Find any solutions (there are 3 solutions...)  Note that one solution needs partial fractions, you can use the Apart[  ] command.

Apart[1/((y+1) (y-1))]

Why does this not violate the uniqueness theorem?

Consider the initial condition y(0)=(1/2).   Find this solution.  Is it unique?

Consider the direction field.  For 0<y(0)=(Subscript[y, 0])<1, what happens as t->∞?  Can you prove this?

points={{0,1}, {0,1/2}, {0,0},{0, -1},{0,3/2},{0,-3/2}};
VectorPlot[
{1, (y+1) (y-1)}, {t, -2, 2}, {y, -2, 2},
AxesLabel->{t,y},
Axes->True,
VectorScale->{Small, Small, None},
VectorStyle->Gray,
StreamPoints->points,
StreamStyle->Directive[Darker[Green],Thick,"Line"],
StreamScale->None,
PlotLabel->"Direction Field Using Particular Solution Curves",
Epilog->{PointSize[Large], Point[points]}
]

Suppose you have two solutions based upon 0<(Subscript[y, 1])(0)<1 and 0<(Subscript[y, 2])(0)<1.  If (Subscript[y, 1])(0) < (Subscript[y, 2])(0), sketch the two curves below.

points={{0,1},  {0,-1}};
VectorPlot[
{1, (y+1) (y-1)}, {t, -2, 2}, {y, -1.5, 1.5},
AxesLabel->{t,y},
Axes->True,
VectorScale->{Small, Small, None},
VectorStyle->Gray,
StreamPoints->points,
StreamStyle->Directive[Darker[Green],Thick,"Line"],
StreamScale->None,
PlotLabel->"Direction Field Using Particular Solution Curves",
Epilog->{PointSize[Large], Point[points]}
]