Add support for the trapezoidal method (RK2)

Error for der(x)=1-x shrinks by halving the step size:
* Euler: x2
* Trapezoidal: x4
* RK4: x16

SimulationRuntime/c/simulation/solver/solver_main.c

@@ -67,14 +67,12 @@
 
 double** work_states;
 
-const int rungekutta_s = 4;
-const double rungekutta_b[4] = { 1.0 / 6.0, 1.0 / 3.0, 1.0 / 3.0, 1.0 / 6.0 };
-const double rungekutta_c[4] = { 0.0, 0.5, 0.5, 1.0 };
-
 typedef struct RK4_DATA
 {
   double** work_states;
   int work_states_ndims;
+  double *b;
+  double *c;
 }RK4_DATA;
 
 
@@ -101,6 +99,7 @@ int solver_main_step(DATA* data, threadData_t *threadData, SOLVER_INFO* solverIn
     retVal = euler_ex_step(data, solverInfo);
     TRACE_POP
     return retVal;
+  case S_TRAPEZOIDAL:
   case S_RUNGEKUTTA:
     retVal = rungekutta_step(data, threadData, solverInfo);
     TRACE_POP
@@ -200,14 +199,36 @@ int initializeSolverData(DATA* data, threadData_t *threadData, SOLVER_INFO* solv
     break;
   }
   case S_RUNGEKUTTA:
+  case S_TRAPEZOIDAL:
   {
     /* Allocate RK work arrays */
 
+    static const int rungekutta_s = 4;
+    static const double rungekutta_b[4] = { 1.0 / 6.0, 1.0 / 3.0, 1.0 / 3.0, 1.0 / 6.0 };
+    static const double rungekutta_c[4] = { 0.0, 0.5, 0.5, 1.0 };
+
+    static const int trapezoidal_s = 2;
+    static const double trapezoidal_b[2] = { 1.0 / 2.0, 1.0 / 2.0 };
+    static const double trapezoidal_c[2] = { 0.0, 1.0 };
+
     RK4_DATA* rungeData = (RK4_DATA*) malloc(sizeof(RK4_DATA));
-    rungeData->work_states_ndims = rungekutta_s;
+
+    if (solverInfo->solverMethod==S_RUNGEKUTTA) {
+      rungeData->work_states_ndims = rungekutta_s;
+      rungeData->b = rungekutta_b;
+      rungeData->c = rungekutta_c;
+    } else if (solverInfo->solverMethod==S_TRAPEZOIDAL) {
+      rungeData->work_states_ndims = trapezoidal_s;
+      rungeData->b = trapezoidal_b;
+      rungeData->c = trapezoidal_c;
+    } else {
+      throwStreamPrint(threadData, "Unknown RK solver");
+    }
+
     rungeData->work_states = (double**) malloc((rungeData->work_states_ndims + 1) * sizeof(double*));
-    for(i = 0; i < rungeData->work_states_ndims + 1; i++)
+    for (i = 0; i < rungeData->work_states_ndims + 1; i++) {
       rungeData->work_states[i] = (double*) calloc(data->modelData->nStates, sizeof(double));
+    }
     solverInfo->solverData = rungeData;
     break;
   }
@@ -315,7 +336,7 @@ int freeSolverData(DATA* data, SOLVER_INFO* solverInfo)
   {
     freeSymEulerImp(solverInfo);
   }
-  else if(solverInfo->solverMethod == S_RUNGEKUTTA)
+  else if(solverInfo->solverMethod == S_RUNGEKUTTA || solverInfo->solverMethod == S_TRAPEZOIDAL)
   {
     /* free RK work arrays */
     for(i = 0; i < ((RK4_DATA*)(solverInfo->solverData))->work_states_ndims + 1; i++)
@@ -740,7 +761,8 @@ static int sym_euler_im_step(DATA* data, threadData_t *threadData, SOLVER_INFO*
 /***************************************    RK4      ***********************************/
 static int rungekutta_step(DATA* data, threadData_t *threadData, SOLVER_INFO* solverInfo)
 {
-  double** k = ((RK4_DATA*)(solverInfo->solverData))->work_states;
+  RK4_DATA *rk = ((RK4_DATA*)(solverInfo->solverData));
+  double** k = rk->work_states;
   double sum;
   int i,j;
   SIMULATION_DATA *sData = (SIMULATION_DATA*)data->localData[0];
@@ -757,13 +779,13 @@ static int rungekutta_step(DATA* data, threadData_t *threadData, SOLVER_INFO* so
     k[0][i] = stateDerOld[i];
   }
 
-  for(j = 1; j < rungekutta_s; j++)
+  for (j = 1; j < rk->work_states_ndims; j++)
   {
     for(i = 0; i < data->modelData->nStates; i++)
     {
-      sData->realVars[i] = sDataOld->realVars[i] + solverInfo->currentStepSize * rungekutta_c[j] * k[j - 1][i];
+      sData->realVars[i] = sDataOld->realVars[i] + solverInfo->currentStepSize * rk->c[j] * k[j - 1][i];
     }
-    sData->timeValue = sDataOld->timeValue + rungekutta_c[j] * solverInfo->currentStepSize;
+    sData->timeValue = sDataOld->timeValue + rk->c[j] * solverInfo->currentStepSize;
     /* read input vars */
     externalInputUpdate(data);
     data->callback->input_function(data, threadData);
@@ -778,9 +800,9 @@ static int rungekutta_step(DATA* data, threadData_t *threadData, SOLVER_INFO* so
   for(i = 0; i < data->modelData->nStates; i++)
   {
     sum = 0;
-    for(j = 0; j < rungekutta_s; j++)
+    for(j = 0; j < rk->work_states_ndims; j++)
     {
-      sum = sum + rungekutta_b[j] * k[j][i];
+      sum = sum + rk->b[j] * k[j][i];
     }
     sData->realVars[i] = sDataOld->realVars[i] + solverInfo->currentStepSize * sum;
   }

SimulationRuntime/c/util/simulation_options.c

@@ -360,13 +360,14 @@ const char *SOLVER_METHOD_NAME[S_MAX] = {
   "lobatto6",
   "symEuler",
   "symEulerSsc",
+  "trapezoidal",
   "qss"
 };
 
 const char *SOLVER_METHOD_DESC[S_MAX] = {
   "unknown",
   "euler",
-  "rungekutta",
+  "rungekutta (fixed step, order 4)",
   "dassl with colored numerical jacobian, with interval root finding - default",
   "optimization",
   "radau5 [sundial/kinsol needed]",
@@ -377,6 +378,7 @@ const char *SOLVER_METHOD_DESC[S_MAX] = {
   "lobatto6 [sundial/kinsol needed]",
   "symbolic implicit euler, [compiler flag +symEuler needed]",
   "symbolic implicit euler with step-size control, [compiler flag +symEuler needed]",
+  "trapezoidal method (fixed step, order 2)",
   "qss"
 };
 

SimulationRuntime/c/util/simulation_options.h

@@ -127,6 +127,7 @@ enum SOLVER_METHOD
   S_LOBATTO6,      /* 10 */
   S_SYM_EULER,     /* 11 */
   S_SYM_IMP_EULER, /* 12 */
+  S_TRAPEZOIDAL,   /* 13 */
   S_QSS,
 
   S_MAX