diff options
| -rw-r--r-- | TODO | 2 | ||||
| -rw-r--r-- | src/Makefile.am | 1 | ||||
| -rw-r--r-- | src/vZ.hpp | 1 | ||||
| -rw-r--r-- | src/vZ/EquationSystem.hpp | 154 | ||||
| -rw-r--r-- | src/vZ/Integrator.hpp | 4 | ||||
| -rw-r--r-- | tests/EquationSystem.cpp | 46 | ||||
| -rw-r--r-- | tests/Makefile.am | 22 |
7 files changed, 217 insertions, 13 deletions
@@ -1,2 +1,2 @@ -- Support systems of ODEs +- Implement vector class - Implement FSAL diff --git a/src/Makefile.am b/src/Makefile.am index bf2b10c..7db29dc 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -23,6 +23,7 @@ nobase_include_HEADERS = vZ.hpp \ vZ/CK45.hpp \ vZ/DP45.hpp \ vZ/Euler.hpp \ + vZ/EquationSystem.hpp \ vZ/HE12.hpp \ vZ/Heun.hpp \ vZ/Integrator.hpp \ @@ -24,6 +24,7 @@ #define VZ_HPP #include <vZ/Traits.hpp> +#include <vZ/EquationSystem.hpp> #include <vZ/Integrator.hpp> #include <vZ/RK.hpp> #include <vZ/Simple.hpp> diff --git a/src/vZ/EquationSystem.hpp b/src/vZ/EquationSystem.hpp new file mode 100644 index 0000000..0963baf --- /dev/null +++ b/src/vZ/EquationSystem.hpp @@ -0,0 +1,154 @@ +/************************************************************************* + * Copyright (C) 2009-2010 Tavian Barnes <tavianator@gmail.com> * + * * + * This file is part of The vZ Library. * + * * + * The vZ Library is free software; you can redistribute it and/or * + * modify it under the terms of the GNU Lesser General Public License as * + * published by the Free Software Foundation; either version 3 of the * + * License, or (at your option) any later version. * + * * + * The vZ Library is distributed in the hope that it will be useful, but * + * WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * + * Lesser General Public License for more details. * + * * + * You should have received a copy of the GNU Lesser General Public * + * License along with this program. If not, see * + * <http://www.gnu.org/licenses/>. * + *************************************************************************/ + +#ifndef VZ_EQUATIONSYSTEM_HPP +#define VZ_EQUATIONSYSTEM_HPP + +#include <cstddef> + +namespace vZ +{ + // A class to easily represent a system of ODEs + template <std::size_t N, typename T = double> + class EquationSystem + { + public: + typedef typename Traits<T>::Scalar Scalar; + + // EquationSystem(); + // ~EquationSystem(); + + T& operator[](std::size_t i) { return m_values[0]; } + const T& operator[](std::size_t i) const { return m_values[0]; } + + EquationSystem& operator+=(const EquationSystem& rhs); + EquationSystem& operator-=(const EquationSystem& rhs); + EquationSystem& operator*=(Scalar rhs); + EquationSystem& operator/=(Scalar rhs); + + private: + T m_values[N]; + }; + + // Traits specialization + template <std::size_t N, typename T> + class Traits<EquationSystem<N, T> > + { + public: + typedef typename Traits<T>::Scalar Scalar; + + private: + Traits(); + }; + + // Binary operators + + template <std::size_t N, typename T> + inline EquationSystem<N, T> + operator+(const EquationSystem<N, T>& lhs, const EquationSystem<N, T>& rhs) + { + EquationSystem<N, T> res = lhs; + res += rhs; + return res; + } + + template <std::size_t N, typename T> + inline EquationSystem<N, T> + operator-(const EquationSystem<N, T>& lhs, const EquationSystem<N, T>& rhs) + { + EquationSystem<N, T> res = lhs; + res -= rhs; + return res; + } + + template <std::size_t N, typename T> + inline EquationSystem<N, T> + operator*(typename EquationSystem<N, T>::Scalar lhs, + const EquationSystem<N, T>& rhs) + { + EquationSystem<N, T> res = rhs; + res *= lhs; + return res; + } + + template <std::size_t N, typename T> + inline EquationSystem<N, T> + operator*(const EquationSystem<N, T>& lhs, + typename EquationSystem<N, T>::Scalar rhs) + { + EquationSystem<N, T> res = lhs; + res *= rhs; + return res; + } + + template <std::size_t N, typename T> + inline EquationSystem<N, T> + operator/(const EquationSystem<N, T>& lhs, + typename EquationSystem<N, T>::Scalar rhs) + { + EquationSystem<N, T> res = lhs; + res /= rhs; + return res; + } + + // Implementation + + template <std::size_t N, typename T> + EquationSystem<N, T>& + EquationSystem<N, T>::operator+=(const EquationSystem<N, T>& rhs) + { + for (std::size_t i = 0; i < N; ++i) { + m_values[i] += rhs.m_values[i]; + } + return *this; + } + + template <std::size_t N, typename T> + EquationSystem<N, T>& + EquationSystem<N, T>::operator-=(const EquationSystem<N, T>& rhs) + { + for (std::size_t i = 0; i < N; ++i) { + m_values[i] -= rhs.m_values[i]; + } + return *this; + } + + template <std::size_t N, typename T> + EquationSystem<N, T>& + EquationSystem<N, T>::operator*=(typename EquationSystem<N, T>::Scalar rhs) + { + for (std::size_t i = 0; i < N; ++i) { + m_values[i] *= rhs; + } + return *this; + } + + template <std::size_t N, typename T> + EquationSystem<N, T>& + EquationSystem<N, T>::operator/=(typename EquationSystem<N, T>::Scalar rhs) + { + for (std::size_t i = 0; i < N; ++i) { + m_values[i] /= rhs; + } + return *this; + } +} + +#endif // VZ_EQUATIONSYSTEM_HPP diff --git a/src/vZ/Integrator.hpp b/src/vZ/Integrator.hpp index c6d6fbe..cb9649a 100644 --- a/src/vZ/Integrator.hpp +++ b/src/vZ/Integrator.hpp @@ -39,9 +39,9 @@ namespace vZ typedef typename Traits<Y>::Scalar Scalar; typedef std::tr1::function<Y (Scalar, Y)> Function; - // By default, y and t start at zero, h starts UNDEFINED + // By default, y, t, and h start UNDEFINED GenericIntegrator(Function f) - : m_f(f), m_y(0), m_x(0), m_iterations(0) { } + : m_f(f),m_iterations(0) { } virtual ~GenericIntegrator() { } GenericIntegrator& y(Y y) { m_y = y; return *this; } diff --git a/tests/EquationSystem.cpp b/tests/EquationSystem.cpp new file mode 100644 index 0000000..973dd30 --- /dev/null +++ b/tests/EquationSystem.cpp @@ -0,0 +1,46 @@ +#include "vZ.hpp" +#include <cmath> +#include <cstdlib> +#include <iostream> +#include <iomanip> + +typedef vZ::EquationSystem<2> Y; + +// y'' = y (y == C*exp(t)) +Y +f(double t, Y y) +{ + Y r; + r[0] = y[1]; + r[1] = y[0]; + return r; +} + +int +main() +{ + Y y; + y[0] = 1.0; + y[1] = 1.0; + vZ::GenericEulerIntegrator<Y> integrator(f); + integrator.y(y).x(0.0).h(0.01); + + integrator.integrate(2.0); + + double actual = integrator.y()[0]; + double expected = std::exp(2.0); + + std::cout << std::setprecision(10) + << "Numerical: " << actual << std::endl + << "Expected: " << expected << std::endl + << "Iterations: " << integrator.iterations() << std::endl; + + double error = std::fabs(expected - actual)/expected; + if (error > 0.01) { + std::cerr << "Error: " << 100.0*error << "%" << std::endl; + return EXIT_FAILURE; + } else { + std::cout << "Error: " << 100.0*error << "%" << std::endl; + return EXIT_SUCCESS; + } +} diff --git a/tests/Makefile.am b/tests/Makefile.am index 5874419..1cd184d 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -27,15 +27,17 @@ check_PROGRAMS = Euler-test \ BS23-test \ RKF45-test \ CK45-test \ - DP45-test + DP45-test \ + EquationSystem-test TESTS = $(check_PROGRAMS) -Euler_test_SOURCES = Euler.cpp -Midpoint_test_SOURCES = Midpoint.cpp -Heun_test_SOURCES = Heun.cpp -RK4_test_SOURCES = RK4.cpp -HE12_test_SOURCES = HE12.cpp -BS23_test_SOURCES = BS23.cpp -RKF45_test_SOURCES = RKF45.cpp -CK45_test_SOURCES = CK45.cpp -DP45_test_SOURCES = DP45.cpp +Euler_test_SOURCES = Euler.cpp +Midpoint_test_SOURCES = Midpoint.cpp +Heun_test_SOURCES = Heun.cpp +RK4_test_SOURCES = RK4.cpp +HE12_test_SOURCES = HE12.cpp +BS23_test_SOURCES = BS23.cpp +RKF45_test_SOURCES = RKF45.cpp +CK45_test_SOURCES = CK45.cpp +DP45_test_SOURCES = DP45.cpp +EquationSystem_test_SOURCES = EquationSystem.cpp |