太阳系模拟项目(速度Verlet帮助) [英] Solar System Simulation Project (velocity verlet help)

问题描述

对于我的建模和模拟课程，我想模拟一个太阳系.我只是从一颗恒星(太阳)和一颗行星(地球)开始，但是我已经遇到了一些问题.我现在花了一些时间来回顾和学习不同的公式和方法，以模拟恒星和周围物体如何影响行星的轨道.我想使用速度Verlet，并最终研究n体问题.我的速度Verlet函数遇到很多问题.有时，它的作用就好像它在正常地绕地球运行一样，然后扭曲"将地球移到某个随机位置.我还注意到我从来没有得到过负"加速度，所以我的x和y坐标一致.总是在增加，所以我看不出地球应该如何环绕太阳.任何帮助是极大的赞赏.我只有AGK :: Prints，所以我可以看到不同的变量如何变化.

For my modelling and simulation class project, I want to simulate a solar system. I'm starting with just a star (sun) and a planet (earth), but I'm running into a few problems already. I've spent some time now just reviewing and learning about different formulas and way to simulate how the planet's orbits will be affected by the star and surrounding objects. I want to use velocity verlet and eventually look into the n-body problem. I'm having numerous issues with my velocity verlet function. Sometimes it acts as if it's making earth orbit normally and then it " warp drives" earth off into some random location. I've also noticed I never get a "negative" acceleration, so my x and y coord. are always increasing, so I don't see how the earth is suppose to wrap back around the sun. Any help is greatly appreciated. The AGK::Prints I have just so I can see how the different variables are changing.

double velocityVerlet(float positionCalc, double position2, 
                      float &velocity, double massCalc, double mass2)
//positionCalc is the position being updated, position 2 is position of 
// other object, same with mass
{
    float force = forceFunc(positionCalc, position2, massCalc, mass2);
    agk::PrintC("Force is: ");
    agk::Print(force);
    float acceleration = accelerationFunc(force,massCalc);
    agk::PrintC("Accel is: ");
    agk::Print(acceleration);`;

    double newAccel = 0;

    positionCalc = positionCalc + velocity*dt + 
                   (.5*acceleration)*pow(dt,2); //calculates new position
    agk::PrintC("New Position is: ");
    agk::Print(positionCalc);
    force = forceFunc(positionCalc,position2,massCalc,mass2);
    newAccel = accelerationFunc(force, massCalc);

    velocity = velocity + .5*(acceleration + newAccel)*dt; //new velocity
    agk::PrintC("Velocity is: ");
    agk::Print(velocity);

    return positionCalc;
}

推荐答案

您的积分器接受标量并且您的问题与二维系统有关，这一事实使我认为您要两次调用积分器，每个组件一次.这根本行不通，因为您的系统将在相空间中进行不切实际的移动.积分器可处理矢量:

The facts that your integrator accepts scalars and that your question is about 2-dimensional system makes me think that you are calling the integrator twice, once for each component. This simply won't work since your system will be taking unrealistic moves through the phase space. The integrator works with vector quantities:

X (t + dt)= X (t)+ V (t)dt +(1/2) A (t)dt ²

X(t+dt) = X(t) + V(t) dt + (1/2) A(t) dt²

V (t + dt)= V (t)+(1/2)( A (t)+ A (t + dt))dt

V(t+dt) = V(t) + (1/2)(A(t) + A(t+dt)) dt

此处 X (t)是由所有粒子的坐标组成的列向量-这是系统相空间的配置子空间. V (t)是所有粒子速度的列向量，从技术上讲代表了动量子空间. A (t)也是如此.那些必须同时进行更新，而不是分别进行更新.

Here X(t) is a column-vector that consists of the coordinates of all particles - this is the configuration subspace of the system's phase space. V(t) is a column-vector of the speeds of all particles, technically representing the momentum subspace. The same applies to A(t). Those have to updated simultaneously, not separately.

对于不依赖于速度(例如经典重力)的力场，整个速度Verlet过程在代码中的解释如下:

The whole velocity Verlet procedure translates as follows in code for force fields that do not depend on the speed (e.g. classical gravity):

Vector forces[num_particles];

// Compute initial forces
forces = computeForces(positions);

for (int ts = 0; ts < num_timesteps; ts++)
{
   // Update positions and half-update velocities
   for (int i = 0; i < num_particles; i++)
   {
      positions[i] += velocities[i]*dt + 0.5*(forces[i] / m[i]) * dt*dt;
      velocities[i] += 0.5*(forces[i] / m[i]) * dt;
   }

   // Compute new forces and half-update velocities
   forces = computeForces(positions);

   for (int i = 0; i < num_particles; i++)
   {
      velocities[i] += 0.5*(forces[i] / m[i]) * dt;
   }
}

请注意，在进行下一轮力评估之前，将首先更新所有位置.另外，由于位置在速度的第二次更新期间不会改变，因此每次迭代仅需要评估一次力.在上面的示例代码中，Vector是一个实现n维向量并包含n分量的类(例如，在2d情况下为2).它还重载了+和+=运算符以实现矢量(逐分量)加法，以及使*和/运算符实现了标量的乘法/除法.这只是为了说明这种情况，可以用每个位置/速度矢量的分量上的内部循环代替.

Note that all positions are updated first before the next round of force evaluation. Also it is only necessary to evaluate the forces once per iteration since the positions do not change during the second update to the velocities. In the example code above Vector is a class that implements an n-dimensional vector and holds n components (e.g. 2 in your 2d-case). It also overloads the + and += operators to implement vector (component-wise) addition, as well as * and / to implement multiplication / division by a scalar. This is just to illustrate the case and can be replaced by inner loops over the components of each position/velocity vector.

正确选择时间步非常重要.时间步太短会大大降低仿真速度.太大的时间步长将导致无法完成的物理操作，例如跳跃的行星.

The correct choice of time step is very important. Too small time step will slow down the simulation significantly. Too large time step will result in impossible physics, e.g. jumping planets.

这篇关于太阳系模拟项目(速度Verlet帮助)的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！