An Approach to Solving Ordinary Differential Equations (I.Urieli, Winter 2002)

In many real life applications where we wish to evaluate a function, we are only able to define the relationship that the derivative of that function must satisfy. One example is the dynamics of a vehicle in which we wish to determine the velocity or acceleration, however the unknowns are specified only in terms of ordinary differential equations (ODEs) and initial conditions. In this note we wish to show one approach to the numerical solution of ODEs using MATLAB and the Classical fourth order Runge-Kutta method. This is a very personal approach and is somewhat different to that normally used, however I will try to justify it further on. The approach will be developed and illustrated in terms of two case studies.

Case Study 1 - a large angle pendulum

Consider a simple pendulum having length L, mass m and instantaneous angular displacement q (theta), as shown below:

Fortunately we can always reduce a n-th order ODE into a n-vector set of first order ODEs. Thus we need only concentrate on developing numerical techniques for solving first order ODEs. Processing vectors by computer simply requires looping sequentially through the indices of the vector, and MATLAB is particularly geared to this type of analysis.

Consider now the practical aspects of writing a function for evaluating this set of derivatives. One problem is that vectors are usually assigned generic names such as y for the set of dependent variables and dy for their derivatives, however we would like to retain their identity in terms of q (theta) and w (omega). This can be done by means of indices, as shown in the MATLAB function below:

   function 
 
 [y,dy] = dpend(t,y) 
 
   %Pendulum differential equations 
 
 
   % G [m/s^2] acceleration due to gravity 
 
 
   % LENGTH [m] pendulum length 
 
 
   % Izzi Urieli - Jan 21, 2002 
 
 
   global 
 
 G LENGTH 
 
  THETA = 1; 
  % index of angle [radians] 
 
 
  OMEGA = 2; 
  % index of angular velocity [rads/sec] 
 
 
  dy(THETA) = y(OMEGA); 
 
  dy(OMEGA) = -(G/LENGTH)*sin(y(THETA)); 
 
  dy = [dy(THETA);dy(OMEGA)]; 
  % Column vector 
 
 

Notice the use of upper case characters for both the indices and the global variables, as is the convention in MATLAB. Notice also that I am returning both y and dy (instead of only dy as is normally done). This is because I find it convenient to load other parameters on these vectors - again I will justify this in more complex example following this case study.

Consider now a typical main function for this case study:

   %Solve the large angle pendulum problem 
 
 
   % Izzi Urieli - Jan 21, 2002 
 
 
  THETA = 1; 
  % index of angle [radians] 
 
 
  OMEGA = 2; 
  % index of angular velocity [rads/sec] 
 
 
   global 
 
 G LENGTH 
 
  G = 9.807; 
  % [m/s/s} accel. due to gravity 
 
 
  LENGTH = input( 
  'Enter pendulum length [m] >' 
 
 ); 
 
  angle0 = input( 
  'Enter initial pendulum angle [degrees] >' 
 
 ); 
 
  period = 2*pi*sqrt(LENGTH/G); 
 
  fprintf(' 
  small angle period is %.3f seconds\n' 
 
 ,period); 
 
  tfinal = 2*period; 
 
  n = input( 
  'Enter the number of solution points >' 
 
 ); 
 
  dt = tfinal/(n - 1); 
 
  t = 0; 
 
  y = [angle0*pi/180;0]; 
  % Initialise y as a Column vector 
 
 
  time(1) = t; 
 
  angle(1) = angle0; 
 
   for 
 
 i = 2:n 
 
  [t, y, dy] = rk4( 
 
   'dpend' 
 
 
  ,2,t,dt,y); 
 
  time(i) = t; 
 
  angle(i) = y(THETA)*180/pi; 
 
  fprintf('%10.3f%10.3f\n',time(i),angle(i)); 
 
   end 
 
 
  plot(time,angle) 
 

As with many MATLAB functions, this one is very short and mostly self-documenting. Notice that after specifying the system parameters we use the for loop to build the ' time ' and ' angle ' arrays in order to plot them. The most interesting aspect of this program is the ' rk4 ' function call, to integrate one step ' dt ' of the differential equation set given in function ' dpend ', as follows:

 [t, y, dy] = rk4( 
  'dpend' 
 
 ,2,t,dt,y); 

