Euler's method for solving ODE's

Name:
Student #:

In this lab, we will solve the simple case of a falling object as a glimpse into numerical solutions of ODE's.

Newton's second law is an ODE

Almost all problems in mechanics can start with Newton's second law of motion, that is

∑F=m a.     	(1)

Equation (1) can be viewed as an ordinary differential equation (ODE). Let's first re-arrange it

a=((∑F)/m).		(2)

Now given that, by definition,

a=(dv/dt)=(((d^2)y)/(dt^2)),

we can write Equation (2) as a first order ODE for velocity v and a first order ODE for position i.e.

(dv/dt)=((∑F)/m),	(3)
	(dy/dt)=v

or as a second order ODE for position y, i.e.,

(((d^2)y)/(dt^2))=((∑F)/m).	(4)

For some problems, Equations (3) or (4), along with proper initial conditions, can be solved analytically.  These cases provide a basis for understanding the behaviour of numerical solutions, as we will see in what follows.

Falling object without air resistance

Consider the classical problem of an object of mass m thrown upward at t=0 with initial velocity V from a height h and falling without air resistance in a gravitational field g. Once the object is released, there is only one force that acts on the object: the force of gravity -mg (so ΣF=-mg in (2) for this situation). Therefore Equation (2) takes the following form,

a=((-mg)/m)=-g		(5)

If we are interested in finding the velocity v and the position y as functions of time t we can re-write (5) as

(dv/dt)=-g		(6)

and from the definition of velocity we can write

(dy/dt)=v			(7)

Equations (6) and (7) are two first-order ODEs with two unknowns v(t) and y(t). The problem requires initial conditions, which are:

y(0)=h		(8)

and

v(0)=V		(9)

Use Mathematica to solve analytically this set of ODEs for v and y. Solve the problem purely symbolically. This can be done with a single line of code.

Remove["Global`*"]

DSolve[{v'[t]==

Equations (6) and (7) can be discretized using Euler's method. As done in class, let's first discretize Equation (6) and solve numerically for v, i.e.,

(dv/dt)≈(((Subscript[v, n+1])-(Subscript[v, n]))/(Δt))=-g	(10)

where n is an index representing time and Δt is the time step that is determined by the person who writes the algorithm. The smaller time step used, the more accurate the solution will be when compared to the exact solution (generally speaking). 
Re-arranging (10) we get an equation for the velocity at a time Δt into the future, i.e. the next time step, represented by the index n+1,

(Subscript[v, n+1])=-g Δt+(Subscript[v, n])	(11)

In this case the initial condition corresponds to the first time index, i.e. n=1. The initial condition (9) is written as

(Subscript[v, 1])=V.		(12)

Now let’s discretize equation (7)

(dy/dt)≈(((Subscript[y, n+1])-(Subscript[y, n]))/(Δt))=(Subscript[v, n])	(13)

where (Subscript[v, n]) appears with the subscript n because it is the velocity corresponding to the same time as (Subscript[y, n]) .
Re-arranging (13) we get an equation for the position at a time Δt into the future, i.e. the next time step, represented by the index n+1,

(Subscript[y, n+1])=(Subscript[v, n]) Δt+(Subscript[y, n])	(14)

In this case the initial condition corresponds to the first time index, i.e. n=1. The initial condition (8) is written as

(Subscript[y, 1])=h.		(15)

This problem described by Equations (11), (12), (14) and (15) can now be solved numerically for y and v as long as the values of the parameters of the problems are provided. As in the case of numerical integration, an ODE cannot be solved symbolically with numerical algorithms. 
For the purpose of this example let's propose that:
	g=9.81 m/(s^2),
	V=30.0 m/s and
	h=1500 m.

For this problem we will need to do some programming and we will need to use the Do loop command. To learn about this command try the following command:
	Do[Print[(k^2)],{k,0,4}]

The Do[ ] function has two arguments separated by a comma:
     1.  The first argument (e.g. Print[(k^2)] ) is the program statement    
           which is calculated repeatedly.
     2.  The second argument (e.g. {k,0,4} ) specifies how many times the 
           expression in the first argument is computed and at what values of k.  
Another example is
	Do[k=(i^2); Print[k],{i,0,4,2}]  
Try this in the input box below

The iterator i takes the values i = 0, i = 2, and  i = 4 since the increment is Δi = 2.  Notice that there are multiple statements within the first argument of Do, and each statement is separated by a semicolon.

The Do command will be used to solve numerically Equations (10) and (13) for y and v with initial conditions (12) and (15).

Below is a code to do this where we have arbitrarily set the time step to Δt=1.0 s and the total time of the simulation to T=10 s.
	g = 9.81; 			(* Gravity *)
	V = 30.0; 			(* Initial velocity *)
	h = 1500;			(* Initial height *)
	T = 10; 			(* Simulation time *)
	Δt = 1.0;			(* Time step *)
	nmax = Round[T/Δt];	(* Max. number of iteration *)
	v[1] = V;			(* Initial velocity *)
	y[1] = h;			(* Initial position *)
	t[1] = 0;			(* Initial time *)
	Do[t[n+1] = n Δt; 
	v[n+1] = -g Δt + v[n]; 
	y[n+1] = v[n] Δt + y[n],
	{n, 1, nmax}]

Examine this code carefully. Try to understand it. Ask questions if you don't. It is important for the rest of this course that you understand this. Notice that in this case, the expression for y[n+1] requires v[n], so it we need to solve for v as we solve for y.

Now, code this here:

This should have created in memory the variables t, y and v which represent tabulated values that resulted from the instructions in the Do loop above. To see the values of these variables type ?t and ?v

If you are interested in plotting the velocity as a function of time, using the ListPlot command, you will need to create a list of pairs {t,v} with the command Table. To do this type the following command
	vdata = Table[{t[n], v[n]}, {n, 1, nmax+1}]

and you can now plot the data using ListPlot with axis labels and units. Do this here:

ListPlot[vdata

It would be illuminating to compare the numerical solution graphically with the exact analytical solution. The difficulty is that while we can plot the exact solution, which is an analytical function, with the command Plot, we can only plot the results of the numerical calculation with the command ListPlot. The two commands Plot and ListPlot are not compatible with each other. 

Fortunately, there is a way to combine graphics together with the command Show. Here's how this could be done.

First create a graphic of the exact position y(t) with the command Plot in the interval from t=0 to t=10 s using the following syntax:
	plot1 = Plot[(1/2)(2 h - g (t^2) + 2 t V ),{t,0,T},AxesLabel->{“t (s)”,”y (m)”}]

Then create the second graph of the numerical solution like this

ydata = Table[{t[n],y[n]},{n,1,nmax+1}];
	plot2 = ListPlot[ydata]

You can then combine these two graphics with the command
	Show[plot1,plot2,PlotRange->All];
Try that.

What do you think of the comparison? Does the result of the numerical algorithm provide an accurate solution compared to the exact solution? Briefly comment.

Comment:

Now redo the same comparative analysis between the exact and the numerical solution to this problem but this time use a time step of Δt=0.1 s instead of Δt=1.0 s. I recommend starting this new problem with a Remove["Global`*"]

Remove["Global`*"]

Finally, briefly comment the result of this analysis as compared to the case above with Δt=1.0 s.

Comment: