Re: [boost] Re: Numerical Integration - continued

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
>obvious.
>
>Yours,
>Martin
>
>--
>Dr. Martin Weiser            Zuse Institute Berlin
>weiser_at_[hidden]                Scientific Computing
>http://www.zib.de/weiser     Numerical Analysis and Modelling
>
>_______________________________________________
>Unsubscribe & other changes: http://lists.boost.org/mailman/listinfo.cgi/boost

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  {
public:
        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);
                swap(temp);
                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);
    }
   
private:
        void print(ostream& s);
};

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

protected:
        object* object_;	//	the parent object

public:
        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);
                swap(temp);
                return *this;
        }

        void swap(control& c)
        {
                std::swap(sn, c.sn);
                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 (sn==c.sn);}
        bool operator<(const control& c) const {return (sn<c.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 
valarray<double>&ret_orientation);
        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.is, 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.m+dmass),
                        (u.is), (u.plan), (u.plan_start);
        }
        state res = state(dpos,dvel,dtime,dmass, (u.is), (u.plan), 
(u.plan_start));
        return res;
}

-- 
---
Harry Erwin, PhD, Senior Lecturer of Computing, University of 
Sunderland. Computational neuroscientist modeling bat bioacoustics 
and behavior. <http://www.cet.sunderland.ac.uk/~cs0her/index.html>