From 5233c01a26329306c9f1d08c1a39733aee2cc518 Mon Sep 17 00:00:00 2001 From: Tim Blechmann Date: Tue, 21 Dec 2004 09:22:10 +0000 Subject: This commit was generated by cvs2svn to compensate for changes in r2423, which included commits to RCS files with non-trunk default branches. svn path=/trunk/externals/tb/; revision=2424 --- chaos/src/ode_base.cpp | 90 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 90 insertions(+) create mode 100644 chaos/src/ode_base.cpp (limited to 'chaos/src/ode_base.cpp') diff --git a/chaos/src/ode_base.cpp b/chaos/src/ode_base.cpp new file mode 100644 index 0000000..cfedeee --- /dev/null +++ b/chaos/src/ode_base.cpp @@ -0,0 +1,90 @@ +// +// +// chaos~ +// Copyright (C) 2004 Tim Blechmann +// +// This program is free software; you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation; either version 2 of the License, or +// (at your option) any later version. +// +// This program 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 General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program; see the file COPYING. If not, write to +// the Free Software Foundation, Inc., 59 Temple Place - Suite 330, +// Boston, MA 02111-1307, USA. + +#include "ode_base.hpp" + +void ode_base::rk1() +{ + m_system (m_k[0], m_data); + for (int i = 0; i != m_num_eq; ++i) + m_data[i] += m_dt * m_k[0][i]; +} + + +void ode_base::rk2() +{ + m_system (m_k[0], m_data); + for (int i = 0; i != m_num_eq; ++i) + m_k[0][i] = m_k[0][i] * 0.5 * m_dt + m_data[i]; + + m_system (m_k[1], m_k[0]); + for (int i = 0; i != m_num_eq; ++i) + m_data[i] += m_dt * m_k[1][i]; +} + + +void ode_base::rk4() +{ + m_system (m_k[0], m_data); + for (int i = 0; i != m_num_eq; ++i) + { + m_k[0][i] *= m_dt; + m_tmp[i] = m_data[i] + 0.5 * m_k[0][i]; + } + + m_system (m_k[1], m_tmp); + for (int i = 0; i != m_num_eq; ++i) + { + m_k[1][i] *= m_dt; + m_tmp[i] = m_data[i] + 0.5 * m_k[1][i]; + } + + m_system (m_k[2], m_tmp); + for (int i = 0; i != m_num_eq; ++i) + { + m_k[2][i] *= m_dt; + m_tmp[i] = m_data[i] + m_k[2][i]; + } + + m_system (m_k[3], m_tmp); + for (int i = 0; i != m_num_eq; ++i) + m_k[3][i] *= m_dt; + + for (int i = 0; i != m_num_eq; ++i) + m_data[i] += (m_k[0][i] + 2. * (m_k[1][i] + m_k[2][i]) + m_k[3][i]) + / 6.; +} + +void ode_base::m_step() +{ + switch (m_method) + { + case 0: + rk1(); + break; + case 1: + rk2(); + break; + case 2: + rk3(); + break; + } + +} -- cgit v1.2.1