Research - Programming

[demonicgerbil]

So, I've been doing research again. Programming.

Right now I'm coding up a Four Box system of interacting populations with Forward Euler. Basically it's a test by Dr. Guidry to see if I can hack it. :p

You can find my program behind the cut. ... How do I get a monospace font in html again? Oh yeah, the pre tag. Anyway, this is mostly done. I just need to fix the "evolve populations" section of code so that it does the math right. This might mean I need to change how I handle the fluxes. It would be easy to code the interactions by hand, but I'm trying to leave it in a more general form so the code has more reusability. Once I get this working and giving the 'right' answers, it'll be time to redo this code to implement 4th order Runge-Kutta, which isn't that much of an extension from what's here now. I could also add a dynamic time-step code to it, but that's getting too fancy for me, I think. :p

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

Filename      : 4box.c
Author        : Elisha D. Feger
Description   : This program implements Forward Euler with flux constraint for a
                simple four box reaction network.

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

#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <iostream.h>
#include <assert.h>
#include <ctype.h>
#include <math.h>
#include <time.h>

int main()
 {
   double curPop[4], newPop[4], fluxes[4][4],k[4][4];
   double numPop;
   int i, m, n;
   double deltaT;
   double maxTime;
   double curTime;


   //Setup Time information
   deltaT=1;
   maxTime=100;
   curTime=0;
   //End setup time information

   //Setup the initial mass fractions
   numPop=1000000;
   curPop[0]=.2*numPop;
   curPop[1]=.8*numPop;
   curPop[2]=.00001*numPop;
   curPop[3]=.0000001*numPop;
   //End setup the initial mass fractions

   //Setup the flux constants kij
   //Set flux constant matrix to 0
   for(m=0;m<4;m++)
     {
       for(n=0;n<4;n++)
	 {
	   k[m][n]=0;
	 }
     }
   //End set flux constant matrix to 0
   //Set non-zero elements of flux constant
   k[0][1]=.002;
   k[0][2]=.0025;
   k[0][3]=.02;
   k[1][0]=.5;
   k[2][0]=.05;
   k[3][2]=.005;
   //End set non-zero elements of the flux constant
   //End setup the flux constants kij

   //Do the Forward Euler method
   while(curTime<maxTime)
     {
   //Start the inner program loop
   //Calculate the derivitives
   for(m=0;m<4;m++)
     {
       for(n=0;n<4;n++)
	 {
	   fluxes[m][n]=k[m][n]*curPop[m];
	 }
     }
   //End calculate the derivitives

   //Apply Flux constraint
   for(m=0;m<4;m++)
     {
       for(n=0;n<4;n++)
         {
           if(fluxes[m][n]<0)
	     {
	       fluxes[m][n]=0;
	     }
         }
     }
   //End apply flux constraint

   //Evolve populations
   for(m=0;m<4;m++)
     {
       newPop[m]=curPop[m];
       for(n=0;n<4;n++)
         {
           newPop[m]=newPop[m]+fluxes[m][n]*deltaT;
         }
     }
   //End evolve populations

   //Output populations
   cout << curTime << " " << curPop[0] << " " << curPop[1] << " " << curPop[2] << " " << curPop[3] << endl;
   //End output populations

   //Setup for next time iteration
   //Setup populations
   for(i=0;i<4;i++)
     {
       curPop[i]=newPop[i];
       newPop[i]=0;
     }
   //End setup populations

   //Increment time
   curTime=curTime+deltaT;
   //End Increment time

   //End the inner program loop
 }
   //End the Forward Euler Method
 }