sandbox/Antoonvh/root-finding.c

Root Finding Using Adaptive Grids

We can find the zeros of a smooth function using adaptive grids. This can be achieved by iteratively refining arround the neighborhood where the function changes sign. Note that this algorithm is not fool proof, e.g. it can not find x=0 for f(x)=x2. But the concept can be easily extended to work for such roots as well. Furtheremore, it migth produce floating point errors if there is a root that is a rational number.

In this example we try to find the non trivial zeros of the first Bessel function of the first kind for x is between zero and 20. I.e.

J1(x)=0, for x(0,20).

We assume that the math library has an accurate enough definitnion of this function.

#include "grid/bitree.h"

#define function(x)  (j1(x))

scalar f[],rf[]; 
int maxlevel = 19; 

int main(){
  FILE * fp1 = fopen("j1.dat","w");
  FILE * fp2 = fopen("roots.dat","w");
  L0=20;
  X0=0.0001;
  N=(128);
  init_grid(N);
  foreach(){
    f[]=function(x);
    rf[]=1/f[];
  }
  int tot = 1;
  while (tot>0) 
    { 
      astats g = adapt_wavelet({rf},(double[]){1.},maxlevel,8); 
      foreach(){
	f[]=function(x);
	rf[]=1/f[];
      }
      tot = g.nf+g.nc;
    }
  int n = 0,nn=0;
  foreach(){
    n++;
    fprintf(fp1,"%g\t%g\n",x,f[]);
    if (level == maxlevel&&(f[]*f[1]<0)){
      fprintf(ferr,"x=%.*g\n",10,(x*f[1] + (x+(L0/(pow(2,maxlevel))))*f[0])/(f[]+f[1]));
      fprintf(fp2,"%g\t%g\n",(x*f[1] + (x+(L0/(pow(2,maxlevel))))*f[0])/(f[]+f[1]) , 0.);
    }
  }
  foreach_cell()
    nn++;
  fprintf(ferr,"with error margin = %.*f\nUsed grid cells: %d or %d if you count non-leaf cells as well.\n",10,L0/pow(2,maxlevel),n,nn);
}

Results

This Produces the following output in the terminal:

x=3.831734287
x=7.015553147
x=10.17344108
x=13.32339081
x=16.47077911
x=19.6158736
with error margin = 0.0000381470
Used grid cells: 2628 or 5565 if you count non-leaf cells as well.

This seems to correspond with online references down to the mentioned accuracy.

We also plot the results: