diff options
| -rw-r--r-- | TODO | 2 | ||||
| -rw-r--r-- | src/vZ.hpp | 1 | ||||
| -rw-r--r-- | src/vZ/EquationSystem.hpp | 44 | ||||
| -rw-r--r-- | src/vZ/Vector.hpp | 228 | ||||
| -rw-r--r-- | tests/Makefile.am | 2 | ||||
| -rw-r--r-- | tests/Vector.cpp | 40 |
6 files changed, 304 insertions, 13 deletions
@@ -1,2 +1,2 @@ -- Implement vector class - Implement FSAL +- Add C++0x API extensions @@ -24,6 +24,7 @@ #define VZ_HPP #include <vZ/Traits.hpp> +#include <vZ/Vector.hpp> #include <vZ/EquationSystem.hpp> #include <vZ/Integrator.hpp> #include <vZ/RK.hpp> diff --git a/src/vZ/EquationSystem.hpp b/src/vZ/EquationSystem.hpp index 790676c..5a50f6d 100644 --- a/src/vZ/EquationSystem.hpp +++ b/src/vZ/EquationSystem.hpp @@ -64,6 +64,38 @@ namespace vZ Traits(); }; + // Unary operators + + template <std::size_t N, typename T> + inline EquationSystem<N, T> + operator+(const EquationSystem<N, T>& rhs) + { + return rhs; + } + + template <std::size_t N, typename T> + inline EquationSystem<N, T> + operator-(const EquationSystem<N, T>& rhs) + { + EquationSystem<N, T> res; + for (std::size_t i = 0; i < N; ++i) { + res[i] = -rhs[i]; + } + return res; + } + + template <std::size_t N, typename T> + typename EquationSystem<N, T>::Scalar + abs(const EquationSystem<N, T>& es) + { + typename EquationSystem<N, T>::Scalar ret(0); + for (std::size_t i = 0; i < N; ++i) { + using std::abs; + ret = std::max(ret, abs(es[i])); + } + return ret; + } + // Binary operators template <std::size_t N, typename T> @@ -155,18 +187,6 @@ namespace vZ } return *this; } - - template <std::size_t N, typename T> - typename EquationSystem<N, T>::Scalar - abs(const EquationSystem<N, T>& es) - { - typename EquationSystem<N, T>::Scalar ret(0); - for (std::size_t i = 0; i < N; ++i) { - using std::abs; - ret = std::max(ret, abs(es[i])); - } - return ret; - } } #endif // VZ_EQUATIONSYSTEM_HPP diff --git a/src/vZ/Vector.hpp b/src/vZ/Vector.hpp new file mode 100644 index 0000000..cf3bf0a --- /dev/null +++ b/src/vZ/Vector.hpp @@ -0,0 +1,228 @@ +/************************************************************************* + * 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_VECTOR_HPP +#define VZ_VECTOR_HPP + +#include <algorithm> +#include <cmath> +#include <cstddef> + +namespace vZ +{ + // An N-dimensional vector + template <std::size_t N, typename T = double> + class Vector + { + public: + typedef typename Traits<T>::Scalar Scalar; + + Vector() { } + explicit Vector(T x) { m_values[0] = x; } + Vector(T x, T y) { m_values[0] = x; m_values[1] = y; } + Vector(T x, T y, T z) { m_values[0] = x; m_values[1] = y; m_values[2] = z; } + // Vector(const Vector& v); + // ~Vector(); + + // Vector& operator=(const Vector& v); + + // Component access + + T& operator[](std::size_t i) { return m_values[i]; } + const T& operator[](std::size_t i) const { return m_values[i]; } + + T x() const { return m_values[0]; } + T y() const { return m_values[1]; } + T z() const { return m_values[2]; } + + // Operators + inline Vector& operator+=(const Vector& rhs); + inline Vector& operator-=(const Vector& rhs); + inline Vector& operator*=(Scalar rhs); + inline Vector& operator/=(Scalar rhs); + + private: + T m_values[N]; + }; + + // Disallow 0-sized Vectors + template <typename T> + class Vector<0, T>; + + // Traits specialization + template <std::size_t N, typename T> + class Traits<Vector<N, T> > + { + public: + typedef typename Traits<T>::Scalar Scalar; + + private: + Traits(); + }; + + // Unary operators + + template <std::size_t N, typename T> + inline Vector<N, T> + operator+(const Vector<N, T>& rhs) + { + return rhs; + } + + template <std::size_t N, typename T> + inline Vector<N, T> + operator-(const Vector<N, T>& rhs) + { + Vector<N, T> res; + for (std::size_t i = 0; i < N; ++i) { + res[i] = -rhs[i]; + } + return res; + } + + template <std::size_t N, typename T> + inline typename Vector<N, T>::Scalar + norm(const Vector<N, T>& v) + { + using std::sqrt; + return sqrt(dot(v, v)); + } + + template <std::size_t N, typename T> + inline typename Vector<N, T>::Scalar + abs(const Vector<N, T>& v) + { + return norm(v); + } + + // Binary operators + + template <std::size_t N, typename T> + inline Vector<N, T> + operator+(const Vector<N, T>& lhs, const Vector<N, T>& rhs) + { + Vector<N, T> res = lhs; + res += rhs; + return res; + } + + template <std::size_t N, typename T> + inline Vector<N, T> + operator-(const Vector<N, T>& lhs, const Vector<N, T>& rhs) + { + Vector<N, T> res = lhs; + res -= rhs; + return res; + } + + template <std::size_t N, typename T> + inline Vector<N, T> + operator*(typename Vector<N, T>::Scalar lhs, + const Vector<N, T>& rhs) + { + Vector<N, T> res = rhs; + res *= lhs; + return res; + } + + template <std::size_t N, typename T> + inline Vector<N, T> + operator*(const Vector<N, T>& lhs, + typename Vector<N, T>::Scalar rhs) + { + Vector<N, T> res = lhs; + res *= rhs; + return res; + } + + template <std::size_t N, typename T> + inline Vector<N, T> + operator/(const Vector<N, T>& lhs, + typename Vector<N, T>::Scalar rhs) + { + Vector<N, T> res = lhs; + res /= rhs; + return res; + } + + template <std::size_t N, typename T> + inline typename Vector<N, T>::Scalar + dot(const Vector<N, T>& lhs, const Vector<N, T>& rhs) + { + typename Vector<N, T>::Scalar res(0); + for (std::size_t i = 0; i < N; ++i) { + res += lhs[i]*rhs[i]; + } + return res; + } + + template <typename T> + inline Vector<3, T> + cross(const Vector<3, T>& lhs, const Vector<3, T>& rhs) + { + return Vector<3, T>(lhs.y()*rhs.z() - lhs.z()*rhs.y(), + lhs.z()*rhs.x() - lhs.x()*rhs.z(), + lhs.x()*rhs.y() - lhs.y()*rhs.x()); + } + + // Implementation + + template <std::size_t N, typename T> + inline Vector<N, T>& + Vector<N, T>::operator+=(const Vector<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> + inline Vector<N, T>& + Vector<N, T>::operator-=(const Vector<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> + inline Vector<N, T>& + Vector<N, T>::operator*=(typename Vector<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> + inline Vector<N, T>& + Vector<N, T>::operator/=(typename Vector<N, T>::Scalar rhs) + { + for (std::size_t i = 0; i < N; ++i) { + m_values[i] /= rhs; + } + return *this; + } +} + +#endif // VZ_VECTOR_HPP diff --git a/tests/Makefile.am b/tests/Makefile.am index 1cd184d..745562c 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -28,6 +28,7 @@ check_PROGRAMS = Euler-test \ RKF45-test \ CK45-test \ DP45-test \ + Vector-test \ EquationSystem-test TESTS = $(check_PROGRAMS) @@ -40,4 +41,5 @@ BS23_test_SOURCES = BS23.cpp RKF45_test_SOURCES = RKF45.cpp CK45_test_SOURCES = CK45.cpp DP45_test_SOURCES = DP45.cpp +Vector_test_SOURCES = Vector.cpp EquationSystem_test_SOURCES = EquationSystem.cpp diff --git a/tests/Vector.cpp b/tests/Vector.cpp new file mode 100644 index 0000000..7f53352 --- /dev/null +++ b/tests/Vector.cpp @@ -0,0 +1,40 @@ +#include "vZ.hpp" +#include <cmath> +#include <cstdlib> +#include <iostream> +#include <iomanip> + +// y' = -y (y == C*exp(-t)) +vZ::Vector<3> +f(double t, vZ::Vector<3> y) +{ + return -y; +} + +int +main() +{ + vZ::GenericDP45Integrator<vZ::Vector<3> > integrator(f); + integrator.tol(1e-6).y(vZ::Vector<3>(1.0, 1.0, 1.0)).x(0.0).h(0.06); + + integrator.integrate(2.0); + + double actual = integrator.y().x(); + double expected = std::exp(-2.0); + + std::cout << std::setprecision(10) + << "Numerical: " << actual << std::endl + << "Expected: " << expected << std::endl + << "h: " << integrator.h() << std::endl + << "Iterations: " << integrator.iterations() << std::endl + << "Rejections: " << integrator.rejections() << std::endl; + + double error = std::fabs(expected - actual)/expected; + if (error > 1.5e-6) { + std::cerr << "Error: " << 100.0*error << "%" << std::endl; + return EXIT_FAILURE; + } else { + std::cout << "Error: " << 100.0*error << "%" << std::endl; + return EXIT_SUCCESS; + } +} |