Romberg Integration (graphing fresnel integrals)

[byrnnryb]

Hello all, I am trying to find a way to graph the Fresnel integrals for a Numerical Analysis class. We have to use Romberg Integration to do this, and I suppose I am confused as to how I can graph this (I am confused when it comes to integrals.)

I believe I have the romberg integration coded correctly below, not sure if I have entered the fresnel integrals (one of them is the function f at the beginning of the code.)

I am using Allegro to try to graph this integral, but suppose I am confused as to how I am supposed to feed the graphing procedure the y values from Romberg Integration.

#include <cstdlib>
#include <iostream>
#include <allegro.h>
#include <cmath>

using namespace std;

// the function
 long double f(long double x)
	{
      return cos((3.14159265/2)*(x*x)) ; 
  
	}


long double Romberg(long double x)
{
        
     long double a,b;
     a=0.0;
     b=5.0;
     int n=4;
     int min = 1;
    int i, j, k;
    long double h, sum, eps;
    eps = 0.001;
 long double R[10][10];
 //STEP 1
    h=b-a;
    R[0][0]=h*(f(a)+f(b))/2.0;
    
    cout<<R[0][0]<<endl;
    
    //STEP 2
                           //cout<<f(a-8);
    for (i=1; i<n; i++)
    {
        sum=0.0;
        
        //STEP 4
        
        for (k=0; k<int(pow(2.0,i-2)-1); k++ )
            {
                  sum +=f(a+(k-0.5)*h);
            }
            
            //STEP 5
        R[1][0]=0.5*R[0][0]+sum*h;
        
        //STEP 6
        for(j=1; j<i; j++)
           {
            R[1][j]=(pow(4.0,(j-1))*R[1][j-1]-(R[0][j-1])/(pow(4.0,j-1)-1));
            }
        for ( int q = 0;q<i;q++)
            {
            cout<<R[1][q]<<"  ";
            } 
        //STEP 7    
            if ((abs(R[1][i]) - R[0][i])<eps and i>=min)
        
            return R[1][i];
            else{
                h=h/2.0;
                for(int j=0; j<i; j++)
                {
                        R[0][j]=R[1][j];
                        }
                       
                        
                        }              
    }  
     return R[1][n];  
}

long double R(long double x)
{
     return Romberg(x);
     }





void graphit(double A, double B, int xsize, int ysize, long double f(long double))
   {
      allegro_init();
      install_keyboard();
    
      int depth = desktop_color_depth();
      if( depth != 0)
         set_color_depth( depth);   
   
      int ret = set_gfx_mode(GFX_AUTODETECT_WINDOWED, xsize,ysize, 0, 0);
      if (ret!=0)
      {
         allegro_message(allegro_error);
         return ;
      }
   
      long double x,  y[SCREEN_W+1], ymin, ymax, scale;
      char title[100];
      y[0] =f( A );
      ymin = y[0];
      ymax = y[0];
      long double dx = (B-A) / SCREEN_W;
   
      for(int i = 1; i <=SCREEN_W ; i++)
      {
         x = A + i * dx;
         y[i] =f( x );
         if( y[i] < ymin) ymin = y[i];
         if( y[i] > ymax) ymax = y[i];
      }
      sprintf( title, "A = %G, B = %G, Minimum = %G, Maximum = %G", A, B, ymin, ymax);
		set_window_title(title);
   
      int white = makecol(255,255,255);
      scale = SCREEN_H / (ymax - ymin);
      for(int i = 0 ; i <= SCREEN_W ; i++)
      {
         int why = SCREEN_H - (int)((y[i] - ymin)*scale);
         putpixel(screen,i,why,white );
      }
      line( screen, 0, (int) (SCREEN_H - (0-ymin)*scale), SCREEN_W,  (int) (SCREEN_H - (0-ymin)*scale), white); //x-axis
      line( screen, (int)((0-A)/(B-A)*SCREEN_W), 0,  (int)((0-A)/(B-A)*SCREEN_W),SCREEN_H, white); //y-axis
   
      clear_keybuf();
      while( !keypressed())
      {
      // sit and spin
      }
   	
      allegro_exit();
   
   } //end graphit




int main(void)
{
   graphit(0.0, 5.0, 640, 480, R);
   Romberg(9);
       system("PAUSE");
    return 0;
}
END_OF_MAIN()  //needed by Allegro