source: trunk/test_more/romberg_one.cpp @ 2439

Last change on this file since 2439 was 2439, checked in by bradbell, 8 years ago

Change BEGIN PROGRAM to BEGIN C++ using followind sed commands:

s| BEGIN PROGRAM| BEGIN C++|
s| END PROGRAM| END C++|

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