小贝子编程

Runge-Kutta (RK4) 衍生品在C++



你好，我在C++中创建了一个小型的运动模拟。有一些物质点在击中球体时会移动和反弹，我想向学生展示欧拉、龙格-库塔和中点方法的差异。

但是当我切换到 Rungy-Kutta 模式时，我在某处犯了一个错误，所有材料点都消失了。我犯的错误在于solveRK4方法。 MinGW 开发者工作室项目的附件中，还有库文件夹可以放入 mingw 编译器目录：http://speedy.sh/CvDHj/LABO3.zip

当您按下"R"按钮时，当"E"到欧拉时，当"M"到中点时，它会切换到 RK4。

问题是我在solveRK4函数中哪里犯了错误？

结果（发布）如下所示：http://speedy.sh/h28VP/zad3.exe 按 R，点将从屏幕上消失。

朋克 - 点韦克托 - 矢量

主要：

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <GLglut.h>
#include <math.h>
#include <time.h>
#include "struktury.h"
#define YS 1.0f
#define PI 3.1415

double rot=0.0f;
int solver=0;
//--------------------------------------------------------------
//
//  to solve
//  {dr/dt=v
//  {d2r/dt2=F/m
//
//  Solvery:
//      solveEuler
//      solveMidpoint
//      solveRK4
//
//  
//--------------------------------------------------------------
void derivatives(Wektor *in, Wektor *out, Punkt *p){
    // dr/dt =v
    // d2r/dt2=F/m
    out[0]=in[1];
    out[1]=p->f*(1.0/p->m);
}

void calcForce(Punkt *p, Wektor v){
//  p->f = p->m * 9.72 -  (1 * v) ;
}
void solveEuler(Punkt  *p,  float dt){
    Wektor in[2] , out [2];
    in[0] = p->r ;
    in[1] = p->v ;
    derivatives(in, out ,p);
    p->v = p->v + out [1] * dt;
    p->r = p->r + out [0] * dt;
}
void solveMidPoint(Punkt *p, float dt){
    //r,v,f ->Wektor
    Wektor k1[2],k2[2];
    Wektor y_k[2];
    Wektor y_h[2];
    Wektor yout[2];
    y_h[0]=p->r;
    y_h[1]=p->v;
    derivatives(y_h,yout,p);
    k1[1]=yout[1];
    k1[0]=yout[0];
    y_k[0]=y_h[0]+k1[0]*0.5*dt;
    y_k[1]=y_h[1]+k1[1]*0.5*dt;
    derivatives(y_k,yout,p);
    k2[1]=yout[1];
    k2[0]=yout[0];
    p->r=p->r+k2[0]*dt;
    p->v=p->v+k2[1]*dt; 

        }
    //MISTAKE HERE  mistake here
    void solveRK4(Punkt *p, float dt){
    //r,v,f ->Wektor
    Wektor k1[2],k2[2],k3[2],k4[2];
    Wektor y_k[2];
    Wektor y_h[2];
    Wektor y_a[2];
    Wektor y_b[2];
    Wektor yout[2];
    y_h[0]=p->r;
    y_h[1]=p->v;
    derivatives(y_h,yout,p);
    k1[0]=yout[0];
    k1[1]=yout[1];
    y_k[0]=y_h[0]+k1[0]*0.5*dt;
    y_k[1]=y_h[1]+k1[1]*0.5*dt;
    derivatives(y_k,yout,p);
    k2[1]=yout[1];
    k2[0]=yout[0];
    y_a[0]=y_k[0]+k2[0]*0.5*dt;
    y_a[1]=y_k[1]+k2[1]*0.5*dt;
    derivatives(y_a,yout,p);
    k3[0]=yout[0];
    k3[1]=yout[1];  
    y_b[0]=y_a[0]+k3[0]*dt;
    y_b[1]=y_a[1]+k3[1]*dt;
    derivatives(y_b,yout,p);
    p->r=p->r+y_b[0]+yout[0]*dt;
    p->v=p->v+y_b[1]+yout[1]*dt;    
}

