00001 00003 // 00004 // -= ODE solver using Heun's Method =- 00005 // 00006 // Definition: "a 2nd order ordinary differential equation solver" 00007 // 00009 // 00010 // $RCSfile: HeunODEsolver.h,v $ 00011 // $Date: 2002/01/10 03:10:47 $ 00012 // $Revision: 1.4 $ 00013 // Copyright (C) 1998, 1999, 2000 Kevin Meinert, kevin@vrsource.org 00014 // 00015 // This library is free software; you can redistribute it and/or 00016 // modify it under the terms of the GNU Library General Public 00017 // License as published by the Free Software Foundation; either 00018 // version 2 of the License, or (at your option) any later version. 00019 // 00020 // This library is distributed in the hope that it will be useful, 00021 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00022 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00023 // Library General Public License for more details. 00024 // 00025 // You should have received a copy of the GNU Library General Public 00026 // License along with this library; if not, write to the Free 00027 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 00028 // 00030 #ifndef SECONDORDER_RUNGEKUTTA_METHOD 00031 #define SECONDORDER_RUNGEKUTTA_METHOD 00032 00033 #include "ani/Dynamics/ODEsolver.h" 00034 00035 namespace ani 00036 { 00037 00038 //: the second-order Runge-Kutta method 00039 template <class _item> 00040 class HeunODEsolver : public ODEsolver<_item> 00041 { 00042 public: 00043 HeunODEsolver() {} 00044 virtual ~HeunODEsolver() {} 00045 00046 //: the second-order Runge-Kutta method 00047 // 00048 // give: 00049 // currentState = current state of the item 00050 // currentTime = current time of the item 00051 // timeDelta = current stepsize (t(n) - t(n+1)) == time delta 00052 // returns: 00053 // currentState = x(t0 + h) == the next state of the item after taking this step... 00054 virtual void exec( _item& currentState, float timeDelta ) 00055 { 00056 // compute next state. 00057 00058 // any particle that is put into f() needs to have its mass and force accums set 00059 // to equal the current _item, since all the operations only act on the phase space 00060 // (pos, vel), this is needed. 00061 temp.copy( currentState ); 00062 00063 // k1 = h * func( x0, t0 ); 00064 // f = func(..., ...) 00065 f.computeDerivative( currentState, timeDelta * 0.0f ); 00066 // k1 = h * f 00067 k1.multiplyPhase( f, timeDelta ); 00068 00069 //: k2 = h * func( x0 + k1, t0 + h ); 00070 // f = func(..., ...) 00071 temp.addPhase( currentState, k1 ); 00072 f.computeDerivative( temp, timeDelta * 1.0f ); 00073 // k4 = h * f 00074 k2.multiplyPhase( f, timeDelta ); 00075 00076 // xt0h = x0 + 1/2 * (k1 + k2); 00077 //currentState = currentState + 1/2 * (k1 + k2); 00078 changeInState.addPhase( k1, k2 ); 00079 changeInState.multiplyPhase( 0.5f ); 00080 00081 currentState.addPhase( changeInState ); 00082 currentState.normalize(); 00083 } 00084 00085 // create the two gradient or ``k'' terms 00086 _item k1, k2; 00087 _item changeInState, f, temp; 00088 }; 00089 }; 00090 00091 #endif