March 20, 2012

How to solve ODEs using octave.

The odepkg package for octave is a great tool to solve ordinary differential equations (ODEs). To solve some problem we only have to bring it into a form octave understands. We can do that following four steps:

  1. determine the variables of the problem
  2. figure out which variables change and how they do it
  3. if higher order derivatives do appear, reduce the order of the system
  4. write a function that can gives the first derivatives for all values considered and solve it with octave's ode solvers
First example: Shooting a paperball using a slingshot (no friction).

problems like the trajectory of a sphere fired by a slingshot can be solved using ODEs

This is a very basic example. In fact, we can give the analytical solution (parabula). The paperball is accelerated by the catapult and flies some distance until it hits the groud. 
1. The variables of the problem are
  • the position of the paperball in space r(t) = ( x(t), y(t) )
  • the velocity    v(t) = ( vx(t), vy(t) )
  • the acceleration a(t) = (ax(t), ay(t) )
Here we use the fact that we can split each of the variables of motion into two independent parts, showing in the x or and in the y-direction.

2. Now we have to consider how our variables change:
  • the change in position is proportional to the velocity at a certain moment:   dr/dt = (vx, vy)
  • the velocity is changed by the acceleration: dv/dt = (ax, ay)
  • for our case, the acceleration is constant. We consider the gravitational force: F=m a where a = g = 9.81 m/s^2 is the gravitational acceleration. It only acts on the y-component of the velocity (perpendicular to the ground)
As we figured out which variables change in what way, we consider how our result should look like. Octave expects the problem to form a vector (1xN components). What we want to know is the position of the paperball and its velocities for all times t.

The result should be of the form   result(t) = [rx, ry, vx, vy]

What we now have to do is to give octave a vector containing the derivative of our result vector. We do this by defining a function.  [drx, dry, dvx, dvy] = deriv_result( time,position)

3. if we did not already do so, we have to formulate the problem as a system of first order ordinary differential equations.

What does this mean? Instead of writing the second derivative of some value, we define some additional variable that is  the first derivative of the initial value and consider its changes.

example: 
  • the acceleration is the second derivative of the position  d^2r / dt^2 = a
  • we introduce the velocity. it is the first derivative of the position dr/dt = v
  • the change of the velocity is given by the acceleration dv/dt = a
  • now we have to solve the problem both for the position and the velocity and such the second order ODE has been reduced to two first order ODE
4. solving the problem in octave

We define derivative function

function dr = dr_gravi(t,r)
  g  = 9.81; %m/s^2
  vx = r(3);
  vy = r(4);
  ax = 0;
  ay = -g;  
  dr = [vx,vy,ax,ay];
endfunction

In the program we have to set some initial value. For our problem this is the initial position r0=(x0,y0) and the initial speed v0 = (vx0,vy0)

X0 = 0
Y0 = 0
VX0 = 2
VY0 = 3

initialVector = [ X0, Y0, VX0, VY0]

Also,  some integration window is necessary: for our problem, this is the starting time and the ending time for time integration:

StartT= 0 %s
StopT = 0.7 %s

then we let octave do the job:

[T,Result] = ode45( @dr_gravi,[StartT, StopT], initialVector)

Here we used the ode45  solver, but there are many more. See the odepkg documentation for details.
Our result looks like the following

T =

   0.00000
   0.07000
   0.14000
   0.21000
   0.28000
   0.35000
   0.42000
   0.49000
   0.56000
   0.63000
   0.70000
   0.70000

Result =

   0.00000   0.00000   2.00000   3.00000
   0.14000   0.18597   2.00000   2.31330
   0.28000   0.32386   2.00000   1.62660
   0.42000   0.41369   2.00000   0.93990
   0.56000   0.45545   2.00000   0.25320
   0.70000   0.44914   2.00000  -0.43350
   0.84000   0.39476   2.00000  -1.12020
   0.98000   0.29231   2.00000  -1.80690
   1.12000   0.14179   2.00000  -2.49360
   1.26000  -0.05679   2.00000  -3.18030
   1.40000  -0.30345   2.00000  -3.86700
   1.40000  -0.30345   2.00000  -3.86700