//--------------------------------------------------------------
// RESZTA KODU
//--------------------------------------------------------------


struct SferaN{
       Wektor r1;
       float r;
       float t;
       float color[3];
       SferaN *right;
};

Wektor g(0,-1.0,0);
 Punkt *root=NULL;
 SferaN *sroot=NULL;
 int loop=0;
SferaN *SphereAlloc(float R, Wektor r1, float t, float c[3])
{
       SferaN *tmp;
       if (!(tmp=(SferaN*)malloc(sizeof(SferaN))))
          return NULL;
        tmp->right=NULL;
        tmp->r=R;
        tmp->r1=r1;
        tmp->t=t;
    tmp->color[0]=c[0];
    tmp->color[1]=c[1];
    tmp->color[2]=c[2];
        return tmp;
}
void SphereTest(SferaN *s, Punkt *p)
{
 float d;
 Wektor n;
 float z;
 d=(s->r1-p->r).len();
 if (d-s->r<0)
 {
    n=(s->r1-p->r);
    n.norm();
    z=d-s->r;
    p->r=p->r+n*z;
    Wektor vs,vn;
    vn=n*(n*p->v);
    vs=p->v-vn;
    p->v=(vs-vn*s->t);
 }

}
void AddSphere(SferaN *ro, float R, Wektor r1, float t, float c[3])
{
     SferaN *tmp;
     for (tmp=ro; tmp->right!=NULL; tmp=tmp->right);
        tmp->right=SphereAlloc(R,r1,t,c);
}

Punkt *PointAlloc(float m, int flaga, Wektor r, Wektor v)
{
      Punkt *tmp;
      if (!(tmp=(Punkt*)malloc(sizeof(Punkt))))
         return NULL;
      tmp->m=m;
      tmp->flag=flaga;
      tmp->r=r;
      tmp->v=v;
      tmp->right=NULL;
      return tmp;
}
void AddPoint(Punkt *ro, float m, int flag, Wektor r, Wektor v)
{
     Punkt *tmp;
     for (tmp=ro; tmp->right!=NULL; tmp=tmp->right);
        tmp->right=PointAlloc(m,flag,r,v);
}
Wektor W_Euler(Wektor f, float h)
{
       return (f*h);
}
void calcDeriv(Punkt *p){
}
void Solver(Punkt *ro, float dt)
{
    //0 euler, 1 midpoint, 2 RK4
     Punkt *tmp;
     for (tmp=ro;tmp!=NULL;tmp=tmp->right)
     {
         switch (solver){
         case 0:
             solveEuler(tmp,dt);
             break;
         case 1:
             solveMidPoint(tmp,dt);
             break;
         case 2:
             solveRK4(tmp,dt);
             break;
        default:
            solveEuler(tmp,dt);
         }
    /*
     tmp->dv=W_Euler(tmp->f*(1/tmp->m),dt);
     tmp->v=tmp->v+tmp->dv;
     tmp->dr=tmp->v*dt;
     tmp->r=tmp->r+tmp->dr;
    */
        /*
    Wektor a=tmp->f*(1.0/tmp->m);        
         //--------------
    Wektor k1=W_Euler(a,dt);
    Wektor k2=W_Euler(a+0.5*k1,dt);
    Wektor k3=W_Euler(a+0.5*k2,dt);
    Wektor k4=W_Euler(a+k3,dt);
    Wektor kk=(1.0/6.0)*(k1+2.0*k2+2.0*k3+k4);
    tmp->v=tmp->v+kk;
    tmp->r=tmp->r+tmp->v*dt;
    */
     }
}
void Sily(Punkt *ro)
{
     Punkt *tmp;
     for (tmp=ro;tmp!=NULL; tmp=tmp->right)
     {
        tmp->f=g*tmp->m;
     }
}


