1 | /* $Id: romberg_one.cpp 1370 2009-05-31 05:31:50Z bradbell $ */ |
---|
2 | /* -------------------------------------------------------------------------- |
---|
3 | CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-06 Bradley M. Bell |
---|
4 | |
---|
5 | CppAD is distributed under multiple licenses. This distribution is under |
---|
6 | the terms of the |
---|
7 | Common Public License Version 1.0. |
---|
8 | |
---|
9 | A copy of this license is included in the COPYING file of this distribution. |
---|
10 | Please 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 PROGRAM%// END PROGRAM%1%$$ |
---|
27 | $$ |
---|
28 | |
---|
29 | $end |
---|
30 | */ |
---|
31 | // BEGIN PROGRAM |
---|
32 | |
---|
33 | # include <cppad/cppad.hpp> |
---|
34 | |
---|
35 | namespace { |
---|
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 | |
---|
101 | bool 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 PROGRAM |
---|