Sacado Package Browser (Single Doxygen Collection) Version of the Day
taylor_example.cpp
Go to the documentation of this file.
00001 // $Id$ 
00002 // $Source$ 
00003 // @HEADER
00004 // ***********************************************************************
00005 // 
00006 //                           Sacado Package
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
00015 // License, or (at your option) any later version.
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 // ***********************************************************************
00030 // @HEADER
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) = 1.0;
00042 //
00043 //     The exact solution is x(t) = tan(t + pi/4)
00044 
00045 #include <iostream>
00046 
00047 #include "Sacado.hpp"
00048 
00049 // Function implementing RHS of ODE
00050 template <typename ScalarT>
00051 void func(ScalarT& f, const ScalarT& x) {
00052   f = 1.0 + x*x;
00053 }
00054 
00055 int main(int argc, char **argv)
00056 {
00057   double x0 = 1.0;                      // Initial condition
00058   unsigned int deg = 40;                // Degree of Taylor series solution
00059 
00060   Sacado::Tay::Taylor<double> x = x0;   // Taylor polynomial for independent
00061   Sacado::Tay::Taylor<double> f;        // Taylor polynomial for dependent
00062 
00063   // Reserve space for degree deg coefficients
00064   x.reserve(deg);
00065 
00066   // Compute Taylor series solution to dx/dt = f(x)
00067   for (unsigned int k=0; k<deg; k++) {
00068     func(f, x);
00069 
00070     // Set next coefficient
00071     x.resize(k+1, true);
00072 
00073     // x_{k+1} = f_k / (k+1)
00074     x.fastAccessCoeff(k+1) = f.coeff(k) / (k+1);
00075   }
00076 
00077   // Print Taylor series solution
00078   std::cout << "Taylor series solution = " << std::endl 
00079       << x << std::endl;
00080 
00081   // Compute Taylor series expansion of solution x(t) = tan(t+pi/4)
00082   double pi = std::atan(1.0)*4.0;
00083   Sacado::Tay::Taylor<double> t(deg);
00084   t.fastAccessCoeff(0) = pi/4.0;
00085   t.fastAccessCoeff(1) = 1.0;
00086   Sacado::Tay::Taylor<double> u = std::tan(t);
00087 
00088   // Print expansion of solution
00089   std::cout << "Exact expansion = " << std::endl
00090       << u << std::endl;
00091 
00092   // Compute maximum relative error
00093   double max_err = 0.0;
00094   double err = 0.0;
00095   for (unsigned int k=0; k<=deg; k++) {
00096     err = std::fabs(x.coeff(k) - u.coeff(k)) / (1.0 + fabs(u.coeff(k)));
00097     if (err > max_err) max_err = err;
00098   }
00099   std::cout << "Maximum relative error = " << max_err << std::endl;
00100 
00101   double tol = 1.0e-12;
00102   if (max_err < tol){
00103     std::cout << "\nExample passed!" << std::endl;
00104     return 0;
00105   }
00106   else {
00107     std::cout <<"\nSomething is wrong, example failed!" << std::endl;
00108     return 1;
00109   }
00110 }
00111 
00112   
00113 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines