I have a predictor method that takes differential equations of 214 biological species and integrates them in time. Right now it is able to successfully integrate the differential equations when the time step is positive, so right now it can integrate forwards in time and I was looking to see if it could be integrated backwards in time. By that I mean I was looking to see if the differential equations could be integrated using a negative time step, based on using the following method for integrating the differential equations.
The predictor method uses the forward Euler equations for each species in order to predict its value at the next time step.
The corrector method uses the trapezoidal rule. So the predictor-corrector method is based on Heun's method. It also uses three terms in the gear integration method which is implemented alongside Heun's method.
So I was wondering, is it possible to use this integration with a negative time step and integrate the equations using a negative time step?
If so, how would I go about doing that?
The code for the predictor method, predic.c, and the corrector method, correc.c, are shown below that is needed too.
predic.c:
#include <stdio.h> #include <math.h> void predic (met,dmet,ddmet,d3met,d4met,d5met,NMET,dt) double *met , *dmet, *ddmet, *d3met,*d4met,*d5met; double dt ; int NMET ; { double c1 , c2 ,c3; int i ; c1 = dt ; c2 = c1*dt/2.0 ; c3 = c2*dt/3.0 ; for(i=1;i<=NMET;i++) { /*|| dmet >= 200000000.000 || dmet <= -100000000.0*/ if(met<=1.0E-20 || dmet >=1.0E300 || dmet <= -1.0E300 ) { met=0.0; dmet=0.0; ddmet=0.0; d3met=0.0; } met = met+c1*(dmet) + c2*(ddmet) + c3*(d3met); dmet = dmet + c1*(ddmet) + c2*(d3met); ddmet = ddmet + c1*(d3met); /*+ c4*(d4met) + c5*(d5met)+ c3*(d4met) + c4*(d5met)+ c2*(d4met)+ c3*(d5met)*/ } }
correc.c:
#include <stdio.h> #include <math.h> void correc(met,dmet,ddmet,d3met,dmetn,NMET,dt) double *met , *dmet, *dmetn ; double *d3met, *ddmet ; double dt; int NMET ; { double gear0 , gear2 ,gear3; double cor; double c1, c2, c3, cmet, cddmet, cd3met; int i ; gear0 = 3.0/8.0; gear2 = 3.0/4.0; gear3 = 1.0/6.0; c1 = dt ; c2 = c1*dt/2.0 ; c3 = c2*dt/3.0 ; cmet = gear0*c1; cddmet = gear2*c1/c2; cd3met = gear3*c1/c3; for(i=1;i<=NMET;i++) { cor = dmetn - dmet; met = met + cmet*cor; dmet = dmetn; ddmet = ddmet + cddmet*cor; d3met = d3met + cd3met*cor; } }
