Simulations Using Integration Algorithms - First Order Systems

The Euler Integration Algorithm

The Euler Integration Algorithm

        In another lesson we have presented the Euler Integration Algorithm.  That algorithm assumes that you have a value of a variable, x(t), a state equation that allows you to compute dx/dt, and that you want to compute x(t + Dt).  The algorithm is simple.  Here it is:

x(t + Dt) ~= x(t) + Dt * (dx/dt)|t.

The algorithm is applied repetitively to compute a solution for the state at equally spaced intervals of time.  The time intervals are Dt apart.

Implementing The Euler Algorithm In Mathcad

        Here is a picture of a Mathcad workspace that computes the response of a linear, first order system and stores the result in an array called EulerResult.  Copy this code into Mathcad and plot EulerResult.  This code is what was used to generate the plots in the discussion of how the Euler algorithm works.

Implementing The Euler Algorithm In C

               Here is a C program that you can use to do simulations using the Euler integration algorithm.  This problem simulates the same example system as the Mathcad example above.  You can copy the entire program from this page and copy it into your compiler.  This program was tested in Borland C++ Builder 5.

 #include <stdio.h>
#include <conio.h>

        /* Variables defined here are global.  */

      float x[11],xdot[11];             /* Define states and state derivatives */
      float t,dt;                               /* Time and integration interval global*/
      int nstate;                             /* Number of states and state equations*/

      /***************** State Equations  *******************/

         state_eqns(void)
             {
                 float Tau = 1.0, Input = 1.0, Gdc = 1.0;
                 xdot[0] = -x[0]/Tau + (Gdc/Tau)*Input;
                 return;
             }

     /* Euler Integration Algorithm **********************************************/
         void euler(void)
             {
                 int count = 0;
                 state_eqns();
                 while (count<nstate) {
                     x[count] = x[count] + dt*xdot[count];
                     count = count+1;
                  }
                  t = t + dt;
                  return;
             }

     /***********************************************************************/

     int main(void)
         {
             int     index,iterations;
             float  Tmax;

     /* Data Input Section *******************************************/
             printf(" Input data for the simulation. \n");
             printf("\n");
             printf("Input the total time for your simulation.\n");
             scanf("%f",&Tmax);
             printf("Input the sample time for your simulation.\n");
             scanf(" %f",&dt);

     /* Initialization section ***************************************/

             t=0.;          /* Initialize time & states*/
             index=0;
             nstate=1;  /* First order system - 1 state */
             while (index<nstate) {
                 x[index]=0.;
                 index = index+1;
              }
             printf("\n");
             printf("Time (sec)         Output\n");
             printf("\n");
             printf("%12f%12f\n",t,x[0]);

             index = 0;
             iterations = Tmax/dt;

     /* Do simulation calculations and print results ************************/

             while (index<=iterations) {
                 euler();
                 printf("%12.4f%12.4f\n",t,x[0]);
                 index = index + 1;
             }

     while (!kbhit()) /* do nothing */ ;
     return 0;
     }

        This program is long enough to be confusing, so we will dissect it a bit.  First, the state equations are implemented in this code.
 /***************** State Equations  *******************/

         state_eqns(void)
             {
                 float Tau = 1.0, Input = 1.0, Gdc = 1.0;
                 xdot[0] = -x[0]/Tau + (Gdc/Tau)*Input;
                 return;
             }

You can change this code to simulate other first order systems by putting the differential equation into xdot[0].  You can add other state equations for systems of higher order by using xdot[1], etc.  Note that the array, x[ ], is dimensioned 11 at the very beginning of the code, so if you want more than 11 states you will have to change that dimensioning.

        The Euler integration algorithm is implemented in separate code, but it calls the state equation function since it needs that information to update the state(s).  Here is the code.

state_eqns(void)
             {
                 float Tau = 1.0, Input = 1.0, Gdc = 1.0;
                 xdot[0] = -x[0]/Tau + (Gdc/Tau)*Input;
                 return;
             }
In addition, in the code above, there is a main program that prompts the user for necessary information, calls the Euler integration algorithm function which in turn calls the state equation function.  When each computation is finished, results are printed to the screen.

        You can modify this program to store data in a file, and you can also modify it for more complex systems.