Introduction - Richardson Extrapolation

Football Fan? Why not check out www.playerrank.co.uk to vote for your favourite ever player?

I'm a final year PhD student, studying Solar Physics, and I decided to start this blog to hopefully help out some fellow students - or anyone else  - who are tackling some of the same problems I've encountered. I don't claim that any code I post is perfect, in fact I'm sure most of it could be improved and I'd actively encourage any suggestions for improvements, or questions, and I'll do my best to respond.

I code in a variety of languages, it's hard to say which one I use most, but it's between Fortran, IDL and C++. My favourite language to code in is C++ because this was the first language I learnt, I use Fortran because my supervisor wrote some code in Fortran that I inherited, and IDL is great for analysing the outputs of my simulations.

This week I'm going to post about some C++ code I wrote to do Richardson Extrapolation, this technique is useful for the integration of differential equations using a Bulirsh-Stoer method. This works by repeatedly integrating the equation, successively halving the value of h used for integration, and then extrapolating to get the value that would be achieved when h=0.  

Above is a poorly drawn illustration of the process. Richardson Extrapolation is ideal for this kind of thing when the error terms in your estimate are all of even powers. For integration of ODE's the modified midpoint method produces values that have only even powers of h, so combining the two methods provides a quick method of integration.

The general formula for Richardson Extrapolation, taken from wikipedia, is presented below:

So now to the code, below is a function that returns a double precision value, corresponding to the estimated value when h=0. It takes two arguments, an array of maximum length 100, containing values of y that are generated from successively halved h values. There is a section of code commented out to print the array "calc" where the working out is done. I suggest uncommenting this to get a feel for how the code works.


double rich_extrap(double y[], int length){


double calc[100][100]={0};
double result;




for(int i=0; i<length; i++){
calc[0][i]=y[i];
}




for(int i=1; i<length; i++){


for(int j=0; j<(length-i); j++){


calc[i][j] = calc[i-1][j+1] + (calc[i-1][j+1] - calc[i-1][j])/pow(2.0,i); 


}


}


/* THE CODE BELOW PRINTS THE CALC ARRAY
for(int i=0; i < 10; i++){
for(int j=0; j < 10; j++){


cout << calc[i][j] << " ";
}
cout << endl;
}
*/


result=calc[length-1][0];
return result;
}


Again, I'd like to point out that I welcome any comments or queries. I'd also love to hear from you if you found this remotely useful. I'm going to try and post weekly, possibly more regularly if I encounter some interesting problems or things that might be useful. I'm thinking that my next post might be about generating a circle on a random plane...