00001 #ifndef GSGL_MATH_SOLVER_H 00002 #define GSGL_MATH_SOLVER_H 00003 00004 // 00005 // $Id: solver.hpp 2 2008-03-01 20:58:50Z kulibali $ 00006 // 00007 // Copyright (c) 2008, The Periapsis Project. All rights reserved. 00008 // 00009 // Redistribution and use in source and binary forms, with or without 00010 // modification, are permitted provided that the following conditions are 00011 // met: 00012 // 00013 // * Redistributions of source code must retain the above copyright notice, 00014 // this list of conditions and the following disclaimer. 00015 // 00016 // * Redistributions in binary form must reproduce the above copyright 00017 // notice, this list of conditions and the following disclaimer in the 00018 // documentation and/or other materials provided with the distribution. 00019 // 00020 // * Neither the name of the The Periapsis Project nor the names of its 00021 // contributors may be used to endorse or promote products derived from 00022 // this software without specific prior written permission. 00023 // 00024 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS 00025 // IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED 00026 // TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A 00027 // PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER 00028 // OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 00029 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 00030 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 00031 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 00032 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 00033 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 00034 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 00035 // 00036 00037 #include "math/math.hpp" 00038 00039 namespace gsgl 00040 { 00041 00042 namespace math 00043 { 00044 00045 00046 /// Base class for numerical equation solvers, where the equation is x' = f(t, x). 00047 /// State classes need to implement operator+(state) const, operator*(double) const, and derivative(double) const. 00048 template <typename S> 00049 class solver 00050 { 00051 public: 00052 virtual S next(const S & x, const double & t, const double & dt) = 0; 00053 }; // class solver 00054 00055 00056 /// Euler's method solver. 00057 template <typename S> 00058 class euler_solver 00059 : public solver<S> 00060 { 00061 public: 00062 S next(const S & x, const double & t, const double & dt) 00063 { 00064 return x + x.derivative(t) * dt; 00065 } 00066 }; // class euler_solver 00067 00068 00069 /// The venerable Runge-Kutta solver. 00070 template <typename S> 00071 class runge_kutta_solver 00072 : public solver<S> 00073 { 00074 public: 00075 S next(const S & x, const double & t, const double & dt) 00076 { 00077 S k1 = x.derivative(t); 00078 S k2 = (x + k1*(dt*0.5)).derivative(t + dt*0.5); 00079 S k3 = (x + k2*(dt*0.5)).derivative(t + dt*0.5); 00080 S k4 = (x + k3*dt).derivative(t + dt); 00081 00082 return x + (k1 + k2*2.0 + k3*2.0 + k4)*(dt/6.0); 00083 } // next() 00084 }; // class runge_kutta_solver 00085 00086 00087 } // namespace math 00088 00089 } // namespace math 00090 00091 #endif
1.5.5