void AnimateScene(void)
{
         glClearColor(0.0f,0.0f,0.0f,0.0f);
     //if (loop==0)
        glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
        Punkt *tmp;
        //Sily(root);
        Solver(root,0.008); 
        Sily(root);
    //  glRotatef(0.1,0.0f,1.0f,0.0f);
rot+=0.01;
//if (rot>36) rot=0.0f;
        glPointSize(5);
  glDisable(GL_LIGHTING); 
  glDisable(GL_LIGHT0); 
        for(tmp=root;tmp!=NULL;tmp=tmp->right)
        {
        //glBegin(GL_POINTS);
               //glColor3f(rand()/(float)RAND_MAX,rand()/(float)RAND_MAX,rand()/(float)RAND_MAX);
               glColor3f(1,1,1);
              // glVertex3f(tmp->r.x/1.0,tmp->r.y/1.0,tmp->r.z/1.0);
            glPushMatrix();
            glTranslatef(tmp->r.x/1.0,tmp->r.y/1.0,tmp->r.z/1.0);
            glutSolidSphere(0.011,5,5);
            glPopMatrix();
               SferaN *stmp;
               for(stmp=sroot;stmp!=NULL;stmp=stmp->right)
               {
                  SphereTest(stmp,tmp);
               }
               //glVertex2f(0,0);
               //printf("%f %fn",tmp->r.x, tmp->r.y);
       // glEnd();
        glBegin(GL_LINES);
                glColor3f(1,0,0);
                double x1,x2,y1,y2,z1,z2;
                x1=tmp->r.x;
                y1=tmp->r.y;
                x2=tmp->v.x;
                y2=tmp->v.y;
        z1=tmp->r.z;
        z2=tmp->v.z;

                glVertex3f(x1,y1,z1);
                glVertex3f(x1+(x2)*0.05,y1+(y2)*0.05,z1+(z2)*0.05);
        glEnd();
        }
  glEnable(GL_LIGHTING); 
  glEnable(GL_LIGHT0);         
        SferaN *stmp;
        for(stmp=sroot;stmp!=NULL;stmp=stmp->right)
        {
         glPushMatrix();
         glTranslatef(stmp->r1.x,stmp->r1.y,stmp->r1.z);
    glColor3fv(stmp->color);
         glutSolidSphere(stmp->r,50,50);
         glPopMatrix();
        }

      //  printf("****n");
        glFlush();
        glutSwapBuffers();
        loop++;

        if (loop>2000)
        {
           double vx,vy,vz;
 double zz=0.01;
           for(tmp=root;tmp!=NULL;tmp=tmp->right)
           {
        vx=0.5-rand()/(float)RAND_MAX;
        vy=1.0-rand()/(float)RAND_MAX;
        vz=0.5-rand()/(float)RAND_MAX;
        vz=0.0f;
        vy*=YS;
            vy+=zz;
            zz+=0.001;
            tmp->r.x=0;
            tmp->r.y=0;
            tmp->v=Wektor(vx,vy,vz);
           }
            loop=0;
        }
}
void InitGraphics()
{
  GLfloat mat_specular[] = { 1.0, 1.0, 1.0, 1.0 }; 
  GLfloat mat_shininess[] = { 90.0 }; 
  GLfloat light_position[] = {1.0, 1.0, 1.0, 0.0}; 
  glMaterialfv(GL_FRONT, GL_SPECULAR, mat_specular); 
  glMaterialfv(GL_FRONT, GL_SHININESS, mat_shininess); 
  glLightfv(GL_LIGHT0, GL_POSITION, light_position); 
  glEnable(GL_LIGHTING); 
  glEnable(GL_LIGHT0); 
  glDepthFunc(GL_LEQUAL); 
  glEnable(GL_DEPTH_TEST); 

        glShadeModel(GL_SMOOTH);


    glEnable(GL_COLOR_MATERIAL);

}

