Solving differential equations numerically with octave.

Object in free fall

Octave is a great tool for solving differential equations. Here I will try to give a simple example of doing so, by showing how to calculate the position and velocity of an object in free fall. I will start with the analytical solution, and move forward to the numerical solution using octave.

Analytical solution

Newtonian physics dictates that, when one object falls freely, its acceleration is g. With that in mind, finding its velocity and position is just a matter of solving two differential equations. These are

\frac{dv}{dt}=a and \frac{d^2v}{dt^2}=a

Solving them analytically, we find

v(t)=g\cdot t+v_0=-10\cdot t
y(t)=\frac{1}{2}gt^2+v_0t+x_0=-\frac{1}{2}10t^2

Substituting with y_0=25 ~m, g=-10 and  v_0=0 ~m/s we end up with

v(t)=-10\cdot t
y(t)=-\frac{1}{2}10t^2

Now lets make some plots in octave, and see how this solution looks like. To do that, we start by defining all the values of interest

y0=25;
v0=0;
g=-10;
t=0:.1:2;

After that, and to make the code more clear, we define y(t) and v(t) as anonymous functions. To find x for a specific time t, we can now use x_t(t)

%define the equations as anonymous functions
y_t= @(t) y0+.5*g*t.**2;
v_t= @(t) g*t;

Now it is time to start plotting the variables of interest. We are going to do three plots. First we plot position, as a function of time

%plot y(t)
subplot(1,3,1)
plot(t ,y_t(t))
xlabel("t (sec)")
ylabel("y (m)")
box()

followed by velocity as a function of time

%plot  v(t)
subplot(1,3,2)
plot(t,v_t(t))
xlabel("t (sec)")
ylabel("v (m/sec)")
box()

and finally, we plot x and y position of the object every one tenth of a second — notice that t was first defined in line 4. To do so, we will need an x vector, equal in size with the y vector. The easiest way to produce it to multiply y_t(t) with zero.

%plot y(x) discretely
subplot(1,3,3)
plot(y_t(t)*0,y_t(t),"b.","markersize",15)
xlabel("x(m)")
ylabel("y(m)")
box()

y0=25;
v0=0;
g=-10;
t=0:.1:2;

%define the equations as anonymous functions
y_t= @(t) y0+.5*g*t.**2;
v_t= @(t) g*t;

%plot y(t)
subplot(1,3,1)
plot(t,y_t(t))
xlabel("t (sec)")
ylabel("y (m)")
box()

%plot v(t)
subplot(1,3,2)
plot(t,v_t(t))
xlabel("t (sec)")
ylabel("v (m/sec)")
box()

%plot y(x) discretely
subplot(1,3,3)
plot(y_t(t)*0,y_t(t),"b.","markersize",15)
xlabel("x(m)")
ylabel("y(m)")
box()

All these plots can be seen in the figure bellow

Analytically derived position and velocity for an object in free fall

Left: Position of object with respect to time Middle: Velocity with respect to time Right: Position of the object every 0.1 sec.

Numerical solution

Now lets try to produce the same results using a numerical solution. We can start again by defining all the values of interest. Only this time we can’t have acceleration defined here.

y0=25;
v0=0;
t=0:.1:2;

Octave can solve equations numerically, but for the solver to be able to call them at once, they should all be added in a single function. This function should have the form xdot=f(x,t), where xdot stands for dx/dt, and thus have the same size as x. For our particular case, we want to define a function where x(1) will represent position, and x(2) will represent velocity. Lets name this vector yv, and its derivative dyv. Notice how we were forced to put the acceleration inside the function (line 9).

%write all the  differential equations
%in a single function
function dyv=dyv(yv,t)
  dyv(1)=yv(2); %dy/dt=v
  dyv(2)=-10;   %dv/dt=-10
endfunction

Now we can call lsode, to compute the solution. Let us first reduce the absolute tolerance though, to get the exact same results as previously.

%change the tolerance to get better results
lsode_options("absolute tolerance",1e-30);
%compute the solution
yv=lsode("dyv",[y0 v0],t);

