/* * odetest.cc * * Numerical ODE (ordinary differential equation) solver tester * with simple Qt & QTX-based GUI. * * Note: This is a quick test program. IF YOU INTEND TO SEE GOOD * PROGRAMMING, DO NOT LOOK AT THIS FILE. * * Copyright (c) 2004 by Wolfgang Wieser (wwieser@gmx.de) * * This file may be distributed and/or modified under the terms of the * GNU General Public License version 2 as published by the Free Software * Foundation. * * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. * */ #include #include #include #include #include #include #include #include #include #include #include #include //gcc -W -Wall -O2 odetest.cc -o odetest -fno-rtti -fno-exceptions -I. -IQt -fmessage-length=$COLUMNS -lqt-mt -L/opt/Qt/lib libqtxlib.a -I../anivision/src/ -I../hlib-build/include/ ../anivision-build/src/numerics/ode/libnumerics_ode.a ../anivision-build/src/numerics/libnumerics.a ../hlib-build/libhlib.a using namespace QTX; #include "../anivision/src/numerics/function.h" #include "../anivision/src/numerics/ode/odesolver.h" //------------------------------------------------------------------------------ const double m1=1.0,m2=1.0; const double l1=100.0,l2=100.0; const double g=9.81; inline double SQR(double x) { return(x*x); } class MyODEFunc : public NUM::Function_ODE { protected: // [overriding virtual] void function(double /*x*/,const double *y,double *result) { register double phi=y[0],phip=y[1],psi=y[2],psip=y[3]; double tmp = m2*l1*l2*sin(phi-psi); double a = -tmp*psip*psip - g*(m1+m2)*l1*sin(phi); double d = tmp*phip*phip - g* m2 *l2*sin(psi); double b = l1*l1*(m1+m2); double c = m2*l1*l2*cos(phi-psi); double e = l2*l2*m2; double det = (b*e-c*c); result[0]=phip; result[1]=(a*e-c*d)/det; result[2]=psip; result[3]=(b*d-c*a)/det; #if 0 // Control: calc energy: double T = 0.5*(m1+m2)*SQR(l1*phip) + m2*( 0.5*SQR(l2*psip) + l1*l2*phip*psip*cos(phi-psi) ); double U = -g*( (m1+m2)*l1*cos(phi) + m2*l2*cos(psi) ); double E = T+U; static int cnt=0; if(!(++cnt%100)) { fprintf(stderr,"\n%g",E); } #endif } public: MyODEFunc() : Function_ODE(4) {} ~MyODEFunc() {} }; //------------------------------------------------------------------------------ class MainWindow : public QWidget { private: XPainter *black,*blue,*red; NUM::ODESolver *odesolve[2]; NUM::Function_ODE *odefunc; double xn[2],h[2],yn[2][4]; int prev_xco[2],prev_yco[2]; bool run[2]; bool running; int iter; void _Redraw(); void timerEvent(QTimerEvent *); void keyPressEvent(QKeyEvent *); void mousePressEvent(QMouseEvent *); void paintEvent(QPaintEvent *); public: MainWindow(QWidget *parent=NULL,const char *name=NULL); ~MainWindow(); }; void MainWindow::_Redraw() { } void MainWindow::mousePressEvent(QMouseEvent * /*ev*/) { //int x0=ev->x(); //int x1=ev->y(); } void MainWindow::keyPressEvent(QKeyEvent *ev) { switch(ev->key()) { case Key_Return: case Key_Enter: break; case Key_Q: close(); break; case Key_P: if(running) killTimers(); else startTimer(10); running=!running; break; case Key_S: timerEvent(NULL); break; case Key_A: run[0]=!run[0]; break; case Key_B: run[1]=!run[1]; break; case Key_T: fprintf(stderr,"time=%g, %g\n",xn[0],xn[1]); break; } } void MainWindow::paintEvent(QPaintEvent *) { _Redraw(); } void MainWindow::timerEvent(QTimerEvent *) { for(int it=0; it<10; it++) { int xco[2],yco[2]; for(int i=0; i<2; i++) { if(!run[i]) { xco[i]=prev_xco[i]; yco[i]=prev_yco[i]; continue; } //odesolve[i]->SingleStep(*odefunc,xn[i],h[i],yn[i]); odesolve[i]->AdaptiveStep(*odefunc,xn[i],h[i],yn[i]); //printf("%g\n",h[i]); double phi=yn[i][0]; double psi=yn[i][2]; xco[i]=-int(l1*sin(phi)+l2*sin(psi) +0.5); yco[i]=-int(l1*cos(phi)+l2*cos(psi) +0.5); } if(xco[0]==xco[1] && yco[0]==yco[1]) { //black->DrawPoint(250+xco[0],250-yco[0]); if(iter) { black->DrawLine( 250+prev_xco[0],250-prev_yco[0], 250+xco[0],250-yco[0]); } } else { //blue->DrawPoint(250+xco[0],250-yco[0]); //red->DrawPoint(250+xco[1],250-yco[1]); if(iter) { blue->DrawLine( 250+prev_xco[0],250-prev_yco[0], 250+xco[0],250-yco[0]); red->DrawLine( 250+prev_xco[1],250-prev_yco[1], 250+xco[1],250-yco[1]); } } prev_xco[0]=xco[0]; prev_yco[0]=yco[0]; prev_xco[1]=xco[1]; prev_yco[1]=yco[1]; ++iter; //if(!(iter%100)) //{ fprintf(stderr,"iter=%d\n",iter); } //if(xn[0]>=200) //{ killTimers(); running=0; } #if 0 static int px=-1000,py=-1; double ss=log10(h[0]); int xx=int(3*xn[0]+0.5); int yy=int(520+40*ss+0.5); red->DrawLine(px,py,xx,yy); px=xx; py=yy; fprintf(stderr,"h=%g<<\n",h[0]); #endif } } MainWindow::MainWindow(QWidget *parent,const char *name) : QWidget(parent,name) { QPainter pnt(this); pnt.setBrush(QBrush(Qt::black,Qt::SolidPattern)); black=new XPainter(&pnt,true); pnt.setBrush(QBrush(Qt::blue,Qt::SolidPattern)); blue=new XPainter(&pnt,true); pnt.setBrush(QBrush(Qt::red,Qt::SolidPattern)); red=new XPainter(&pnt,true); assert(black && blue && red); blue->SetFunction(GXand); red->SetFunction(GXand); odesolve[0]=new NUM::ODESolver_RungeKutta(); odesolve[1]=new NUM::ODESolver_RungeKutta(); odefunc=new MyODEFunc(); assert(odesolve[0] && odesolve[1] && odefunc); // Start conditions: for(int i=0; i<2; i++) { xn[i]=0; h[i]=0.01; yn[i][0]=0+0.9*M_PI/1.0+i*1e-12; yn[i][1]=0; yn[i][2]=0+M_PI/2.0; yn[i][3]=0; odesolve[i]->SetAdaptivePars(i ? 1e-10 : 1e-10); run[i]=1; } resize(500,500); move(100,100); setBackgroundColor(Qt::white); setCaption("ODE test"); show(); startTimer(10); running=1; iter=0; } MainWindow::~MainWindow() { delete odesolve[0]; delete odesolve[1]; delete odefunc; delete black; delete blue; delete red; } char *prg_name="odetest"; int main(int argc,char **arg) { QApplication qapp(argc,arg); qapp.setStyle(new QCDEStyle(TRUE)); //qapp.setStyle(new QSGIStyle(TRUE)); MainWindow mainwin(NULL); qapp.setMainWidget(&mainwin); int r=qapp.exec(); return(r); }