Thus the input arguments are the string constant representing the derivative function 'dpend', the number of differential equations in the set (2), the independent variable and increment ( t,dt ), and the dependent variable vector ( y ). The output arguments are the solution variables and derivatives ( t,y,dy ) integrated over one time step dt .

MATLAB provides five separate routines for solving ODEs of the form dy = f(x,y) where x is the independent variable and y and dy are the solution and derivative vectors. They are extremely simple to use, however they do not provide the versatility that I am looking for, including choosing the step size dx and overloading the y -vector (see case study following) so I decided to write my own routine. One of the most widely used of all the single step Runge-Kutta methods (which include the Euler and Modified Euler methods) is the so called 'Classical' fourth order method with Runge's coefficients. Being a fourth order method (equivalent in accuracy to including the first five terms of the Taylor series expansion of the solution) it requires four evaluations of the derivatives over each increment dx , denoted respectively by dy1 , dy2 , dy3 , dy4 . These are then weighted and summed in a very specific manner to obtain the final derivative dy .

   function 
 
 [x, y, dy] = rk4(deriv,n,x,dx,y) 
 
   %Classical fourth order Runge-Kutta method 
 
 
   %Integrates n first order differential equations 
 
 
   %dy(x,y) over interval x to x+dx 
 
 
   %Izzi Urieli - Jan 21, 2002 
 
 
  x0 = x; 
 
  y0 = y; 
 
  [y,dy1] = feval(deriv,x0,y); 
 
   for 
 
 i = 1:n 
 
  y(i) = y0(i) + 0.5*dx*dy1(i); 
 
   end 
 
 
  xm = x0 + 0.5*dx; 
 
  [y,dy2] = feval(deriv,xm,y); 
 
   for 
 
 i = 1:n 
 
  y(i) = y0(i) + 0.5*dx*dy2(i); 
 
   end 
 
 
  [y,dy3] = feval(deriv,xm,y); 
 
   for 
 
 i = 1:n 
 
  y(i) = y0(i) + dx*dy3(i); 
 
   end 
 
 
  x = x0 + dx; 
 
  [y,dy] = feval(deriv,x,y); 
 
   for 
 
 i = 1:n 
 
  dy(i) = (dy1(i) + 2*(dy2(i) + dy3(i)) + dy(i))/6; 
 
  y(i) = y0(i) + dx*dy(i); 
 
   end 
 
 

The most interesting aspect of this function is the use of the MATLAB function ' feval ' to evaluate the function argument which is passed as a string. This is a fundamental aspect of MATLAB programming, and understanding this is essential to developing MATLAB programs.

The plotting capabilities of MATLAB are extremely simple to use and very versatile. I show a typical output plot for the pendulum case study as follows:

We see that 20 degrees is about the limit of the small angle period (2.006 s), however at 80 degrees the period is seen to be much larger (2.293 s). This could only have been determined by numerically solving the differential equation set.

Case Study 2 - an electric bicycle

Since September 2001 I have been designing an electric bicycle (refer: www.eBiciBikes.com ). In this case study I would like to indicate how I am using the above method to do a dynamic simulation of the eBici bike in order to understand the behaviour of the vehicle under acceleration from a complete stop. Consider the main function as follows:

   %trialrun.m - of eBici 
 
 
   % y(DIS) = distance covered (m) 
 
 
   % y(VEL) = velocity (m/s) 
 
 
   % y(POW) = instantaneous power (watts) 
 
 
   % y(AMP) = motor current (amps) 
 
 
   % Izzi Urieli - Jan 25, 2002 
 
 
  DIS = 1; 
 
  VEL = 2; 
 
  POW = 3; 
 
  AMP = 4; 
 
   global 
 
 VBAT SLOPE 
 
  tfinal = 30; 
  % (s) time span to integrate over 
 
 
  VBAT = 36; 
  % Volts 
 
 
  SLOPE = 0; 
  % Zero slope for this trial 
 
 
  npoints = 100; 
  % number of solution points 
 
 
  dt = tfinal/(npoints - 1); 
 
  t = 0; 
  % Initialise independent variable time 
 
 
  y = [0;0;0;0]; 
  % initial [dist;velo;power;amps] 
 
 
  time(1) = t; 
 
  dist(1) = y(DIS); 
 
  velo(1) = y(VEL)*9/4; 
  % meters/sec to mph 
 
 
  power(1) = y(POW)/20; 
  % watts/20 - for plot scale 
 
 
  amps(1) = y(AMP); 
 
   for 
 
 i = 2:npoints 
 
  [t, y, dy] = rk4( 
 
   'deBici' 
 
 
  ,2,t,dt,y); 
 
  time(i) = t; 
 
  dist(i) = y(DIS); 
 
  velo(i) = y(VEL)*9/4; 
 
  power(i) = y(POW)/20; 
 
  amps(i) = y(AMP); 
 
   end 
 
 
  plot(time,velo, 
  'k' 
 
 ) 
 
  axis([0,30,0,40]) 
 
  grid on 
 
  hold on 
 
  plot(time,power, 
  'r' 
 
 ) 
 
  plot(time,amps, 
  'b' 
 
 ) 
 

