Go to the documentation of this file.
```00001 // \$Id\$
00002 // \$Source\$
00004 // ***********************************************************************
00005 //
00007 //                 Copyright (2006) Sandia Corporation
00008 //
00009 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00010 // the U.S. Government retains certain rights in this software.
00011 //
00012 // This library is free software; you can redistribute it and/or modify
00013 // it under the terms of the GNU Lesser General Public License as
00014 // published by the Free Software Foundation; either version 2.1 of the
00016 //
00017 // This library is distributed in the hope that it will be useful, but
00018 // WITHOUT ANY WARRANTY; without even the implied warranty of
00019 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00020 // Lesser General Public License for more details.
00021 //
00022 // You should have received a copy of the GNU Lesser General Public
00023 // License along with this library; if not, write to the Free Software
00024 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00025 // USA
00026 // Questions? Contact David M. Gay (dmgay@sandia.gov) or Eric T. Phipps
00027 // (etphipp@sandia.gov).
00028 //
00029 // ***********************************************************************
00031
00032 // taylor_example
00033 //
00034 //  usage:
00035 //     taylor_example
00036 //
00037 //  output:
00038 //     prints the results of computing a single Taylor series expansion of
00039 //     the solution to:
00040 //
00041 //           dx/dt = 1 + x^2,    x(0) = x0 = 1.0;
00042 //
00043 //     The exact solution is x(t) = tan(t + atan(x0)) = tan(t + pi/4)
00044 //     Also computes the derivative of the Taylor series solution with
00045 //     respect to x0.  The exact series derivative is
00046 //     dx/dx0(t) = 1/2 * sec^2(t + atan(x0))
00047
00048 #include <iostream>
00049
00051
00052 // Function implementing RHS of ODE
00053 template <typename ScalarT>
00054 void func(ScalarT& f, const ScalarT& x) {
00055   f = 1.0 + x*x;
00056 }
00057
00058 int main(int argc, char **argv)
00059 {
00060   double x0 = 1.0;                      // Initial condition
00061   unsigned int deg = 40;                // Degree of Taylor series solution
00062
00063   Sacado::Tay::Taylor<double> x = x0;   // Taylor polynomial for independent
00064   Sacado::Tay::Taylor<double> f;        // Taylor polynomial for dependent
00065
00066   // Reserve space for degree deg coefficients
00067   x.reserve(deg);
00068
00069   // Compute Taylor series solution to dx/dt = f(x)
00070   for (unsigned int k=0; k<deg; k++) {
00071     func(f, x);
00072
00073     // Set next coefficient
00074     x.resize(k+1, true);
00075
00076     // x_{k+1} = f_k / (k+1)
00077     x.fastAccessCoeff(k+1) = f.coeff(k) / (k+1);
00078   }
00079
00080   // Compute derivative w.r.t. x0
00085   dxdx0.fastAccessCoeff(0) = 1.0;
00086   for (unsigned int k=0; k<deg; k++) {
00087     dxdx0.fastAccessCoeff(k+1) = 0.0;
00088     for (unsigned int j=0; j<=k; j++)
00089       dxdx0.fastAccessCoeff(k+1) += f_fad.dx(0).coeff(k-j) * dxdx0.coeff(j);
00090     dxdx0.fastAccessCoeff(k+1) /= k+1;
00091   }
00092
00093   // Print Taylor series solution
00094   std::cout << "Taylor series solution = " << std::endl
00095       << x << std::endl;
00096
00097   // Print Taylor series solution derivative
00098   std::cout << "Taylor series solution derivative= " << std::endl
00099       << dxdx0 << std::endl;
00100
00101   // Compute Taylor series expansion of solution x(t) = tan(t+pi/4)
00102   double pi = std::atan(1.0)*4.0;
00104   t.fastAccessCoeff(0) = pi/4.0;
00105   t.fastAccessCoeff(1) = 1.0;
00107
00108   // Compute Taylor series expansion of derivative
00110
00111   // Print expansion of solution
00112   std::cout << "Exact expansion = " << std::endl
00113       << u << std::endl;
00114
00115   // Print expansion of solution
00116   std::cout << "Exact expansion derivative = " << std::endl
00117       << dudx0 << std::endl;
00118
00119   // Compute maximum relative error
00120   double max_err = 0.0;
00121   double err = 0.0;
00122   for (unsigned int k=0; k<=deg; k++) {
00123     err = std::fabs(x.coeff(k) - u.coeff(k)) / (1.0 + fabs(u.coeff(k)));
00124     if (err > max_err) max_err = err;
00125   }
00126   std::cout << "Maximum relative error = " << max_err << std::endl;
00127
00128   // Compute maximum derivative relative error
00129   double deriv_max_err = 0.0;
00130   double deriv_err = 0.0;
00131   for (unsigned int k=0; k<=deg; k++) {
00132     deriv_err = std::fabs(dxdx0.coeff(k) - dudx0.coeff(k)) /
00133       (1.0 + fabs(dudx0.coeff(k)));
00134     if (deriv_err > deriv_max_err) deriv_max_err = deriv_err;
00135   }
00136   std::cout << "Maximum derivative relative error = " << deriv_max_err
00137       << std::endl;
00138
00139   double tol = 1.0e-12;
00140   if ((max_err < tol) && (deriv_max_err < tol)) {
00141     std::cout << "\nExample passed!" << std::endl;
00142     return 0;
00143   }
00144   else {
00145     std::cout <<"\nSomething is wrong, example failed!" << std::endl;
00146     return 1;
00147   }
00148
00149   return 0;
00150 }
00151
00152
00153
```
Generated on Wed Apr 13 10:19:30 2011 for Sacado Package Browser (Single Doxygen Collection) by  1.6.3