T is the vector that gives us the time in seconds for some row. Results contains the position and velocity components for the times stored in T.  Using this data, we can plot the trajectory and other dynamic information.

The following code plots the trajectory of the ball and the time dependence of the velocity components:

figure(1)
  plot(Result(:,1),Result(:,2),'o');
  xlabel 'x position / m'
  ylabel 'y position / m'

figure(2)
  plot(T, Result(:,3),'.',
     T, Result(:,4),'x');
  xlabel 'time / s'
  ylabel 'velocity / (m /s)'
  legend('vx','vy')



Is there a way to influence the calculation? Of course there is. We can do that using odeoptions

options = odeset( 'RelTol',1e-4,
'AbsTol',1e-4,
'InitialStep',StopT/1e3,
'MaxStep',StopT/1e3)
[T,Result] = ode45( @dr_gravi,[StartT, StopT], initialVector , options )

There are many options to play around with, of which most are self-explanatory . You have to find a tradeoff between accuracy and calculation time. Again, see the odepkg documentation.

A more complex example: slingshot with friction.
Additional to the gravity, let's consider the drag force influencing the flight of our paperball. It is a force proportional to the velocity  v of an object.

F = -1/2 rho Cd A v^2

Additional parameters appear, which we set here: 
  • rho, the density of air ( rho = 1.2 kg / m^3)
  • Cd - the drag coefficient . We consider a sphere here, which has Cd = 0.45
  • A is the crosssection of the sphere
Again, we can split the force into two components: Fx and Fy. To calculate the acceleration components, we also need the mass of our object, as F=ma which leads to a = F/m.

We can give the derivative function as follows with additional parameters (area and mass of the sphere):


function dr = dr_gravi_friction(t,r,area,mass)
  g = 9.81; %m / s^2
  cdSphere = 0.45;
  rhoAir = 1.20; %kg / m^3
  frictioncoefficient = 1/2 * rhoAir * cdSphere * area / mass;  
  
  vx = r(3);
  vy = r(4);
  ax = - frictioncoefficient * vx^2;          % only friction

  ay = -( g + frictioncoefficient * vy^2 ) ;  % friction and  
                                              %  gravitation
  dr = [vx,vy,ax,ay];
endfunction

Now we can try out how friction affects the trajectory.  We calculate this for a relative heavy sphere, a light sphere and without friction. The radius of the spheres shall be r=2 cm


mass1 = 0.1 %kg
mass2 = 0.001 %kg
area1 = pi * 0.02^2   %m^2
area2 = area1;
options = odeset( 'RelTol',1e-4,
  'AbsTol',1e-4,
  'InitialStep',StopT/1e3,
  'MaxStep',StopT/1e3)

% trajectory of heavy sphere
[T1,Result1] = ode45( @dr_gravi_friction,[StartT, StopT], \
    initialVector , options, area1, mass1  )

% trajectory of light sphere
[T2,Result2] = ode45( @dr_gravi_friction,[StartT, StopT], \
    initialVector , options,area2, mass2 )
% trajectory without friction
[T3,Result3] = ode45( @dr_gravi,[StartT, StopT], \
    initialVector , options)

When we look at the results, we see the following: The relative impact of the friction force on the light sphere is quite strong. It considerably gets slowed down and hits the ground earlier. The trajectory deviates more and more from the parabula shape while propagating. Contrarily, the heavy sphere almost acts as if  friction was absent.  This basically is what we expected.

trajectories of spheres

2 comments:

  1. I am trying to solve an ODE with octave using lsode; I can do for the simple cases and I know how to setup the .m file. But I have a case in which the function need to be changed. For example, I want to reset some variables to zero everytime it crosses one. Once this happens a fixed number of time, I want the lsode to solve a different set of equations.

    In case of a chemical reaction, we have propagation phase but once the reactants are over, we need to go to the termination phase. I am trying to simulate something like that.

    Any help will be appreciated.

    CKM

    ReplyDelete
  2. Well basically there is a way to do something like that. You are free to include if - statements in your derivative function code, e.g

    if [condition]
    rhs = x1
    else
    rhs = x2
    endif

    (when rhs is the derivative you want to integrate ...)

    ReplyDelete