Notice that there are only two differential equations, velocity - dy(DIS) , and acceleration - dy(VEL) , however I have loaded the y -vector with two more variables of interest, power - y(POW) , and motor current - y(AMP) . Since these values are evaluated within the derivative function ' deBici ', it is convenient to return them to the main function via the y -vector. The function ' deBici ' follows:

   function 
 
 [y,dy] = deBici(t,y) 
 
   %deBici - eBici differential equations of motion 
 
 
   % y(DIS) = distance covered (m) 
 
 
   % y(VEL) = velocity (m/s) 
 
 
   % y(POW) = instantaneous power (watts) 
 
 
   % y(AMP) = motor current (amps) 
 
 
   % Izzi Urieli - Jan 21, 2002 
 
 
  DIS = 1; 
 
  VEL = 2; 
 
  POW = 3; 
 
  AMP = 4; 
 
   global 
 
 VBAT SLOPE 
 
   % eBici basic parameters: 
 
 
  mt = 100.0; 
 
   % [kg] total mass of bike + rider 
 
 
  mw = 1.0; 
 
   % [kg] mass of each wheel 
 
 
  rw = 0.2115; 
 
   % [m] radius of wheel 
 
 
  jw = mw*(rw^2); 
 
   % [kg m^2] moment of inertia of each wheel (kg m^2) 
 
 
   % equivalent inertial mass of bike + rider: 
 
 
  mi = mt + 2*jw/(rw^2); 
 
   % motor/generator mg62 specs: 
 
 
  winding = 0.395; 
 
   % [ohms] resistance 
 
 
  kt = 0.62; 
 
   % [Nm/A] torque constant 
 
 
  imax = 32; 
 
   % [amps] current limit setting on controller 
 
 
   % wheel/motor gear ratio: 
 
 
  n = 30/22; 
 
   % wheel_torque / motor_torque 
 
 
   % applied torque at rear wheel tw: 
 
 
  tw = (n*kt/winding)*(VBAT - kt*n*y(VEL)/rw); 
 
  twmax = n*kt*imax; % current limit setting on controller 
 
   if 
 
 
  tw > twmax % stall torque by current limit 
 
  tw = twmax; 
 
   end 
 
 
   % current drain 
 
 
  y(AMP) = tw/(n*kt); 
 
   % mechanical power (force*velocity (watts)): 
 
 
  y(POW) = (tw/rw)*y(VEL); 
 
   % resistive forces: drag, slope, rolling resistance 
 
 
  cda = 0.4; 
 
   % [m^2] coef. of drag * frontal area 
 
 
  density = 1.18; 
 
   % [kg/m^3] air density 
 
 
  drag = cda*density*y(2).^2/2; 
 
  cr = 0.005; 
 
   % coefficient of rolling resistance 
 
 
  g = 9.807; 
 
   % [m/s^2] gravity 
 
 
  resist = mt*g*(cr + SLOPE); 
 
  dy(DIS) = y(VEL); 
 
   % velocity 
 
 
  dy(VEL) = (tw/rw - (drag + resist))/mi; 
 
   % acceleration 
 
 
  dy = [dy(DIS); dy(VEL)]; 
 
   % must be a column vector 
 
 

The ability to load the y -vector with the motor current and power values allowed me to develop the following transient performance plot from the main function: