如何使用函数更新c ++类成员？ 错误样式 Verlet

问题描述

为了描述问题，我尝试使用代码中的对象简化三体问题。我有以下对象的代码：

Facades\App\Model\Structures::get()

现在，我相信速度Verlet算法是正确的，但是当我使用此主文件运行代码时：

#include <stdlib.h>
#include <cstdio>
#include <iostream>
#include <cmath>
#include <vector>
#include "star.h"

using namespace std;

Star::Star( double m,double x_p,double y_p,double x_v,double y_v )
{
    init( m,x_p,y_p,x_v,y_v);
}

void Star::init( double m,double y_v )
{
    Mass = m;
    X_Position = x_p;
    Y_Position = y_p;
    X_VeLocity = x_v;
    Y_VeLocity = y_v;
    R_Position[0] = X_Position;
    R_Position[1] = Y_Position;
    R_VeLocity[0] = X_VeLocity;
    R_VeLocity[1] = Y_VeLocity;
}

double Star::potential( Star star2,double dx,double dy )
{
    double G = 3.0548e34;
    double Potential;

    double x_component = X_Position - star2.X_Position + dx;
    double y_component = Y_Position - star2.Y_Position + dy;

    double R = sqrt(x_component*x_component + y_component*y_component);

    Potential = G* Mass* star2.Mass / R;
    return Potential;
}

double * Star::compute_forces( Star star2 )
{
    double h_x = ( X_Position - star2.X_Position )/1000;
    double h_y = ( Y_Position - star2.Y_Position )/1000;

    double *F = new double[2];

    F[0] = ( potential( star2,h_x,0.0 ) - potential( star2,-h_x,0.0 ) )/2*h_x;
    F[1] = ( potential( star2,0.0,h_y ) - potential( star2,-h_y ) )/2*h_y;

    return F;
}

void Star::verlet( Star star2,double h )
{
    double *Force = compute_forces( star2 );

    X_Position += h*X_VeLocity + 0.5*h*h*Force[ 0 ];
    Y_Position += h*Y_VeLocity + 0.5*h*h*Force[ 1 ];

    double *Force_new = compute_forces( star2 );

    X_VeLocity += 0.5*h*(Force[ 0 ] + Force_new[ 0 ] );
    Y_VeLocity += 0.5*h*(Force[ 1 ] + Force_new[ 1 ] );
}

Star对象的成员的值未更新：/。我有什么想念的吗？ 判罚的结果是这样的：

#include <iostream>
#include <fstream>
#include <cmath>
#include <cstdio>
#include "star.h"

using namespace std;

int main()
{
    Star star1( 50,0.0 );
    Star star2( 1.00,1.0,-1.0,1.0 );
    Star star3( 1.00,1.0 );

    Star arr[3] = { star1,star2,star3 };

    double h = 10/1000;

    //for ( double time = 0.0; time <= 10.0; )
    //{
        for ( int inst = 0 ; inst< 3; ++inst )
        {
            for ( int jnst = 0; jnst < 3; ++jnst )
            {
                if ( inst != jnst )
                {
                    arr[ inst ].verlet( arr[ jnst ],h );
                    double *pos = arr[ inst ].get_positions();
                    cout << " " << pos[ 0 ] << " " << pos[ 1 ] << endl;

                }
            }

                    }
        //time += h;
    //}

    return 0;
}

提前谢谢！

编辑：

我尝试为自己的部队实施 0 0 0 0 0 1 0 1 0 -1 0 -1 ，但最终出现了分段错误。

编辑2： 检查我的std::vector<double>方法后，我注意到它只返回初始化值。因此，我尝试实现了这一点：

get_positions()

它奏效了，所以我在我的主要代码中实现了以下内容。

std::vector<double> get_positions(){  std::vector<double> temp = { X_Position,Y_Position }; return temp; }

但是，现在我陷入了一个全新的问题...现在，我得到以下用于算法更新的数字！

        std::vector<double> p1 = star1.get_positions();
        std::vector<double> p2 = star2.get_positions();
        std::vector<double> p3 = star3.get_positions();

        cout << p1[ 0 ] << " " << p1[ 1 ] << " "  << p2[ 0 ] << " " << p2[ 1 ] << " " << p3[ 0 ] << " " << p3[ 1 ] << endl;

这意味着我在代码中某处乘以零。问题是我无法一生看到哪里。谢谢您的帮助！

解决方法

错误

如果要除以2*h_x，则需要将其写为/(2*h_x)，否则要除以2再乘以h_x，以得到极小的力值，因此不会移动系统。

作为补充，您在主程序中将步长定义为

double h = 10/1000;

右边的值被标识为整数除法的结果，0。有了这个步长，一切都不会改变。

