source: trunk/test_more/romberg_one.cpp @ 678

Last change on this file since 678 was 678, checked in by bradbell, 14 years ago

trunk: Change CppAD/*.h to lower case names and .hpp extension.

File size: 2.2 KB
Line 
1/* --------------------------------------------------------------------------
2CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-06 Bradley M. Bell
3
4CppAD is distributed under multiple licenses. This distribution is under
5the terms of the
6                    Common Public License Version 1.0.
7
8A copy of this license is included in the COPYING file of this distribution.
9Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
10-------------------------------------------------------------------------- */
11
12/*
13$begin RombergOne.cpp$$
14$spell
15        Romberg
16$$
17
18$section One Dimensional Romberg Integration: Example and Test$$
19
20$index Romberg, example$$
21$index example, Romberg$$
22$index test, Romberg$$
23
24$code
25$verbatim%example/romberg_one.cpp%0%// BEGIN PROGRAM%// END PROGRAM%1%$$
26$$
27
28$end
29*/
30// BEGIN PROGRAM
31
32# include <CppAD/cppad.hpp>
33
34namespace {
35        class Fun {
36        private:
37                const size_t degree;
38        public:
39                // constructor
40                Fun(size_t degree_) : degree(degree_) 
41                { }
42
43                // function F(x) = x^degree
44                template <class Float>
45                Float operator () (const Float &x)
46                {       size_t i;
47                        Float   f = 1;
48                        for(i = 0; i < degree; i++)
49                                f *= x;
50                        return f;
51                }
52        };
53
54        template <class Float>
55        bool RombergOneCase(void)
56        {       bool ok = true;
57                size_t i;
58
59                size_t degree = 4;
60                Fun F(degree);
61
62                // arguments to RombergOne
63                Float a = 0.;
64                Float b = 1.;
65                Float r;
66                size_t n = 4;
67                Float e;
68                size_t p;
69
70                // int_a^b F(x) dx =
71                //      [ b^(degree+1) - a^(degree+1) ] / (degree+1)
72                Float bpow = 1.;
73                Float apow = 1.;
74                for(i = 0; i <= degree; i++)
75                {       bpow *= b;
76                        apow *= a;
77                } 
78                Float check = (bpow - apow) / Float(degree+1);
79
80                // step size corresponding to r
81                Float step = (b - a) / exp(log(Float(2.))*Float(n-1));
82                // step size corresponding to error estimate
83                step *= Float(2.);
84                // step size raised to a power
85                Float spow = Float(1.);
86
87                for(p = 0; p < n; p++)
88                {       spow = spow * step * step;
89
90                        r = CppAD::RombergOne(F, a, b, n, p, e);
91
92                        ok  &= e < (degree+1) * spow;
93                        ok  &= CppAD::NearEqual(check, r, Float(0.), e);       
94                }
95
96                return ok;
97        }
98}
99
100bool RombergOne(void)
101{       bool ok = true;
102        using CppAD::AD;
103
104        ok     &= RombergOneCase<double>();
105        ok     &= RombergOneCase< AD<double> >();
106        ok     &= RombergOneCase< AD< AD<double> > >();
107        return ok;
108}
109
110// END PROGRAM
Note: See TracBrowser for help on using the repository browser.