Initialization Cells (Code Needed for Most Plot3D Statements)

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}
];

First Order Differential Equations

Modeling

What Is Modeling

Modeling in this class will be the use of derivatives and equations to approximate a “real” object.  To accomplish this you should
	1.	Clearly state all assumptions on which the model will be based.  These assumptions should describe the relationships among the variables to be studied.
	2.	Completely describe (give names to) the independent variables (usualy time) or dependent variables (number of animals at any given time or the amount of a reactant at any given time) and parameters (these should be know - like a death or a birth rate or a chance of infection) to be used in the model.
	3.	Derive equation(s) which relate the variables together in terms of variables and parameters and most importantly derivatives (i.e. rate of change of... or velocity, or acceleration, or rate of increase of..., etc).

The next interesting goal is to run some sort of experiment which gathers data and allows us to somehow test the usefulness of the model.

We might also try to “solve” the differential equation.  Good luck with that!  However, there are three techniques we’ll follow
	A.	Analytical methods allow us (sometimes) to find formulas for the dependent variable(s) in terms of the independent variable(s) and the parameters.  Note that we do not need to have estimates for the parameters.  This is purely symbolic manipulation.
	B.	Qualitative methods allow us to gain information by creating a rough sketch of the graph in terms of the independent variable(s).  It gives us a rough idea as to the end behavior of the model.  We probably don’t need estimates of the parameters.
	C.	Numerical methods involves arithmetic (either by hand or by computer) to approximate the behavior of the model.  We need to have estimates of the parameters.

Definitions

An ordinary differential equation is an equation that relates an unknown function, one independent variable, and one or more of the functions derivatives.

A partial differential equation is a differential equation which has two or more independent variables.  (Different class!)