lsode needs the function name, the initial state vector, and the time to solve for. Time t should be a vector, and the first position of this vector — t(1) — should correspond to the time where the system is in the initial state. Here t was defined in line 3 to be t=0:.1:2. The solution — yv — is a matrix with two columns, corresponding to y and v, and with as many rows as the size of t, or in this time 21. Having solved our system, we can proceed in plotting position and velocity as before.

%plot y(t)
subplot(1,3,1)
plot(t,yv(:,1))
xlabel("t (sec)")
ylabel("y (m)")
box()

%plot v(t)
subplot(1,5,3)
plot(t,yv(:,2))
xlabel("t (sec)")
ylabel("v (m/sec)")
box()

%plot y(x) discretely
subplot(1,5,5)
plot(yv(:,1)*0,yv(:,1),"b.","markersize",15);
xlabel("x(m)")
ylabel("y(m)")
axis([-1 1])
box()

 

y0=25;
v0=0;
t=0:.1:2;

%write all the  differential equations
%in a single function
function dyv=dyv(yv,t)
  dyv(1)=yv(2); %dy/dt=v
  dyv(2)=-10;   %dv/dt=-10
endfunction

%change the tolerance to get better results
lsode_options("absolute tolerance",1e-30);
%compute the solution
yv=lsode("dyv",[y0 v0],t);

%plot y(t)
subplot(1,3,1)
plot(t,yv(:,1))
xlabel("t (sec)")
ylabel("y (m)")
box()

%plot v(t)
subplot(1,5,3)
plot(t,yv(:,2))
xlabel("t (sec)")
ylabel("v (m/sec)")
box()

%plot y(x) discretely
subplot(1,5,5)
plot(yv(:,1)*0,yv(:,1),"b.","markersize",15);
xlabel("x(m)")
ylabel("y(m)")
axis([-1 1])
box()

As expected, the results are the same as before.

Numerically computed position and velocity for an object in free fall

As above, only this time position and velocity were computed numerically

Code

You can find the full code of the examples in the text, under freefall_analytical.m and freefall_numerical.m. Saving to an image file requires a little bit of effort though. You can find the code I used to create these two plots bellow.

%%%%
% Analytical solution
%%%%
figure(1)

x0=25;
v0=0;
g=-10;
t=0:.1:2;

%define the equations as anonymous functions
y_t= @(t) x0+.5*g*t.**2;
v_t= @(t) g*t;

%plot y(t)
subplot(1,5,1)
plot(t,y_t(t),"linewidth",4)
box()
xlabel("t (sec)")
ylabel("y (m)")
set(gca,'xtick',[0 1 2])

%plot v(t)
subplot(1,5,3)
plot(t,v_t(t),"linewidth",4)
box()
xlabel("t (sec)")
ylabel("v (m/sec)")
set(gca,'xtick',[0 1 2])

%plot y(x) discretely
subplot(1,5,5)
plot(y_t(t)*0,y_t(t),"b.","markersize",15)
xlabel("x(m)")
ylabel("y(m)")
box()
axis([-1 1])
set(gca,'xtick',[0])

%save plots to a png file
print -dpng -color "-S600,600" freefallAnalytical.png

%%%%
% Numerical solution
%%%%
figure(2)
x0=25;
v0=0;
t=0:.1:2;

%write all the differential equations
%in a single function
function dyv=dyv(xv,t)
 dyv(1)=xv(2); %dy/dt=v
 dyv(2)=-10; %dv/dt=-10
endfunction

%change the tolerance to get better results
lsode_options("absolute tolerance",1e-30);
%compute the solution
yv=lsode("dyv",[x0 v0],t);

%plot y(t)
subplot(1,5,1)
plot(t,yv(:,1),"linewidth",4)
xlabel("t (sec)")
ylabel("y (m)")
set(gca,'xtick',[0 1 2])
box()

%plot v(t)
subplot(1,5,3)
plot(t,yv(:,2),"linewidth",4)
xlabel("t (sec)")
ylabel("v (m/sec)")
set(gca,'xtick',[0 1 2])
box()

%plot y(x) discretely
subplot(1,5,5)
plot(yv(:,1)*0,yv(:,1),"b.","markersize",15);
xlabel("x(m)")
ylabel("y(m)")
axis([-1 1])
set(gca,'xtick',[0])
box()

%save plots to a png file
print -dpng -color "-S600,600" freefallNumerical.png

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s