Euler's Method 
for a single equation

The basic idea here is to create an approximation to the solution of the single differential equation

y'[t]==f[t,y].

at a particular initial value

y[(Subscript[t, 0])]==(Subscript[y, 0.])

We'll start with a simple code that uses a given step size h and makes k iterations.  The form of a “For” loop in Mathematica is 

For[start, test, increment, body]  

Each of the four arguments of the “For” loop may contain multiple statements separated by semicolon (:).

We can easily make the step size and number of iterations arguments of our code, but for now let’s define them externally:

h=0.1;
k=4;

We will also use a specific example function f for the right side  along with initial values so we can test our code.

f[tvariable_,yvariable_]=yvariable^2+tvariable^2;

tzero=1;
yzero = 1;

What does the existence and uniqueness theorem tell you about the ODE with this f as the right hand side?  Do you know how to solve this ODE?  Is it separable, linear, exact (anything you know how to solve)?

Let's create variables in which to store the time values and the dependent variable.

Tvalues=Array[T,k+1,0]

Yvalues=Array[Y,k+1,0]

Let's make a preliminary loop to assign the T values, just to make sure we know how to do a “For” loop and to simplify the main code.  (We can easily include this step in the main “For” loop later.

For[   j=0;T[0]=tzero,
	j<k,
	j++,
	T[j+1]=T[j]+h
       ]

Tvalues

(That seems to have worked.)

For[   j=0; Y[0]=yzero, 
	j<k,
	j++,
	Y[j+1]=Y[j]+h f[T[j],Y[j]]
       ]

Yvalues

output=Table[{T[j],Y[j]},{j,0,k}]

outputplot=Show[ListPlot[output,PlotStyle->{Black,AbsolutePointSize[7]},AxesOrigin->{0,0}],ListPlot[output,PlotStyle->Black,Joined->True]]

Now, I’d like to compare the output of our simple Euler method code with something that represents the “actual solution.”  In this case, I mean the result of a better numerical approximation.

soln=NDSolve[{odey'[odet]==f[odet,odey[odet]],odey[tzero]== yzero},{odey},{odet,1,1+k h}]

y[t_]=odey[t]/.soln[[1]]

Show[outputplot,Plot[y[t],{t,1,1+k h},PlotStyle->Green]]

startplot=Show[outputplot,Plot[y[t],{t,1,1+k h},PlotStyle->Green],PlotRange->All]

Here is a more general code with some arguments.

firsteuler[stepsize_,iterations_]:=For[
	j=0; Clear[T,Y,Tvalues,Yvalues];
	Tvalues=Array[T,iterations+1,0];
	Yvalues=Array[Y,iterations+1,0];
	output=Table[{T[j],Y[j]},{j,0,iterations}];
	T[0]=tzero;Y[0]=yzero, 
	j<iterations,
	j++,
	T[j+1]=T[j]+stepsize;
	Y[j+1]=Y[j]+stepsize f[T[j],Y[j]]
       ]

firsteuler[0.1,4]

ListPlot[output]

firsteuler[0.4/10,10]

secondoutputplot=Show[ListPlot[output,PlotStyle->{Black,AbsolutePointSize[6]},AxesOrigin->{0,0}],ListPlot[output,PlotStyle->Black,Joined->True]]

Show[startplot,secondoutputplot]

Here is a manual test designed to bring out the third important point about numerical methods.

test[n_]:=firsteuler[0.4/n,n]
testplot[n_]:=Show[startplot,ListPlot[output,PlotStyle->Black,Joined->True]]

test[100]
testplot[100]

test[20]
testplot[20]

test[1000000]
testplot[1000000]

Lesson 1:  It’s takes a lot of time to do many iterations, and it takes a lot of memory.

ls -l simpleeuler.nb
15633286 Feb 25 11:21 simpleeuler.nb

Let’s mitigate the memory problem.

streamlinedfirsteuler[stepsize_,iterations_]:=For[
	j=0; Clear[finalvalue];
	time=tzero;val=yzero, 
	j<iterations,
	j++,
	newtime=time+stepsize;
	newval=val+stepsize f[time,val];
	time=newtime;
	val=newval
       ]

streamlinedfirsteuler[0.4/10^6,10^6]

val
y[1.4]

streamlinedfirsteuler[0.4/10^7,10^7]

val
y[1.4]

10^10 seems to be too much.