如何使用函数更新c ++类成员? [英] How can I update c++ class members with a function?
问题描述
为了描述问题，我尝试使用代码中的对象简化三体问题.我为该对象提供了以下代码:
#include< stdlib.h>#include< cstdio>#include< iostream>#include< cmath>#include< vector>#include"star.h"使用命名空间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);}无效Star :: init(double m，double x_p，double y_p，double x_v，double y_v){质量= 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){双倍G = 3.0548e34;双电势；double x_component = X_Position-star2.X_Position + dx;双重y_component = Y_Position-star2.Y_Position + dy;双R = sqrt(x_component * x_component + y_component * y_component);电位= G *质量*星形2.质量/R;回报潜力;}double * Star :: compute_forces(Star star2){双h_x =(X_Position-star2.X_Position)/1000;double h_y =(Y_Position-star2.Y_Position)/1000;double * F =新的double [2];F [0] =(电势(star2，h_x，0.0)-电势(star2，-h_x，0.0))/2 * h_x;F [1] =(电势(star2，0.0，h_y)-电势(star2，0.0，-h_y))/2 * h_y;返回F;}无效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];双* 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]);} 现在，我相信速度Verlet算法是正确的，但是当我使用此主文件运行代码时:
#include< iostream>#include< fstream>#include< cmath>#include< cstdio>#include"star.h"使用命名空间std;int main(){星星star1(50，0.0，0.0，0.0，0.0);星星star2(1.00，0.0，1.0，-1.0，1.0);星星star3(1.00，0.0，-1.0，1.0，1.0);Star arr [3] = {star1，star2，star3};双h = 10/1000;//for(double time = 0.0; time< = 10.0;)//{for(int inst = 0; inst< 3; ++ inst){对于(int jnst = 0; jnst< 3; ++ jnst){如果(inst！= jnst){arr [inst] .verlet(arr [jnst]，h);双* pos = arr [inst] .get_positions();cout<<""<<pos [0]<<""<<pos [1]<<恩德尔}}}//时间+ = h;//}返回0;} Star对象的成员的值未更新:/.我有什么想念的吗?判罚的结果是这样的:
0 00 00 10 10 -10 -1 提前谢谢！
我尝试为自己的部队实施 std :: vector< double> ，但最终出现了分段错误.
在检查我的 get_positions()方法后，我注意到它仅返回初始化的值.因此，我尝试实现了这一点:
std :: vector< double>get_positions(){std :: vector< double>temp = {X_Position，Y_Position};返回温度} 它起作用了，所以我在我的主要代码中实现了以下内容.
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]<<恩德尔 但是，现在我陷入了一个全新的问题...现在，我得到以下用于算法更新的数字！
<代码> 5.66002e-320 2.31834e-3161.132e-316 4.63669e-3131.698e-319 6.95503e-3161.132e-316 4.63669e-3135.66002e-320 2.31834e-3161.132e-316 4.63669e-3131.698e-319 6.95503e-3161.132e-316 4.63669e-3135.66002e-320 2.31834e-3161.132e-316 4.63669e-3131.698e-319 6.95503e-3161.132e-316 4.63669e-313
这意味着我在代码中某处乘以零.问题是我无法一生看到哪里.谢谢您的帮助！
错误
如果要除以 2 * h_x ，则需要将其写为/(2 * h_x)，否则，除以2并乘以 h_x ，给出的力值很小，因此不会使系统移动太多.
作为补充，您在主程序中将步长定义为
double h = 10/1000; 右侧的值被标识为整数除法的结果，该整数除法是 0 .有了这个步长，一切都不会改变.
样式
不要为相同的值构造两个数据字段，您必须确保这些字段始终同步.使用getter方法以其他格式显示数据.
对于科学而言，最好使用既定的向量类，然后再提供向量算术，例如boost/Eigen之一.
在构造函数中使用初始化列表语法，不需要 init 函数仅分配值.
Verlet
Verlet方法不能以这种方式工作.即使一切都在编码方面正确进行，结果仍是一种既不保留能量也不动量的一阶方法.
- 有关将Verlet用于重力模拟的信息，请参见
To describe the problem, I am trying to use objects in my code to stream line a three body problem. I have the following code for the object:
#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 ] ); }Now I believe that the velocity verlet algorithm is correct, but when I run the code using this main file:
#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; }The values of members of the Star object are not updating :/. Is there something I am missing? the out put of the cout is this:
0 0 0 0 0 1 0 1 0 -1 0 -1Thank you in advance!
Edit:
I tried implementing a
std::vector<double>for my forces, but I ended up with a segmentation fault.Edit 2: After checking my
get_positions()method I noticed it was returning only the initialized values. So I tried implementing this:std::vector<double> get_positions(){ std::vector<double> temp = { X_Position , Y_Position }; return temp; }And it worked so i implemented the following into my main code.
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;However now I am stuck on a completely new problem... Now I am getting the following numbers for the algorithm updates!
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-313Which means some where I am multiplying by zeros somewhere in my code. The problem is I cant for the life of me see where. Thanks if there is any help!
解决方案Error
If you want to divide by
2*h_x, you need to write this as/(2*h_x), else you divide by 2 and multiply byh_x, giving miniscule values for forces and thus not moving the system by much.To complement this, you defined the step size in the main program as
double h = 10/1000;The value on the right is identified as result of an integer division, which is
0. With this step size nothing will change.Style
Do not construct two data fields for the same value, you would have to ensure that these fields are always synchronized. Use getter methods to present data in a different format.
For science it would be better to use an established vector class that then also provides vector arithmetic, like the one of boost/Eigen.
Use initialization list syntax in the constructor, you do not need an
initfunction to just assign the values.Verlet
The Verlet method does not work this way. Even if everything goes right coding-wise, the result is a first order method that neither preserves energy nor momentum.
- For using Verlet with gravity simulations, see N-Body Gravity Simulation in JavaScript.
- The same in a slightly different context, Lennard-Jones potential simulation.
- Perhaps also this discussion about general properties of the Verlet method, Velocity verlet algorithm not conserving energy, can be helpful.
The short version is, the stages of the Verlet method are the outer frame. In each stage, all computations have to be carried out for all objects before changing to the next stage. That is, all velocities change, then all positions chance, then all forces are computed and accumulated, then all velocities change with the new forces/accelerations for all objects.
Mixing these steps destroys the order of the method and all conservation properties. (The first two stages can be interleaved, as there is no interaction between objects.)
I implemented some of the suggested changes, using the data of the Pleiades IVP test suite example, as the provided data lead to a rapid explosion of the system.
The main program
solarsystem.cwith the main Verlet loop#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; }and the implementation of the
Starclass 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); };and corpus
#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; }This results in the trajectories below. The picture is very depending on the step size
dtas near singularities of the potential function, that is, if bodies come very close together, the promise of symplectic methods of near conservation of energy and momentums breaks apart.这篇关于如何使用函数更新c ++类成员?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