样式

不要为相同的值构造两个数据字段，您必须确保这些字段始终同步。使用getter方法以其他格式显示数据。

对于科学而言，最好使用既定的向量类，然后再提供向量算术，例如boost / Eigen。

在构造函数中使用初始化列表语法，不需要init函数即可分配值。

Verlet

Verlet方法无法以这种方式工作。即使一切都在编码方面正确进行，结果还是一阶方法，既不保留能量也不动量。

• 有关在重力模拟中使用Verlet的信息，请参见N-Body Gravity Simulation in JavaScript。
• 在稍有不同的情况下，Lennard-Jones potential simulation相同。
• 也许有关Verlet方法的一般属性Velocity verlet algorithm not conserving energy的讨论可能会有所帮助。

简短的版本是，Verlet方法的阶段是外部框架。在每个阶段中，必须先对所有对象执行所有计算，然后才能进入下一个阶段。也就是说，所有速度发生变化，然后所有位置发生变化，然后计算并累积所有力，然后所有物体的速度随新的力/加速度发生变化。

混合这些步骤会破坏方法的顺序和所有保留属性。 （前两个阶段可以交错，因为对象之间没有交互。）

---

我使用Pleiades IVP测试套件示例的数据实施了一些建议的更改，因为所提供的数据导致系统迅速爆炸。

带有主Verlet循环的主程序solarsystem.c

#include <iostream>
#include <cstdio>
#include <vector>
#include "star.h"

using namespace std;

int main()
{
    vector<Star> arr = {
        Star( 1,3.0,0.0,0.0 ),Star( 2,-3.0,Star( 3,-1.0,2.0,Star( 4,-1.25 ),Star( 5,1.0 ),Star( 6,-2.0,-4.0,1.75,Star( 7,4.0,-1.5,0.0 )
    };
    int N = arr.size();
    double dt = 0.001;
    int count = 10;
    for ( double time = 0.0; time <= 3.0; time += dt)
    {
        for ( int inst = 0 ; inst< N; ++inst ) {
            arr[inst].Verlet_stage1(dt);
        }
        for ( int inst = 0 ; inst< N; ++inst ) {
            for ( int jnst = inst+1; jnst < N; ++jnst ) {
                arr[inst].acceleration(arr[jnst]);
            }
        }
        for ( int inst = 0 ; inst< N; ++inst ) {
            arr[inst].Verlet_stage2(dt);
        }
        if( 10 == count) {
            count = 0;
            for ( int inst = 0 ; inst< N; ++inst ) {
                cout << " " << arr[inst].Position[1] << " " << arr[inst].Position[0];
            }
            cout << "\n";
        }
        count++;
    }

    return 0;
}

以及带有标题的Star类的实现

#pragma once
#include <eigen3/Eigen/Dense>

typedef Eigen::Vector2d Vec2D;

const double G = 1;

class Star {
    public:
        Star( double m,double x_p,double y_p,double x_v,double y_v )
             :Mass(m),Position(x_p,y_p),Velocity(x_v,y_v) {};
        double Mass;
        Vec2D Position,Velocity,Acceleration;
        
        void Verlet_stage1(double dt);
        void Verlet_stage2(double dt);
        
        double potential(Star other);
        void acceleration(Star &other);
};

和语料

#include "star.h"

double Star::potential( Star other )
{

    Vec2D diff = Position-other.Position;
    double R = diff.norm();
    return G * Mass * other.Mass / R;
}

void Star::acceleration( Star &other )
{
    Vec2D diff = Position-other.Position;
    double R = diff.norm();
    Vec2D acc = (-G / (R*R*R)) * diff;
    Acceleration += other.Mass * acc;
    other.Acceleration -= Mass * acc;
}

void Star::Verlet_stage1( double dt )
{
    Velocity += (0.5*dt) * Acceleration;
    Position += dt*Velocity;
    Acceleration *= 0;
}
void Star::Verlet_stage2( double dt )
{
    Velocity += (0.5*dt) * Acceleration;
}

这将导致下面的轨迹。图像非常取决于步长dt，因为势函数的近奇点非常接近，也就是说，如果物体非常靠近，则能量和动量近守恒的辛方法的前景就会破裂。

,

我个人并不反对使用原始指针，但是如果对它们的管理不当，则会出现复杂的情况。我不知道这段代码是做什么的，更何况它是怎么做的！但是，我尝试改善了一些我可以观察到的错误，但是显然，此代码需要认真检查。我想这段代码的缺点仅仅是由于经验不足，可以理解。

https://gcc.godbolt.org/z/5zT5o9请记住，由于各种功能体中原始指针的使用（非管理），该代码仍在泄漏。