*** empty log message ***

Makefile

@@ -7,7 +7,7 @@ PROJNAME = miscmaths
 USRINCFLAGS = -I${INC_NEWMAT} -I${INC_PROB} -I${INC_ZLIB}
 USRLDFLAGS = -L${LIB_NEWMAT} -L${LIB_PROB} -L${LIB_ZLIB}
 
-OBJS = miscmaths.o optimise.o miscprob.o kernel.o histogram.o base2z.o t2z.o f2z.o volume.o volumeseries.o minimize.o cspline.o sparse_matrix.o sparsefn.o
+OBJS = miscmaths.o optimise.o miscprob.o kernel.o histogram.o base2z.o t2z.o f2z.o volume.o volumeseries.o minimize.o cspline.o sparse_matrix.o sparsefn.o rungekutta.o
 #OBJS = miscmaths.o optimise.o miscprob.o kernel.o histogram.o base2z.o t2z.o f2z.o volume.o volumeseries.o minimize.o cspline.o
 
 LIBS = -lutils -lnewmat -lprob -lm

rungekutta.cc

+/*  rungekutta.cc
+
+    Mark Woolrich, FMRIB Image Analysis Group
+
+    Copyright (C) 2002 University of Oxford  */
+
+/*  CCOPYRIGHT  */
+
+#include "rungekutta.h"
+
+namespace MISCMATHS {
+
+
+void rk(ColumnVector& ret, const ColumnVector& y, const ColumnVector& dy, float x, float h, const Derivative& deriv,const ColumnVector& paramvalues)
+{ 
+  Tracer tr("rk"); 
+
+  float hh=h*0.5;
+  float xh=x+hh;
+  
+  //first step
+  ColumnVector yt=y+hh*dy;
+  
+  //second step
+  ColumnVector dyt = deriv.evaluate(xh,yt,paramvalues);
+  yt=y+hh*dyt;
+  
+  //third step
+  ColumnVector dym = deriv.evaluate(xh,yt,paramvalues);
+  yt=y+h*dym;
+  dym=dym+dyt;
+  
+  //fourth step
+  dyt = deriv.evaluate(x+h,yt,paramvalues);
+  
+  //addup
+  ret = y+h*(dy+dyt+2*dym)/6;
+}
+
+void rkqc(ColumnVector& y, float& x, float& hnext, ColumnVector& dy, float htry, float eps, const Derivative& deriv,const ColumnVector& paramvalues)
+{
+  Tracer tr("rkqc"); 
+
+  float xsav = x;
+  ColumnVector dysav = dy;
+  ColumnVector ysav = y;
+  float h = htry;
+  float hdid;
+  ColumnVector ytemp;
+
+  while(true)
+    {
+      // take 2 1/2 step sizes
+  
+      // first 1/2 step
+      float hh=h*0.5;
+
+      rk(ytemp,ysav,dysav,xsav,hh,deriv,paramvalues);
+  
+      // second 1/2 step
+      x=xsav+hh;
+      dy = deriv.evaluate(x,ytemp,paramvalues);
+      rk(y,ytemp,dysav,xsav,hh,deriv,paramvalues);
+      x=xsav+h;
+      if(x==xsav) cerr << "step size too small" << endl;
+
+      // take large step size
+      rk(ytemp,ysav,dysav,xsav,h,deriv,paramvalues);
+   
+      // eval accuracy
+      float errmax = 0.0;
+      for(int i=1; i<=y.Nrows(); i++)
+	{
+	  //errmax=max(abs((y-ytemp)./y));
+	  
+	  float tmp = fabs((y(i)-ytemp(i))/y(i));
+	  if(tmp > errmax) errmax = tmp;
+	}
+
+      errmax=errmax/eps;
+      
+      if(errmax <=1.0) 
+	{
+	  // step OK, compute step size for next step
+	  hdid=h;
+	  
+	  if(errmax>6e-4)
+	    hnext=h*std::exp(-0.2*std::log(errmax));
+	  else
+	    hnext=4*h;
+	  
+	  break;
+      }
+      else 
+	{
+	  // step too large,
+	  h=h*std::exp(-0.25*std::log(errmax));
+	}
+    }
+
+  y = y+(y-ytemp)/15;
+}
+
+void runge_kutta(Matrix& yp, ColumnVector& xp, ColumnVector& hp, const ColumnVector& ystart, float x1, float x2, float eps, float hmin, const Derivative& deriv,const ColumnVector& paramvalues)
+{
+  Tracer tr("runge_kutta"); 
+
+  int MAXSTEP=1000;
+
+  ColumnVector y = ystart;
+
+  float x=x1;
+  xp.ReSize(MAXSTEP,1);
+  xp = 0;
+  xp(1) =x1;
+
+  float h=hp(1);
+  hp.ReSize(MAXSTEP,1);
+  hp = 0;
+
+  yp.ReSize(MAXSTEP,y.Nrows());
+  yp = 0;
+
+  int kout=1;
+  
+  ColumnVector dy;
+
+  for(int k=1; k <= MAXSTEP; k++)
+    { 
+      dy = deriv.evaluate(x,y,paramvalues);
+
+      // store results:
+      xp(kout)=x;
+      yp.Row(kout)=y;
+      hp(kout)=h;
+  
+      kout=kout+1;
+    
+      // stop overshoot of step past x2:
+      if((x+h-x2)*(x+h-x1)>0) h=x2-x;
+
+      float hnext = 0.0;
+      rkqc(y,x,hnext,dy,h,eps,deriv,paramvalues);
+
+      if((x-x2)*(x2-x1) >= 0.0)
+      {
+	xp(kout)=x;
+	yp.Row(kout)=y;
+	hp(kout)=h;
+	//kout=kout+1;
+	
+        xp = xp.Rows(1,kout);
+	yp = yp.Rows(1,kout);
+
+	return;
+      }
+      else
+      {
+	if(hnext<=hmin) cerr << "step size too small" << endl;
+	h=hnext;
+      } 
+      
+    }
+  cerr << "too many steps" << endl;
+}
+
+}

rungekutta.h

+/*  rungekutta.h
+
+    Mark Woolrich - FMRIB Image Analysis Group
+
+    Copyright (C) 2002 University of Oxford  */
+
+/*  COPYRIGHT  */
+
+#if !defined(rungekutta_h)
+#define rungekutta_h
+
+#include <iostream>
+#include <fstream>
+#include <string>
+#include <vector>
+
+#include "newimage/newimageall.h"
+
+using namespace NEWIMAGE;
+
+namespace MISCMATHS {
+
+class Derivative
+{
+public:
+  Derivative(int pny) : ny(pny), dy(pny) {}
+
+  // x is time point to evaluate at
+  // y is state variables
+  // paramvalues are "constants" in the diff eqn
+  virtual const ColumnVector& evaluate(float x,const ColumnVector& y,const ColumnVector& paramvalues) const = 0;
+
+  virtual ~Derivative(){};
+
+protected:
+  int ny;
+  mutable ColumnVector dy;
+};
+
+void rk(ColumnVector& ret, const ColumnVector& y, const ColumnVector& dy, float x, float h, const Derivative& deriv,const ColumnVector& paramvalues);
+
+void rkqc(ColumnVector& y, float& x, float& hnext, ColumnVector& dy, float htry, float eps, const Derivative& deriv,const ColumnVector& paramvalues);
+
+void runge_kutta(Matrix& yp, ColumnVector& xp, ColumnVector& hp, const ColumnVector& ystart, float x1, float x2, float eps, float hmin, const Derivative& deriv,const ColumnVector& paramvalues);
+
+}
+
+#endif