Notation:

	F(x, y(x), y'(x), y''(x), ..., (y^n)(x)) = 0
	This represents an (n^th) order differential equation with independent variable x.
	
	(y^((n)))(x) = G(x, y(x), y'(x), ..., (y^((n-1)))(x))		
	This is called the normal form for a (n^th) order differential equation with independent variable x.
	
An ordinary differential equation is called autonomous if it does not explicitly depend upon time.  Otherwise it is called nonautonomous.
	
Initial Value Problems are differential equations with extra information or initial conditions:
	y(a)=(Subscript[b, 0]),   y'(a)=(Subscript[b, 1]),   y''(a)=(Subscript[b, 2]), ..., (y^((n-1)))(a)=(Subscript[b, n-1]).

A solution to F(x, y(x), y'(x), y''(x), ..., (y^n)(x)) = 0 is a function u(x) that is continuous on an interval I such that u', u'', ..., (u^((n))) exist on I and F(x, u(x), u'(x), u''(x), ..., (u^n)(x)) = 0 is true for all x∈I.

A particular solution is found by using the initial values in an initial value problem.

A general solution is one that gives us a “family” of solutions when we fail to have initial conditions.

An equilibrium solution occurs when the rate(s) of change are zero.

Examples

Classify the following differential equations:
	
	1.	y'(x) = cos(x)
	2.	y''(x)=cos(x)
	3.	y''(x) + 3 y'(x) = cos(x)
	4.	y'(x)=cos(y)

Example

Verify that (Subscript[y, 1])(x) = (1/(x^2)) and (Subscript[y, 2])(x)=((ln(x))/(x^2)) are solutions to (x^2) y'' + 5 x y' + 4 y = 0.

Example

Verify that  (Subscript[y, 2])(x)=((ln(x))/(x^2)) is a solution to (x^2) y'' + 5 x y' + 4 y = 0 using Mathematica.

We need to define a function that we can then differentiate.  Notice the underscore to the left of the equal sign after the independent variable x and the complete lack of any underscore to the right of the equal sign.  Also, the natural logarithm is Log[x].  For a different base, use Log[x, b].

(Subscript[y, 2])[x_]=Log[x]/x^2

For differentiation, you can use (Subscript[y, 2])'[x] and (Subscript[y, 2])''[x].  This is equivalent to D[(Subscript[y, 2])[x], x] and D[(Subscript[y, 2])[x], {x, 2}], respectively.  In general, the (n^th) partial derivative of the function f is given by D[f, {x, n}].

(Subscript[y, 2])'[x]
(Subscript[y, 2])''[x]
(Subscript[dy, 2])[x_]=D[(Subscript[y, 2])[x],x]
D[D[(Subscript[y, 2])[x],x], x]
(Subscript[ddy, 2])[x_]=D[(Subscript[y, 2])[x],{x, 2}]
D[(Subscript[y, 2])[x],{x, 20}]

Finally, plug them into the equation and see if it is a solution.

x^2 (Subscript[ddy, 2])[x] + 5 x (Subscript[dy, 2])[x]+4(Subscript[y, 2])[x]
Simplify[%]

Doesn't look like it worked...so let's use the magical Simplify[] command!

Simplify[x^2 (Subscript[ddy, 2])[x] + 5 x (Subscript[dy, 2])[x]+4(Subscript[y, 2])[x]]

((4 Log[x])/(x^2))+(x^2) (-(5/(x^4))+((6 Log[x])/(x^4)))+5 x ((1/(x^3))-((2 Log[x])/(x^3)))//Simplify

FullSimplify[x^2 (Subscript[ddy, 2])[x] + 5 x (Subscript[dy, 2])[x]+4(Subscript[y, 2])[x]]

Example

Find all r such that y=(e^rx) is a solution to the differential equation 3 y'' + 3y' -4 y= 0.

Malthusian Growth

A very (naive) basic model for population growth is the following:

	The rate of growth of the population is proportional to the size of the population.

	A.	Write a model using this assumption and k as the parameter and P as the dependent variable.
	B.	Show that a family of solutions is P(t)=C (e^(k t)) where C is a constant.
	C.	Sketch a qualitative graph of the population over the  interval [0, ∞).
	D.	Why is this an unreasonable model?  When might it be reasonable?
	E.	Solve the initial value problem if P(0) = 10.
	F.	Graph several typical solutions for varying initial populations (i.e. C=10, 2, or 15).  For all such solutions, assume that k=2 and 0≤t≤1.
	G.	Find the equilibrium solution (when the rate of change is zero).

Plot[
{
10 Exp[2 t],
2 Exp[2 t],
15 Exp[2 t]
},
{t,0,1},
Ticks->Automatic, (* This overides my default behavior for tickmarks *)
PlotStyle->{{Blue, Thick}, {Red, Dotted, Thick}, {Darker[Green], Dashed,Thick}},
AxesLabel->{"t", "P(t)"},
PlotLabel->"Malthusian Growth",
LabelStyle->Directive[Bold, Purple]
]

Example

The line tangent to the graph of g at the point (x, y) intersects the x-axis at the point ((x/2), 0).  Find a differential equation that models this situation.

Limited Growth

In this type of model, we assume that there is some sort of limiting factor to the growth of the population, such as limited resources.  Thus when the population is below a certain threshold, the derivative (the growth rate) is positive, when the population is above a certain threshold, the deriative (the growth rate) is negative.

The assumptions
	1.	(dP/dt)≈k P when P is small
	2.	(dP/dt)<0 when P>N, the population threshold.

One method of modeling this behavior is (dP/dt)=k (1 - (P/N)) P.  This is called the logistic population model.
	A.	Find the equilibria.
	B.	Perform a qualitative analysis.  That is, sketch the equilibria and typical solution curves.

Homework - For Software, Use Mathematica.

Pages: 	14 - 21
Problems: 	1, 2, 3, 5, 7, 8, 11, 12, 17

Separation of Variables

Basic Definitions

The first order differential equation (dy/dx)=H(x,y) is called separable if H(x,y) = g(x) h(y).

Recall:  If y=f(x), then the differential of y is dy=f'(x) dx, so we can represent the derivative as 
		f'(x)=(dy/dx)=Quotient of dy and dx =(dy/dx)=Derivative of y with respect to x.

Theorem

Suppose (dy/dx)=g(x) h(y) (i.e. we have a separable differential equation).
Then

({{({{(dy/(h(y)))=g(x) dx   ⟹   }})∫(dy/(h(y))) = ∫g(x)dx + C.}})

Example

Consider the differential equation (dy/dx)=3 y.
	A.	Solve the initial value problem when y(0)=1.
	B.	Solve the initial value problem when y(0)=-1.

Example

Consider the differential equation (dy/dx)=3 (y^2).
	A.	Solve the initial value problem when y(0)=c.
	B.	Solve the initial value problem when y(0)=0.  (Oopsie!  What just happened?)
Remark:  When writing a general solution for this problem, you must include the constant part (i.e. the one in A) as well as the singular solution where no value of c gives zero (i.e. the one in B).

Partial Fractions Example

Remember Partial Fraction Decomposition?

	(1/(1-(x^2)))=(1/((1-x)(1+x)))=(A/(1-x))+(B/(1+x)) 
⟹ Clear Fractions by multiplying by (1-(x^2))
({{1=A(1+x)+B(1-x)}, {1=A+B + (A-B)x}})  
⟹Solve the system of equations
({{       1=A+B}, {       0=A-B}})
⟹A=(1/2),  B=(1/2).
As luck would have it, Mathematica has a built-in function for this...

Apart[(1/(1-x^2))]

Example

Solve the differential equation (1-(x^2)) (dy/dx)=2y.  
	A.	Solve the differential equation.
	B.	Consider the initial value problems 
			({{y(0)=1}, {y(0)=-1}, {y(2)=1}, {y(2)=-1}, {y(-2)=1}, {y(-2)=-1}})
		Determine the intervals of existence.

Using Mathematica to solve DEs....

Simply use the DSolve[ ] command...

Consider the differential equation (dy/dx)=3 y.
	A.	Solve the differential equation.
	B.	Solve the initial value problem when y(0)=1.
	C.	Plot the solution of the IVP above.

sol1=DSolve[{y'[x]==3 y[x]}, y[x], x]

To obtain a function, you need to use the replace all command...Notice the constant is included with this function.

y[x] /. sol1

ReplaceAll[y[x],{{y[x]->(E^(3 x)) C[1]}}]

sol2=DSolve[{y'[x]==3 y[x], y[0]==1}, y[x], x]

Plot[y[x] /. sol2, {x, -1,1}, PlotStyle->{Blue, Thick}]

Solve the differential equation (dy/dx)=-(x/y), where y(2)=-1  Illustrate several solution curves.

Notice that we are warned that there may be more than one "part" of the solution.

DSolve[{y'[x]==-x/y[x]}, y[x], x]

Go get the top part...

sol3A=DSolve[{y'[x]==-x/y[x], y[2]==1}, y[x], x]

Now for the bottom part...

sol3B=DSolve[{y'[x]==-x/y[x], y[2]==-1}, y[x], x]

Now, plot them both....

Plot[{y[x] /. sol3A, y[x] /. sol3B}, {x, -3, 3}, 
PlotStyle->{{Blue, Thick}, {Red, Thick}},
AspectRatio->Automatic]

The Mixing Problem

These mysterious problems are easy to model if you just hold your mouth right...the solution is:
	Rate of Change = Rate In - Rate Out.
	
You sometimes model the rate of change of an amount = concentration × volume.

You sometimes model the rate of change of a concentration = (amount/volume).

Example

Consider a large 1200 gallon tank initially contains 600 gallons of water with 7 pounds of salt dissolved in it.  Water enters the tank at a rate 10 gallons per minute with a concentration of 0.3 pounds of salt per gallon from pipe #1, and brine enters from pipe #2 at a rate of 12 gallons per minute with 0.5 pounds of salt per gallon.  The tank is always perfectly mixed and a third pipe ensures that the number of gallons never changes.  
	A.	Find an expression for the total amount of salt in pounds.
	B.	Find an expression for the concentration of salt in pounds per gallon.

Homework - For Software, Use Mathematica.

Pages:  	33 - 36
Problems:	1, 5, 7, 9, 11, 13, 20, 23, 24, 25, 27, 33, 39, 40

Slope Fields

Basic Definitions

Consider the differential equation (dy/dx)=f(x,y).  This represents the slope through the point (x, y).

A direction (or slope) field is generated by drawing a short line segment through each point (x, y) having slope given by 
f(x,y).

Direction Fields with Mathematica

In Mathematica, you can use the VectorPlot command to generate direction fields.  The basic definition of the command is VectorPlot[{(Subscript[v, x]), (Subscript[v, y])}, {x, (Subscript[x, min]), (Subscript[x, max])}, {y, (Subscript[y, min]), (Subscript[y, max])}]  A bit of explanation is in order.  We are creating a plot of a vector field {(Subscript[v, x]), (Subscript[v, y])}, where (Subscript[v, x]) and (Subscript[v, y]) are both functions of x and y!  You can think of each vector being plotted as having a length of (Subscript[v, x]) in the x direction and a length of (Subscript[v, y]) in the y direction.  

Since we only have one differential equation, we assume that (Subscript[v, x])=1.  So, for the typical differential equation (dy/dx)=f(x,y), we have (Subscript[v, x])=1 (at least for a while...) and (Subscript[v, y])=f(x,y).

Thus, the typical command will be VectorPlot[{1, f(x,y)}, {x, (Subscript[x, min]), (Subscript[x, max])}, {y, (Subscript[y, min]), (Subscript[y, max])}].

Recall that the creation of direction fields is a qualitative analysis.  Finding the solution of a differential equation is an analytic process.  It is really nice to be able to find solutions to differential equations, however, it is rarely possible!

Example

Create a slope field for the differential equation (dy/dx)=k y,  k=2, 0.5, -1, -3.

VectorPlot[
{1, 2 y}, {x, -4, 4}, {y, -4, 4},
FrameLabel->{x,y},
Axes->True,
VectorScale->{Small, Small, None},
PlotLabel->"Direction Field"
]

Here is the "cute" way to create these plots quickly...use the Table command.

Table[
VectorPlot[
{1,k y}, {x, -4, 4}, {y, -4, 4},
FrameLabel->{x,y},
Axes->True,
VectorScale->{Small, Small, None},
PlotLabel->StringForm["Direction Field, k = ``.", k],
ImageSize->300
],
{k, {2, 0.5, -1, -3}}
]

We can also add sample solution curves to a direction field...Mathematica calls these StreamPoints

VectorPlot[
{1, 2 y}, {x, -4, 4}, {y, -4, 4},
AxesLabel->{x,y},
Axes->True,
VectorScale->{Small, Small, None},
VectorStyle->Gray,
StreamPoints->Coarse,
StreamStyle->{Blue, Thick, "Line"},
StreamScale->Full,
PlotLabel->"Direction Field with Solution Curves"
]

If you wish to have particular solution curves through given initial conditions, then...

points={{0,1}, {0,-1}, {-2, 1}, {2, 1}, {2, -1}};
VectorPlot[
{1, 2 y}, {x, -4, 4}, {y, -4, 4},
FrameLabel->{x,y},
Axes->True,
VectorScale->{Small, Small, None},
VectorStyle->Gray,
StreamPoints->points,
StreamStyle->{Blue, Thick, "Line"},
StreamScale->Full,
PlotLabel->"Direction Field Using Particular Solution Curves",
Epilog->{PointSize[Large],Red, Point[points]}
]

Special Cases:  Slope Fields for Autonomous ((dy/dt)=f(y)) and Equations of the form (dy/dt)=f(t).

For differential equations of the form (dy/dt)=f(t), the slope field does not change in the vertical direction.

For autonomous differential equations, i.e. (dy/dt)=f(y), the slope field does not change in the horizontal direction.  As you will quickly notice, you get translated copies of a single curve, so it is easy to construct infinitely many solutions just by shifting one graph left or right.

Integrate[1/(y^2+1), y]

Example

A.	Solve and then create a slope field for (dy/dx)=(x^2)+1.
B.	Solve and then create a slope field for (dy/dx)=(y^2)+1.

GraphicsGrid[{
{VectorPlot[
{1, x^2+1}, {x, -4, 4}, {y, -4, 4},
FrameLabel->{x,y},
Axes->True,
VectorScale->{Small, Small, None},
VectorStyle->Gray,
StreamPoints->Coarse,
StreamStyle->{Blue, Thick, "Line"},
StreamScale->Full,
PlotLabel->Style["Direction Field, (dy/dx) = (x^2)+1", Bold, 16]
],
VectorPlot[
{1, y^2+1}, {x, -4, 4}, {y, -4, 4},
FrameLabel->{x,y},
Axes->True,
VectorScale->{Small, Small, None},
VectorStyle->Gray,
StreamPoints->Coarse,
StreamStyle->{Blue, Thick, "Line"},
StreamScale->Full,
PlotLabel->Style["Direction Field, (dy/dx) = (y^2)+1", Bold, 16]
]}
},
ImageSize->600]

Homework - For Software, Use Mathematica.

Pages:	 47 - 51
Problems:	   1, 2, 4, 6, 7, 8, 11, 12, 13, 15, 16

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])).

Hene, 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:  Case h = 1

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

(* The same basic code will be used over and over again...the only change is the step size, h.  So only this segment of code will be commented. *)
beginPoint=0;  (* Starting point of the approximation *)
endPoint=4;  (* Ending point of the approximation *)
h=1;   (* Step size used in Euler's Method *)
x0=0;  (* Initial x condition *)
y0=1;   (* Initial y condition *)
MaxN=(endPoint-beginPoint)/h+1;  (* Formula for computing the number of steps *)
xList={x0};  (* Initialize our list of x-values *)
yList={y0};  (* Initialize our list of y-values *)
xyList1={{x0, y0}};  (* Initialize our lis of ordered pairs *)

(* Here we loop from 1 to MaxN - 1 performing Euler's method and appending the results in the lists *)
For[n=1,n<MaxN,n++,
xList=Append[xList,xList[[n]]+h];
yList=Append[yList,yList[[n]]+h*f[xList[[n]], yList[[n]]]];
xyList1=Append[xyList1, {xList[[n]]+h, yList[[n]]+h*f[xList[[n]], yList[[n]]]}]
]


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

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

N[xyList1]

Here is the Code to Perform Euler’s Method:  Case h = 0.5

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

(* The same basic code will be used over and over again...the only change is the step size, h.  So only this segment of code will be commented. *)
beginPoint=0;  (* Starting point of the approximation *)
endPoint=4;  (* Ending point of the approximation *)
h=0.5;   (* Step size used in Euler's Method *)
x0=0;  (* Initial x condition *)
y0=1;   (* Initial y condition *)
MaxN=(endPoint-beginPoint)/h+1;  (* Formula for computing the number of steps *)
xList={x0};  (* Initialize our list of x-values *)
yList={y0};  (* Initialize our list of y-values *)
xyList1={{x0, y0}};  (* Initialize our lis of ordered pairs *)

(* Here we loop from 1 to MaxN - 1 performing Euler's method and appending the results in the lists *)
For[n=1,n<MaxN,n++,
xList=Append[xList,xList[[n]]+h];
yList=Append[yList,yList[[n]]+h*f[xList[[n]], yList[[n]]]];
xyList1=Append[xyList1, {xList[[n]]+h, yList[[n]]+h*f[xList[[n]], yList[[n]]]}]
]

eulerHisOneHalf=(* Note that you can save plots for later! *)
ListPlot[xyList1,
Joined->True,
Mesh->All,PlotMarkers->{Automatic,Medium},
PlotStyle->{Red,Thick},
AxesOrigin->{0,0}];

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

Several Euler’s Method Plots Compared

(* This is the legend for the plot *)
Needs["PlotLegends`"];
myLegend={
{Graphics[{Blue, Thick, Line[{{0.05,.7},{0.75,.7}}]}], "h = 1"},
{Graphics[{Red, Thick, Line[{{0.05,.7},{0.75,.7}}]}], "h = 0.5"}
};

ShowLegend[ (* This places a legend on the graphical output of the Show command *)
Show[
eulerHisOne, eulerHisOneHalf,
PlotRange->{-0.5,2.6},
AxesLabel->{x, y},
PlotLabel->Style[Text["Euler's Method Comparison"], 16]
],
{myLegend, LegendPosition->{0.2, -0.3}, LegendSize->{0.62, 0.3}}
]

Direction Field for this Problem

f[x_, y_]=Exp[-x^2];  

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->{Darker[Green], Thick, "Line"},
StreamScale->Full
];  (* 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.

GraphicsGrid[{
{
Show[dField,
PlotLabel->Style["Direction Field, (dy/dx) = (e^(-(x^2)))", Bold, 16]],
Show[eulerHisOne, PlotRange->{-2,2.6},
PlotLabel->Style["Euler's Method, h = 1", Bold, 16]]
},
{
Show[eulerHisOneHalf, PlotRange->{-2,2.6},
PlotLabel->Style["Euler's Method, h = 0.5", Bold, 16]],
Show[dField,eulerHisOneHalf,eulerHisOne,
PlotLabel->Style["Direction Field  and Euler's Method", Bold, 16]
]
}
},
ImageSize->600]

Errors in Numerical Methods

We are approximating solutions, so there can be error.  Consider this example

DSolve[y'[x]==Exp[x] Sin[y[x]], y[x],x]

Integrate[Csc[y],y]

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

(* The same basic code will be used over and over again...the only change is the step size, h.  So only this segment of code will be commented. *)
beginPoint=0;  (* Starting point of the approximation *)
endPoint=5.3;  (* Ending point of the approximation *)
h=0.1;   (* Step size used in Euler's Method *)
x0=0;  (* Initial x condition *)
y0=5;   (* Initial y condition *)
MaxN=(endPoint-beginPoint)/h+1;  (* Formula for computing the number of steps *)
xList={x0};  (* Initialize our list of x-values *)
yList={y0};  (* Initialize our list of y-values *)
xyList1={{x0, y0}};  (* Initialize our lis of ordered pairs *)

(* Here we loop from 1 to MaxN - 1 performing Euler's method and appending the results in the lists *)
For[n=1,n<MaxN,n++,
xList=Append[xList,xList[[n]]+h];
yList=Append[yList,yList[[n]]+h*f[xList[[n]], yList[[n]]]];
xyList1=Append[xyList1, {xList[[n]]+h, yList[[n]]+h*f[xList[[n]], yList[[n]]]}]
]

oopsie1=(* Note that you can save plots for later! *)
ListPlot[xyList1,
Joined->True,
Mesh->All,PlotMarkers->None,
PlotStyle->{Blue,Thick},
AxesOrigin->{0,0}];

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

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

(* The same basic code will be used over and over again...the only change is the step size, h.  So only this segment of code will be commented. *)
beginPoint=0;  (* Starting point of the approximation *)
endPoint=5.3;  (* Ending point of the approximation *)
h=0.05;   (* Step size used in Euler's Method *)
x0=0;  (* Initial x condition *)
y0=5;   (* Initial y condition *)
MaxN=(endPoint-beginPoint)/h+1;  (* Formula for computing the number of steps *)
xList={x0};  (* Initialize our list of x-values *)
yList={y0};  (* Initialize our list of y-values *)
xyList1={{x0, y0}};  (* Initialize our lis of ordered pairs *)

(* Here we loop from 1 to MaxN - 1 performing Euler's method and appending the results in the lists *)
For[n=1,n<MaxN,n++,
xList=Append[xList,xList[[n]]+h];
yList=Append[yList,yList[[n]]+h*f[xList[[n]], yList[[n]]]];
xyList1=Append[xyList1, {xList[[n]]+h, yList[[n]]+h*f[xList[[n]], yList[[n]]]}]
]

oopsie2=(* Note that you can save plots for later! *)
ListPlot[xyList1,
Joined->True,
Mesh->All,PlotMarkers->None,
PlotStyle->{Red,Thick},
AxesOrigin->{0,0}];

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

Look how closely these agree...and then suddenly, not so much!

(* This is the legend for the plot *)
Needs["PlotLegends`"];
myLegend={
{Graphics[{Blue, Thick, Line[{{0.05,.7},{0.75,.7}}]}], "h = 0.1"},
{Graphics[{Red, Thick, Line[{{0.05,.7},{0.75,.7}}]}], "h = 0.05"}
};

ShowLegend[ (* This places a legend on the graphical output of the Show command *)
Show[
oopsie1, oopsie2,
PlotRange->{-0.5,5.2},
AxesLabel->{x, y},
PlotLabel->Style[Text["Euler's Method Comparison"], 16]
],
{myLegend, LegendPosition->{-0.5, -0.4}, LegendSize->{0.62, 0.3}}
]

Homework - For Software, Use Mathematica.  You need to learn some programming, so use the code I have provided!  In particular, stay away from DSolve and NDSolve.

Pages:		61 - 63
Problems:	    1, 2, 5, 6, 11