Sacado Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
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
53template <typename ScalarT>
54void func(ScalarT& f, const ScalarT& x) {
55 f = 1.0 + x*x;
56}
57
58int 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);
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;
106 Sacado::Tay::Taylor<double> u = std::tan(t);
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
fabs(expr.val())
int main()
Fad specializations for Teuchos::BLAS wrappers.
Taylor polynomial class.
void reserve(int d)
Reserve space for a degree d polynomial.
T & fastAccessCoeff(int i)
Returns degree i term without bounds checking.
const T * coeff() const
Returns Taylor coefficient array.
void func(ScalarT &f, const ScalarT &x)
const double tol