如何使用函数更新 c++ class 成员？

Question

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

#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 x_p, double y_p, double x_v, 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, 0.0, -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 ] );
}

现在我相信velocity verlet算法是正确的，但是当我运行使用这个主文件的代码时：

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

using namespace std;

int main()
{
    Star star1( 50, 0.0, 0.0, 0.0, 0.0 );
    Star star2( 1.00, 0.0, 1.0, -1.0, 1.0 );
    Star star3( 1.00, 0.0, -1.0, 1.0, 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;
}

Star 对象的成员值未更新 :/。有什么我想念的吗？ cout 的输出是这样的：

提前致谢！

编辑：

我尝试为我的部队实施 std::vector<double>，但我最终遇到了分段错误。

编辑 2：检查我的 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;

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

5.66002e-320 2.31834e-316
1.132e-316 4.63669e-313
1.698e-319 6.95503e-316
1.132e-316 4.63669e-313
5.66002e-320 2.31834e-316
1.132e-316 4.63669e-313
1.698e-319 6.95503e-316
1.132e-316 4.63669e-313
5.66002e-320 2.31834e-316
1.132e-316 4.63669e-313
1.698e-319 6.95503e-316
1.132e-316 4.63669e-313

这意味着我在代码中某处乘以零。问题是我一辈子都看不到哪里。如有帮助，谢谢！

Answer 1

我个人并不反对使用原始指针，但如果管理不当，就会出现复杂情况。我不知道这段代码是做什么的，更不知道它是怎么做的！尽管如此，我已经尝试改善一些我可以观察到的错误，但显然这段代码需要认真修改。我想这段代码的不足之处只是经验不足，可以理解。

https://gcc.godbolt.org/z/5zT5o9请记住，由于在各种函数体中使用（non-manage）原始指针，此代码仍然存在泄漏。

Answer 2

错误

如果你想除以 2*h_x，你需要写成 /(2*h_x)，否则你除以 2 再乘以 h_x，给力的值很小，因此没有太多移动系统。

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

double h = 10/1000;

右边的值被识别为整数除法的结果，即0。有了这个步长，什么都不会改变。

风格

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

对于科学来说，最好使用已建立的向量 class，它还提供向量运算，例如 boost/Eigen.

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

维莱特

Verlet 方法不是这样工作的。即使一切顺利 coding-wise，结果也是既不保留能量也不保留动量的一阶方法。

有关将 Verlet 与重力模拟结合使用的信息，请参阅 N-Body Gravity Simulation in JavaScript。
在略有不同的上下文中相同，Lennard-Jones potential simulation。
也许这个关于 Verlet 方法的一般属性的讨论也会有所帮助。

简而言之，Verlet 方法的阶段是外部框架。在每个阶段，必须对所有对象进行所有计算，然后才能进入下一阶段。也就是说，所有速度发生变化，然后所有位置发生变化，然后所有力被计算和累加，然后所有速度随着所有对象的新forces/accelerations而变化。

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

我使用 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, 3.0, 0.0, 0.0 ),
        Star( 2, 3.0,-3.0, 0.0, 0.0 ),
        Star( 3,-1.0, 2.0, 0.0, 0.0 ),
        Star( 4,-3.0, 0.0, 0.0,-1.25 ),
        Star( 5, 2.0, 0.0, 0.0, 1.0 ),
        Star( 6,-2.0,-4.0, 1.75, 0.0 ),
        Star( 7, 2.0, 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 class with header

的实现

#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 作为势函数的近奇点，也就是说，如果物体靠得很近，那么近能量和动量守恒的辛方法的承诺就会破裂。

如何使用函数更新 c++ class 成员？

How can I update c++ class members with a function?

c++

oop

object

numerical-methods

错误

风格

维莱特