Sacado Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
dfad_taylor_example.cpp
Go to the documentation of this file.
1 // $Id$
2 // $Source$
3 // @HEADER
4 // ***********************************************************************
5 //
6 // Sacado Package
7 // Copyright (2006) Sandia Corporation
8 //
9 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
10 // the U.S. Government retains certain rights in this software.
11 //
12 // This library is free software; you can redistribute it and/or modify
13 // it under the terms of the GNU Lesser General Public License as
14 // published by the Free Software Foundation; either version 2.1 of the
15 // License, or (at your option) any later version.
16 //
17 // This library is distributed in the hope that it will be useful, but
18 // WITHOUT ANY WARRANTY; without even the implied warranty of
19 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 // Lesser General Public License for more details.
21 //
22 // You should have received a copy of the GNU Lesser General Public
23 // License along with this library; if not, write to the Free Software
24 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
25 // USA
26 // Questions? Contact David M. Gay (dmgay@sandia.gov) or Eric T. Phipps
27 // (etphipp@sandia.gov).
28 //
29 // ***********************************************************************
30 // @HEADER
31 
32 // taylor_example
33 //
34 // usage:
35 // taylor_example
36 //
37 // output:
38 // prints the results of computing a single Taylor series expansion of
39 // the solution to:
40 //
41 // dx/dt = 1 + x^2, x(0) = x0 = 1.0;
42 //
43 // The exact solution is x(t) = tan(t + atan(x0)) = tan(t + pi/4)
44 // Also computes the derivative of the Taylor series solution with
45 // respect to x0. The exact series derivative is
46 // dx/dx0(t) = 1/2 * sec^2(t + atan(x0))
47 
48 #include <iostream>
49 
50 #include "Sacado_No_Kokkos.hpp"
51 
52 // Function implementing RHS of ODE
53 template <typename ScalarT>
54 void func(ScalarT& f, const ScalarT& x) {
55  f = 1.0 + x*x;
56 }
57 
58 int main(int argc, char **argv)
59 {
60  double x0 = 1.0; // Initial condition
61  int deg = 40; // Degree of Taylor series solution
62 
63  Sacado::Tay::Taylor<double> x = x0; // Taylor polynomial for independent
64  Sacado::Tay::Taylor<double> f; // Taylor polynomial for dependent
65 
66  // Reserve space for degree deg coefficients
67  x.reserve(deg);
68 
69  // Compute Taylor series solution to dx/dt = f(x)
70  for (int k=0; k<deg; k++) {
71  func(f, x);
72 
73  // Set next coefficient
74  x.resize(k+1, true);
75 
76  // x_{k+1} = f_k / (k+1)
77  x.fastAccessCoeff(k+1) = f.coeff(k) / (k+1);
78  }
79 
80  // Compute derivative w.r.t. x0
83  func(f_fad, x_fad);
84  Sacado::Tay::Taylor<double> dxdx0(deg);
85  dxdx0.fastAccessCoeff(0) = 1.0;
86  for (int k=0; k<deg; k++) {
87  dxdx0.fastAccessCoeff(k+1) = 0.0;
88  for (int j=0; j<=k; j++)
89  dxdx0.fastAccessCoeff(k+1) += f_fad.dx(0).coeff(k-j) * dxdx0.coeff(j);
90  dxdx0.fastAccessCoeff(k+1) /= k+1;
91  }
92 
93  // Print Taylor series solution
94  std::cout << "Taylor series solution = " << std::endl
95  << x << std::endl;
96 
97  // Print Taylor series solution derivative
98  std::cout << "Taylor series solution derivative= " << std::endl
99  << dxdx0 << std::endl;
100 
101  // Compute Taylor series expansion of solution x(t) = tan(t+pi/4)
102  double pi = std::atan(1.0)*4.0;
104  t.fastAccessCoeff(0) = pi/4.0;
105  t.fastAccessCoeff(1) = 1.0;
107 
108  // Compute Taylor series expansion of derivative
109  Sacado::Tay::Taylor<double> dudx0 = 0.5*(1.0+u*u);
110 
111  // Print expansion of solution
112  std::cout << "Exact expansion = " << std::endl
113  << u << std::endl;
114 
115  // Print expansion of solution
116  std::cout << "Exact expansion derivative = " << std::endl
117  << dudx0 << std::endl;
118 
119  // Compute maximum relative error
120  double max_err = 0.0;
121  double err = 0.0;
122  for (int k=0; k<=deg; k++) {
123  err = std::fabs(x.coeff(k) - u.coeff(k)) / (1.0 + fabs(u.coeff(k)));
124  if (err > max_err) max_err = err;
125  }
126  std::cout << "Maximum relative error = " << max_err << std::endl;
127 
128  // Compute maximum derivative relative error
129  double deriv_max_err = 0.0;
130  double deriv_err = 0.0;
131  for (int k=0; k<=deg; k++) {
132  deriv_err = std::fabs(dxdx0.coeff(k) - dudx0.coeff(k)) /
133  (1.0 + fabs(dudx0.coeff(k)));
134  if (deriv_err > deriv_max_err) deriv_max_err = deriv_err;
135  }
136  std::cout << "Maximum derivative relative error = " << deriv_max_err
137  << std::endl;
138 
139  double tol = 1.0e-12;
140  if ((max_err < tol) && (deriv_max_err < tol)) {
141  std::cout << "\nExample passed!" << std::endl;
142  return 0;
143  }
144  else {
145  std::cout <<"\nSomething is wrong, example failed!" << std::endl;
146  return 1;
147  }
148 
149  return 0;
150 }
151 
152 
153 
void resize(int d, bool keep_coeffs)
Resize polynomial to degree d.
void f()
atan(expr.val())
T & fastAccessCoeff(int i)
Returns degree i term without bounds checking.
int main()
Definition: ad_example.cpp:191
tan(expr.val())
const T * coeff() const
Returns Taylor coefficient array.
void reserve(int d)
Reserve space for a degree d polynomial.
const double tol
const T func(int n, T *x)
Definition: ad_example.cpp:49
fabs(expr.val())