void myReshape(GLsizei w, GLsizei h) 
{ 
  glViewport(0, 0, w, h); 
  glMatrixMode(GL_PROJECTION); 
  glLoadIdentity(); 
    double p1=1.0f;
  if(w <= h) 
     glOrtho(-p1,p1,-p1*(GLfloat)h/(GLfloat)w,p1*(GLfloat)h/(GLfloat)w,-10.0,10.0); 
  else 
     glOrtho(-p1*(GLfloat)w/(GLfloat)h,p1*(GLfloat)w/(GLfloat)h,-p1,p1,-10.0,10.0); 
  glMatrixMode(GL_MODELVIEW); 
  glLoadIdentity(); 
} 
void keyboard(unsigned char key, int x, int y){
    switch (key){
    case 'e':
        solver=0;
        printf("Solver --> Eulern");
        break;
    case 'm':
        solver=1;
        printf("Solver --> MidPointn");
        break;
    case 'r':
        solver=2;
        printf("Solver --> RK4n");
        break;
    }
}
void idle(void)
{
        glutPostRedisplay();
}
int main(int argc, char *argv[])
{
 double vx,vy,vz;
 int i;
 srand(time(NULL));
 double zz=0.01;
 for (i=1; i<620; i++)
 {
        vx=0.5-rand()/(float)RAND_MAX;
        vy=1.0-rand()/(float)RAND_MAX;
        vz=0.5-rand()/(float)RAND_MAX;
    vz=0.0f;
        vy*=YS;
        vy+=zz;
        zz+=0.001;
    if (!root)
       root=PointAlloc(1,0,Wektor(0,0,0),Wektor(vx,vy,vz));
    else
        AddPoint(root,1,0,Wektor(0,0,0),Wektor(vx,vy,vz));
 }
 zz=-0.2;
 for (i=1; i<140; i++)
 {
       AddPoint(root,10,0,Wektor(zz,1,0),Wektor(0,0,0));
       zz+=0.01;
 }
    float c1[] = {0,0,1};
    float c2[] = {0,1,0};
    float c3[] = {1,1,0};
    float c4[] = {1,0,1};
    float c5[] = {0.6,0.5,1};
    float c6[] = {0.6,0.5,0.3};
  sroot=SphereAlloc(0.5,Wektor(0,-0.5,0),0.9,c1);
  AddSphere(sroot,0.1,Wektor(-0.5,0.5,0),0.8,c2);
  AddSphere(sroot,0.1,Wektor(0.5,0.5,0),0.3,c3);
  AddSphere(sroot,0.1,Wektor(1,0.0),0.3,c4);
  AddSphere(sroot,0.15,Wektor(0,0.7,0),0.9,c5);
  AddSphere(sroot,0.2,Wektor(-1,-0.2,0),1.0,c6);
        Sily(root);
        glutInit(&argc, argv);
        glutInitWindowSize(600,600);
        glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE | GLUT_DEPTH);
        glutCreateWindow("GLUT example");

        InitGraphics();
        glutDisplayFunc(AnimateScene);
        glutIdleFunc(idle);
    glutKeyboardFunc(keyboard);
    glutReshapeFunc(myReshape); 
        glutMainLoop();
return 0;
}

结构：

#ifndef __STRUKTURY_H
#define __STRUKTURY_H

class Wektor{
public:
       Wektor(double _x=0, double _y=0, double _z=0):
       x(_x),y(_y),z(_z){}
       Wektor operator+(Wektor const &);
       Wektor operator-(Wektor const &);
       double operator*(Wektor const &);
       Wektor operator*(double);
       double len();
       void norm();
       double x,y,z;
};
Wektor operator*(double, Wektor const &);

