Lindo Systems

! Estimate the parameters of a differential equation
  to describe the velocity in ft/sec and distance in ft
  of an object falling through the atmosphere,
  using Runge-Kutta approximation and minimization
  of squared errors;
! Keywords: Differential equations, Dynamic system, LINGO, ODE, Runge Kutta;
Sets:
 TIME: t, v, y, yobs, h,
     k1v, k2v, k3v, k4v,
     k1y, k2y, k3y, k4y;
ENDSETS
DATA:
  TIME = 1..21;
 ! Time and observed position;
  t, yobs = 
   0    0 
   1   16
   2   62 
   3  134 
   4  222 
   5  320 
   6  421 
   7  525 
   8  630 
   9  734 
  10  839
  11  944 
  12 1049 
  13 1154 
  14 1259 
  15 1364 
  16 1470
  17 1575
  18 1680
  19 1785 
  20 1890;
 ENDDATA

  MIN = sumerrs;! Min sum of errors of observed vs. predicted;
   sumerrs = @SUM( TIME(i): ( y(i)- yobs(i))^2);
  ! This is a completely determined, triangular system
   if p0, p1, p2, p3, and v(1) are specified;
!   @free(p0);
!   @free( p1);
!   @free( p2);
!   @free( p3);
!   Turning on the bounds helps bring the solution near the optimum;
   @bnd(30, p0, 33);
   @bnd(-.04, p1, -.00005);
   @bnd(-.005, p2, -.000001);
   @bnd(-.000005, p3, .0000001);

  @for( TIME(i) | i #GT# 1:
    h(i) = t(i)- t(i-1);
 ! Our model of velocity change is:
     dv/dt = p0 + p1* v + p2* v^2 + p3* v^3;
    ! Four (Runge-Kutta) estimates of dv/dt;
    @free( k1v(i)); @free( k2v(i));@free( k3v(i));@free( k4v(i));
    k1v(i) = p0 + p1*   v(i-1)                   + p2*  v(i-1)^2                  + p3*v(i-1)^3;
    k2v(i) = p0 + p1* ( v(i-1)+0.5* h(i)*k1v(i)) + p2*( v(i-1)+0.5* h(i)* k1v(i))^2 + p3*( v(i-1)+0.5* h(i)* k1v(i))^3;
    k3v(i) = p0 + p1* ( v(i-1)+0.5* h(i)*k2v(i)) + p2*( v(i-1)+0.5* h(i)* k2v(i))^2 + p3*( v(i-1)+0.5* h(i)* k2v(i))^3;
    k4v(i) = p0 + p1* ( v(i-1)+     h(i)*k3v(i)) + p2*( v(i-1)+     h(i)* k3v(i))^2 + p3*( v(i-1)    + h(i)* k3v(i))^3;
 ! Our model of distance change is:
     dy/dt = v;
   !  Four(Runge-Kutta) estimates of dy/dt;
    @free( k1v(i));@free( k2v(i));@free( k3v(i));@free( k4v(i));
    k1y(i) = v(i-1);
    k2y(i) = v(i-1) + 0.5* h(i)* k1v(i);
    k3y(i) = v(i-1) + 0.5* h(i)* k2v(i);
    k4y(i) = v(i-1) +      h(i)* k3v(i);

   ! Finally, take weighted average to move to next step;
    v(i) = v(i-1) + ( h(i)/6)*( k1v(i) +2* k2v(i) +2* k3v(i) + k4v(i));
    y(i) = y(i-1) + ( h(i)/6)*( k1y(i) +2* k2y(i) +2* k3y(i) + k4y(i));
      );