1 

/*

2 

* $Id: dfqromb.cpp 1029 20130531 18:27:30Z jsibert $

3 

*

4 

* Author: David Fournier

5 

* Copyright (c) 20092012 ADMB Foundation

6 

*/

7 

/**

8 

\file dfqromb.cpp

9 

Contains routines for numerical integration

10 

*/

11 


12 

#include <admodel.h>

13 

#include "trace.h"

14 


15 

//#define EPS 1.0e4

16 

#define JMAX 50

17 

//#define JMAXP JMAX+1

18 

//#define K 5

19 


20 

class model_parameters;

21 


22 

dvariable trapzd(dvariable (model_parameters::*func)(const dvariable&),double a,double b,int n);

23 

dvariable trapzd(dvariable (model_parameters::*func)(const dvariable&),double a, const dvariable& b,int n);

24 

dvariable trapzd(dvariable (model_parameters::*func)(const dvariable&), const dvariable& a, const dvariable& b, int n);

25 

dvariable trapzd(dvariable (model_parameters::*func)(const dvariable&), const dvariable& a, double b, int n);

26 


27 

void polint(const dvector& xa, const dvar_vector& ya,int n,double x,

28 

const dvariable& y, const dvariable& dy);

29 


30 

/** Romberg integration.

31 

\param func Pointer to a member function of class model_parameters.

32 

\param a Lower limit of integration.

33 

\param b Upper limit of inegration.

34 

\param ns

35 

\return The integral of the function from a to b using Romberg's method

36 

*/

37 

dvariable function_minimizer::adromb(dvariable (model_parameters::*func)(const dvariable&),double a,double b,int ns)

38 

{

39 

const double base = 4;

40 

int MAXN = min(JMAX, ns);

41 

dvar_vector s(1,MAXN+1);

42 


43 

for(int j=1; j<=MAXN+1; j++)

44 

{

45 

s[j] = trapzd(func,a,b,j);

46 

}

47 


48 

for(int iter=1; iter<=MAXN+1; iter++)

49 

{

50 

for(int j=1; j<=MAXN+1iter; j++)

51 

{

52 

s[j] = (pow(base,iter)*s[j+1]s[j])/(pow(base,iter)1);

53 

}
