Eulers method to approximate gravity in one dimension

Question

I am going to program Eulers method in Octave to approximate gravity in 1-dimension. I understand the formula for Eulers method, which is equal to:

What I don't understand is what my function $f(t,y)$ is in this case. What do I have to insert into the formula to get my next $y$-point?

Answer 1

Depending on the scale, you have $$ \ddot x=a(t,x)=-g $$ or $$ \ddot x=a(t,x)=-\frac{GM}{(R+x)^2} $$ where $G$ is the gravitational constant, $M$ the mass of Earth and $R$ its radius.

As this is second order, a first order system would have the form $$ \binom{\dot x}{\dot v}=f(t,(x,v))=\binom{v}{a(t,x)}. $$

Answer 2

Euler's method is usually used to discretize systems of explicit first order ordinary differential equations. Such a system is given by $$y'(x) = f(x,y(x))$$ where $f$ is in general a vector-valued, continuous function and we are looking for the function $y(x)$. In general, enough initial conditions are provided and we call this an initial value problem. In classical mechanics, the systems originate from Newton's second law, that the sum of the acting forces is equal to mass times acceleration of your object in consideration. In your case we have $$-mg = ma = mx''(t) \Leftrightarrow x''(t) = -g$$ However this is a second order equation and we must therefore write it as a system os first order to solve it using Euler. We get $$\begin{bmatrix} x_1(t)\\x_2(t)\end{bmatrix}'= \begin{bmatrix} x_2(t)\\-g\end{bmatrix}$$ so we have $$f(t,x(t)) = \begin{bmatrix} x_2(t)\\-g\end{bmatrix}$$ Here is an example code:

y0 = [10;-1];
t0 = 0;
tN = 1;
N = 10;
g = 9.81;
f = @(t,x) [x(2);-g];
[t,x] = euler(f,t0,tN,y0,N);
plot(t,x(1,:));
grid on;

using this implementation of the euler method:

function [ t,y ] = euler( f,t0,tN,y0,N )
h = (tN - t0)/N;
t = t0:h:tN;
y = zeros(length(y0),N+1);
y(:,1) = y0;
for k = 1:N
    y(:,k+1) = y(:,k) + h * f(t(k),y(:,k));
end
end

implemented in MATLAB. I used initial velocity $-1$ and initial place $10$. It produces:

Analytical Solution. Consider the IVP

$$\begin{bmatrix} x_1(t)\\x_2(t)\end{bmatrix}'= \begin{bmatrix} x_2(t)\\-g\end{bmatrix} = \begin{bmatrix} x_2(t)\\0\end{bmatrix}+\begin{bmatrix} 0\\-g\end{bmatrix} \qquad x(0) = \begin{bmatrix} x_0\\v_0\end{bmatrix}$$ From the last coordinate it is obvious, that $$x_2(t) = -gt + c$$ for some real constant $c$. Furthermore, from the first coordinate we get $$x'_1(t) = -gt + c \Leftrightarrow x_1(t) = -\frac{1}{2}gt^2 + ct + d$$ for some real constant $d$. by the initial conditions, we get the linear system $$x_0 = d \qquad v_0 = c$$ Hence our solution to the IVP is the vector-valued function $$x(t) = \begin{bmatrix} -\frac{1}{2}gt^2 + v_0t + x_0\\-gt + v_0 \end{bmatrix}$$ where the first coordinate corresponds to the place and the second to the velocity. Adding this in the plot above yields

Where red is the euler method, blue is the analytical solution and black my implementation of the step-size controlled DOPRI5 algorithm for reference (which can be found here).