https://1drv.ms/w/s!Ato0Wsph7iDRjr4t-FYt5WfRkFeeiw
We had also partly implemented a C++ code to test the equations, but it was never completely implemented and properly debugged.
https://1drv.ms/f/s!Ato0Wsph7iDRjr08UXUdh4I3l9JLpA
We had this idea of using predictor corrector numerical methods to increase the accuracy of calculations. I have written the following code for the predictor section of the code.
Code: Select all
#include <iostream>
/*
This is a predictor corrector numerical solver:
here we suppose that dx/dt = f, so f is the slope
Knowing the slope of a function at several points,
it would be possible to predicts its value at a next point
The value of slope (f), and the value of x are stored in circular arrays
this array helps us to access these parameters
at this iteration by x[ 0] or f[ 0]
the previous iteration by x[-1] or f[-1]
the iteration before previous one by x[-2] or f[-2]
and so on
knowing these values it would be possible to predict the value of next x
with different accuracies: 2nd, 3rd, or 4th order accuracies are implemented
in this class.
note initialization of predictor is important.
predictor-corrector solvers are multistep methods.
To achieve a 4th order accuracy, only the initial values of x should be
updated by at least a single step 3rd order accurate (or more) numerical method.
note: Adams-Moulton method is more stable than Adams-Bashforth's
*/
template <typename T>
class predictor
{
public:
// initialized at first the values of f and x to steady state
predictor(T initial_value)
{
x.initialize(initial_value);
f.initialize(0);
}
~predictor() {}
void update(T current_value, float deltaT)
{
f.update((current_value - x[0])/deltaT);
x.update(current_value);
}
T predict2ndOrder(float deltaT)
{
return x[0] + deltaT/2*(3*f[0] - f[-1]);
}
T predict3rdOrder(float deltaT)
{
return x[0] + deltaT/12*(23*f[0] - 16*f[-1] + 5*f[-2]);
}
// 4th order predictor Adams-Bashforth
T predict4thOrder(float deltaT)
{
return x[0] + deltaT/24*(55*f[0] - 59*f[-1] + 37*f[-2] - 9*f[-3]);
}
// 4th order predictor-corrector Adams-Moulton
T pc4thOrder(float deltaT)
{
return x[0] + deltaT/24*(9*(predict4thOrder(deltaT)-x[0])/deltaT + 19*f[0] - 5*f[-1] + f[-2]);
}
inline T operator[](int index)
{
return x[index];
}
private:
// size should be in power of 2 = (2,4,8,16,..,128)
// since we are going to use high performance bit-wise modulus
#define size 4
#define SizeNegOne 3
template <typename S>
class carray
{
public:
carray()
{
head = SizeNegOne;
}
// copy constructor
carray(const carray& that)
{
head = that.head;
for (int i=0; i<size; i++)
array[i] = that.array[i];
}
// copy assignment operator
carray& operator=(const carray& that)
{
head = that.head;
for (int i=0; i<size; i++)
array[i] = that.array[i];
return *this;
}
// fill all the elements of the array with
// a specific value
void initialize(S n)
{
for (head=0; head<size; head++)
array[head]=n;
head = SizeNegOne;
}
// add a new value to head of this array
// and delete one element from its tail
// we use array for this purpose so,
// In the other words,
// if you have reached the end of the array,
// start again from its head and therefore
void update(S value)
{
// head index = (head+1)%size
head = (head+1) & SizeNegOne;
array[head]=value;
}
// access elements of this array with negative indexes
// which shoes how many iterations before they have updated
inline S operator[](int index)
{
return array[(head + index) & SizeNegOne];
}
private:
S array[size];
signed char head;
#undef size
#undef SizeNegOne
};
carray<T> f, x; // f=dx/dt
};
int main()
{
// create our objects and initialize them with a fixed value
predictor<float> phi(-64), ph(0);
// suppose that each deltaT = 1 ms, we update the phi value
phi.update(10, 1);
phi.update(11, 1);
phi.update(12, 1);
phi.update(13, 1);
phi.update(14, 1);
// after initialization, we ask the predictor
// its predictions for the next iteration that would take
// deltaT = 1 ms
std::cout<<phi.predict2ndOrder(1)<<std::endl;
std::cout<<phi.predict3rdOrder(1)<<std::endl;
std::cout<<phi.predict4thOrder(1)<<std::endl;
std::cout<<phi.pc4thOrder(1) <<std::endl;
// here I show that it is possible to use this object like
// a regular carray to access previous values of phi
// we need this feature for the calculation of phi and Concentration
std::cout<<phi[0] <<std::endl;
std::cout<<phi[-1]<<std::endl;
std::cout<<phi[-2]<<std::endl;
std::cout<<phi[-3]<<std::endl;
// to make sure that the rule of three is satisfied
ph=phi;
// here we test the copied object
std::cout<<ph.predict2ndOrder(1)<<std::endl;
std::cout<<ph.predict3rdOrder(1)<<std::endl;
std::cout<<ph.predict4thOrder(1)<<std::endl;
std::cout<<ph[0] <<std::endl;
std::cout<<ph[-1]<<std::endl;
std::cout<<ph[-2]<<std::endl;
std::cout<<ph[-3]<<std::endl;
return 0;
}