Wektor Wektor::operator-(Wektor const &w)
{
       return Wektor(x-w.x,y-w.y,z-w.z);
}
double Wektor::operator*(Wektor const &w)
{
       return (x*w.x + y*w.y + z*w.z);
}
double Wektor::len()
{
       return sqrt((*this)*(*this));
}
Wektor Wektor::operator+(Wektor const &w)
{
       return Wektor(x+w.x,y+w.y,z+w.z);
}
Wektor Wektor::operator*(double l)
{
       return Wektor(x*l, y*l, z*l);
}

Wektor operator*(double s, Wektor const &w)
{
       return Wektor(s*w.x, s*w.y, s*w.z);
}
void Wektor::norm()
{
     double d = len();
     if (d)
        (*this)=(*this)*(1/d);
}
struct Punkt{
        int flag;
        float m;
        Wektor f;
        Wektor r;
        Wektor v;
        Wektor dr;
        Wektor dv;
        Punkt *right;
};
#endif

这些行中的某些索引似乎是错误的：

y_b[1]=y_a[0]+k3[0]*dt;
p->v=p->v+y_b[0]+yout[1]*dt;

请仔细阅读您的 RK4 并检查所有索引。甚至可能有更多的复制/粘贴错误。

我已经修复了它，现在工作正常！ ：-）

0.在开始之前 - 非常仔细地阅读 WIKI 段落！

1.你在实际集成部分做错了，你写了

 p->r = p->r + y_b[0] + yout[0] * dt;

这与维基百科版本的经典 4 阶龙格-库塔 Mtehod 无关。 这是一条直肠线：

 p->r = p->r + (yout[0][0] + 2* yout[1][0] + 2* yout[2][0] + yout[3][0] )*  dt * (1.0/6);

所以，再看看维基并进行比较。

2.后来我发现了一系列的台词 

y_a[...] = y_k[...] + k2[...] * 0.5*dt;
y_b[...] = y_a[...] + k3[...] * dt;

你应该写： 

y_a[0] = y_h[0] + k2[0] * 0.5*dt;
y_b[0] = y_h[0] + k2[0] * 0.5*dt;
使用

INPUT 但不使用新计算的代理（临时）值。你惹了那个y_h，y_k，y_a y_b...y_qwertyuiop!!!!不要这样命名

）））

3.你也有太多的变量...我使它更加简单易读，也使用更少的内存。

反正少字多处理：这是工作方法代码，复制粘贴享受：

void solveRK4(Punkt *p, float dt) {
    Wektor y_in[2];
    Wektor y_temp[2];
    Wektor k[4][2];    
    y_in[0] = p->r;
    y_in[1] = p->v;
    derivatives(y_in, k[0], p);
    y_temp[0] = y_in[0] + k[0][0] * 0.5*dt;
    y_temp[1] = y_in[1] + k[0][1] * 0.5*dt;
    derivatives(y_temp, k[1], p);
    y_temp[0] = y_in[0] + k[1][0] * 0.5*dt;
    y_temp[1] = y_in[1] + k[1][1] * 0.5*dt;
    derivatives(y_temp, k[2], p);
    y_temp[0] = y_in[0] + k[2][0] * dt;
    y_temp[1] = y_in[1] + k[2][1] * dt;
    derivatives(y_temp, k[3], p);
    p->r = p->r + (k[0][0] + 2* k[1][0] + 2* k[2][0] + k[3][0] )*  dt * (1.0/6);
    p->v = p->v + (k[0][1] + 2 * k[1][1] + 2 * k[2][1] + k[3][1]) * dt * (1.0/6);
}

4.不要使用GLUT;使用GLFW！ ;-）

相关内容

  • 没有找到相关文章

最新更新


热门标签：

javascript python java c# php android html jquery c++ css ios sql mysql arrays asp.net json python-3.x ruby-on-rails .net sql-server django objective-c excel regex ruby linux ajax iphone xml vba spring asp.net-mvc database wordpress string postgresql wpf windows xcode bash git oracle list vb.net multithreading eclipse algorithm macos powershell visual-studio image forms numpy scala function api selenium