Boost logo

Boost :

From: Harry Erwin (harry.erwin_at_[hidden])
Date: 2002-10-24 05:29:34

>On Mittwoch, 23. Oktober 2002 21:52, Matthias Troyer wrote:
>> I am confused about the scope of a numerical integration library. It
>> seems that three things are being mixed up here:
>> i) 1D integration of a function with x-points determined by the
>> algorithm (e.g. Romberg integration)
>> ii) Integration of first order ODEs by some algorithm (e.g.
>> Runge-Kutta)
>> iii) Forward Euler integration of an ODE with fixed step witdh delta_x
>> and values of the derivative, given by some external function.
>> While I agree that both i) and ii) would be interesting libraries, the
>> usage case iii) which seems much discussed does not warrant a separate
>> library in my opinion.
>Agreed. Let me add that I think the interface to i) is more or less
>straightforward, while the interface to ii) is much more complex and less
>Dr. Martin Weiser Zuse Institute Berlin
>weiser_at_[hidden] Scientific Computing
> Numerical Analysis and Modelling
>Unsubscribe & other changes:

Significant portions of the sample header files. The control type is
a functor. The state type provides a state vector. The fourth-order
general, adaptive step-size, and Bulirsch-Stoer algorithms based on
the yellow book, so I'm leaving them out.

const state rungekutta(double deltat, state& u, control& c, int order);
// arbitrary order

const state rk4(double deltat, state& u, state& dvdt, control& c);
// fourth-order general

const state rkqc(const double htry, const state& u,
        const state& dudt, control& c, const double eps,
        const state& yscale, double& hdid, double& hnext);
// adaptive step-size 5th order RK

const state bsstep(const double htry, const state& u, const state& dvdt,
        const double eps, const state& yscale, double& hdid, double& hnext);

class state {
        valarray<double> p; // tvector<T> p; // position
        valarray<double> dp; // tvector<T> dp; // velocity
        double t; // time
        double m; // mass

        state(void): p(0.0,3), dp(0.0,3), t(0.0),
        m(0.0) {} // default initialization

        state(const valarray<double>& pos, const valarray<double>& vel,
                double time, double mass):
                p(pos), dp(vel), t(time), m(mass) {}

        state(double x1, double y1, double z1, double dx1, double
dy1, double dz1,
                double time, double mass);

        state(const state & source);

        ~state() {;}

        friend ostream& operator<<(ostream& s, state& v);
        friend ostream& operator>>(istream& s, state& v);
           friend state operator + ( const state&, const state& );
           friend state operator - ( const state&, const state& );
           friend state operator * ( const state&, const double );
           friend state operator / ( const state&, const double );
           friend state operator * ( const double, const state & );
           friend double operator * (const state &, const state &);
           friend state operator / (const state &, const state &);
        friend bool operator != ( const state&, const state& );
           friend bool operator == ( const state&, const state& );

        state& operator = ( const state& rold){
                if(this==&rold) return *this;
                state temp(rold);
                return *this;

    void operator += ( const state& );
    void operator -= ( const state& );
    void operator *= ( const double );
    void operator /= ( const double );
    void operator - ();
    void swap(state& b)
                   std::swap(p, b.p);
                   std::swap(dp, b.dp);
                   std::swap(t, b.t);
                   std::swap(m, b.m);
        void print(ostream& s);

class control { // wrappable, not a wrapper
// note this defines a functor
        size_t sn;
        static size_t snm;
        static size_t count;

        object* object_; // the parent object

        control(): sn(++snm), object_(0) { ++count; }
        control(const control& c): sn(++snm), object_(c.object_) { ++count; }
        control(object* the_object): sn(++snm), object_(the_object) {
++count; }

        control& operator=(const control& c)
                if(this == &c)return *this;
                control temp(c);
                return *this;

        void swap(control& c)
                std::swap(object_, c.object_);

        void set_object(object* an_object) { object_ = an_object;}

        object* get_object() { return object_; }

        virtual ~control() {--count;} // note object_ is not owned

        virtual state operator()(state& t){
                // the function call. returns acceleration
                state temp (0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
                return temp; // override!
        virtual control* clone() const { return new control(*this); }
        bool operator==(const control& c) const {return (;}
        bool operator<(const control& c) const {return (sn<;}
        virtual void perform_update(const state& where, const double ret_time,
                const valarray<double>&cry_position, const
valarray<double>&cry_orientation, const double cry_time,
                const valarray<double>&ret_position, const
        virtual void perform_update(const afference& datum);
        virtual void next_action();

The following is not compatible with the forgoing, but provides
insight into why I defined state the way I did.

const state rungekutta(double deltat, state& u, const intent& in,
control& c, int order) // arbitrary order
        double alphadt;
        valarray<double> dpos(0.0,3);
        valarray<double> dvel(0.0,3);
        double dtime;
        double dmass;
        state interu = state(u.p, u.dp, u.t, u.m,, u.plan, u.plan_start);
        for(int i=0;i<order;i++) {
                alphadt = deltat/(order - i);
                state deltau = c(interu, in);
                dpos = (deltau.p)*alphadt;
                dvel = (deltau.dp)*alphadt;
                dtime = alphadt;
                dmass = 0.0;
                interu = state((u.p+dpos), (u.dp+dvel), (u.t+dtime),
                        (, (u.plan), (u.plan_start);
        state res = state(dpos,dvel,dtime,dmass, (, (u.plan),
        return res;

Harry Erwin, PhD, Senior Lecturer of Computing, University of 
Sunderland. Computational neuroscientist modeling bat bioacoustics 
and behavior. <>

Boost list run by bdawes at, gregod at, cpdaniel at, john at