verlet积分的质点运动仿真算例-小球平抛
vertlet积分常用于质点牛顿动力学仿真,像分子动力学仿真,原理数值分析和博客都比较详细了,下面的MassPt这个类有
mS-位置,mV-速度,mF-力这几个牛顿运动方程中的质点的几个物理状态,测试时仿真的高中小球平抛运动5秒中,每0.1秒
一个仿真步,初始是在0高度,根据 h = 0.5*g*t^2, g取9.8,则可见s在t=5时mS.y=-122.5
y
|
* ---->x
*
*
*
下面算例在控制台下打印每个仿真步的位置,没有做多余的可视化。
#include <iostream> #include "vec3d.h" using namespace osg; class MassPt { public: MassPt() { mId = -1; mM = 10.0; mS.set(0, 0, 0); mV.set(0, 0, 0); mF.set(0, 0, 0); } MassPt(int id, double m, Vec3d s, Vec3d v, Vec3d f): mId(id), mM(m), mS(s), mV(v), mF(f) { } int mId; // id double mM; Vec3d mS; // position Vec3d mV; // velocity Vec3d mF; // force }; void verlet_step(MassPt& p, double dt) { MassPt t; p.mS = p.mS + p.mV * dt + p.mF * ((dt * dt) / (2.0*p.mM)); p.mV = p.mV + p.mF * (dt / p.mM); } void cout_state(MassPt& p, double t) { std::cout << "time:" << t << ", "; std::cout << "pos: (" << p.mS.x() << ", " << p.mS.y() << ", " << p.mS.z() << ")" << ";\n"; } void vertlet_test0() { MassPt a(0, 10.0, Vec3d(0, 0, 0), Vec3d(20, 0, 0), Vec3d(0, -98, 0)); double dt = 0.1; double t = 0.0; cout_state(a, t); for (size_t i = 0; i < 50; i++) { verlet_step(a, dt); t += dt; cout_state(a, t); } } int main() { vertlet_test0(); system("pause"); return 0; }
打印结果
time:0, pos: (0, 0, 0); time:0.1, pos: (2, -0.049, 0); time:0.2, pos: (4, -0.196, 0); time:0.3, pos: (6, -0.441, 0); time:0.4, pos: (8, -0.784, 0); time:0.5, pos: (10, -1.225, 0); time:0.6, pos: (12, -1.764, 0); time:0.7, pos: (14, -2.401, 0); time:0.8, pos: (16, -3.136, 0); time:0.9, pos: (18, -3.969, 0); time:1, pos: (20, -4.9, 0); time:1.1, pos: (22, -5.929, 0); time:1.2, pos: (24, -7.056, 0); time:1.3, pos: (26, -8.281, 0); time:1.4, pos: (28, -9.604, 0); time:1.5, pos: (30, -11.025, 0); time:1.6, pos: (32, -12.544, 0); time:1.7, pos: (34, -14.161, 0); time:1.8, pos: (36, -15.876, 0); time:1.9, pos: (38, -17.689, 0); time:2, pos: (40, -19.6, 0); time:2.1, pos: (42, -21.609, 0); time:2.2, pos: (44, -23.716, 0); time:2.3, pos: (46, -25.921, 0); time:2.4, pos: (48, -28.224, 0); time:2.5, pos: (50, -30.625, 0); time:2.6, pos: (52, -33.124, 0); time:2.7, pos: (54, -35.721, 0); time:2.8, pos: (56, -38.416, 0); time:2.9, pos: (58, -41.209, 0); time:3, pos: (60, -44.1, 0); time:3.1, pos: (62, -47.089, 0); time:3.2, pos: (64, -50.176, 0); time:3.3, pos: (66, -53.361, 0); time:3.4, pos: (68, -56.644, 0); time:3.5, pos: (70, -60.025, 0); time:3.6, pos: (72, -63.504, 0); time:3.7, pos: (74, -67.081, 0); time:3.8, pos: (76, -70.756, 0); time:3.9, pos: (78, -74.529, 0); time:4, pos: (80, -78.4, 0); time:4.1, pos: (82, -82.369, 0); time:4.2, pos: (84, -86.436, 0); time:4.3, pos: (86, -90.601, 0); time:4.4, pos: (88, -94.864, 0); time:4.5, pos: (90, -99.225, 0); time:4.6, pos: (92, -103.684, 0); time:4.7, pos: (94, -108.241, 0); time:4.8, pos: (96, -112.896, 0); time:4.9, pos: (98, -117.649, 0); time:5, pos: (100, -122.5, 0);
改变MassPt的初始的mV, 可以仿真抛体运动,像炮弹,这个mF一直没变的,如果mF也变化,或者是距离或者和
其它masspt有关,可以仿真多质点运动,像质点弹簧,分子间作用仿真,当然verlet精度有缺陷,有一些改进数值
方法。
下面是Vec3d.h,其中verlet_step中的3d向量运算表达向量运算的迭代公式很直接,是简单从osg的vec3d头文件拷贝过来的,注释了几个函数
/* -*-c++-*- OpenSceneGraph - Copyright (C) 1998-2006 Robert Osfield * * This library is open source and may be redistributed and/or modified under * the terms of the OpenSceneGraph Public License (OSGPL) version 0.0 or * (at your option) any later version. The full license is in LICENSE file * included with this distribution, and on the openscenegraph.org website. * * This library is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * OpenSceneGraph Public License for more details. */ #ifndef OSG_VEC3D #define OSG_VEC3D 1 //#include <osg/Vec2d> //#include <osg/Vec3f> #include <math.h> namespace osg { /** General purpose double triple for use as vertices, vectors and normals. * Provides general math operations from addition through to cross products. * No support yet added for double * Vec3d - is it necessary? * Need to define a non-member non-friend operator* etc. * Vec3d * double is okay */ class Vec3d { public: /** Data type of vector components.*/ typedef double value_type; /** Number of vector components. */ enum { num_components = 3 }; value_type _v[3]; /** Constructor that sets all components of the vector to zero */ Vec3d() { _v[0] = 0.0; _v[1] = 0.0; _v[2] = 0.0; } //inline Vec3d(const Vec3f& vec) { _v[0] = vec._v[0]; _v[1] = vec._v[1]; _v[2] = vec._v[2]; } //inline operator Vec3f() const { return Vec3f(static_cast<float>(_v[0]), static_cast<float>(_v[1]), static_cast<float>(_v[2])); } Vec3d(value_type x, value_type y, value_type z) { _v[0] = x; _v[1] = y; _v[2] = z; } //Vec3d(const Vec2d& v2, value_type zz) //{ // _v[0] = v2[0]; // _v[1] = v2[1]; // _v[2] = zz; //} inline bool operator == (const Vec3d& v) const { return _v[0] == v._v[0] && _v[1] == v._v[1] && _v[2] == v._v[2]; } inline bool operator != (const Vec3d& v) const { return _v[0] != v._v[0] || _v[1] != v._v[1] || _v[2] != v._v[2]; } inline bool operator < (const Vec3d& v) const { if (_v[0] < v._v[0]) return true; else if (_v[0] > v._v[0]) return false; else if (_v[1] < v._v[1]) return true; else if (_v[1] > v._v[1]) return false; else return (_v[2] < v._v[2]); } inline value_type* ptr() { return _v; } inline const value_type* ptr() const { return _v; } inline void set(value_type x, value_type y, value_type z) { _v[0] = x; _v[1] = y; _v[2] = z; } inline void set(const Vec3d& rhs) { _v[0] = rhs._v[0]; _v[1] = rhs._v[1]; _v[2] = rhs._v[2]; } inline value_type& operator [] (int i) { return _v[i]; } inline value_type operator [] (int i) const { return _v[i]; } inline value_type& x() { return _v[0]; } inline value_type& y() { return _v[1]; } inline value_type& z() { return _v[2]; } inline value_type x() const { return _v[0]; } inline value_type y() const { return _v[1]; } inline value_type z() const { return _v[2]; } /** Returns true if all components have values that are not NaN. */ //inline bool valid() const { return !isNaN(); } /** Returns true if at least one component has value NaN. */ //inline bool isNaN() const { return osg::isNaN(_v[0]) || osg::isNaN(_v[1]) || osg::isNaN(_v[2]); } /** Dot product. */ inline value_type operator * (const Vec3d& rhs) const { return _v[0] * rhs._v[0] + _v[1] * rhs._v[1] + _v[2] * rhs._v[2]; } /** Cross product. */ inline const Vec3d operator ^ (const Vec3d& rhs) const { return Vec3d(_v[1] * rhs._v[2] - _v[2] * rhs._v[1], _v[2] * rhs._v[0] - _v[0] * rhs._v[2], _v[0] * rhs._v[1] - _v[1] * rhs._v[0]); } /** Multiply by scalar. */ inline const Vec3d operator * (value_type rhs) const { return Vec3d(_v[0] * rhs, _v[1] * rhs, _v[2] * rhs); } /** Unary multiply by scalar. */ inline Vec3d& operator *= (value_type rhs) { _v[0] *= rhs; _v[1] *= rhs; _v[2] *= rhs; return *this; } /** Divide by scalar. */ inline const Vec3d operator / (value_type rhs) const { return Vec3d(_v[0] / rhs, _v[1] / rhs, _v[2] / rhs); } /** Unary divide by scalar. */ inline Vec3d& operator /= (value_type rhs) { _v[0] /= rhs; _v[1] /= rhs; _v[2] /= rhs; return *this; } /** Binary vector add. */ inline const Vec3d operator + (const Vec3d& rhs) const { return Vec3d(_v[0] + rhs._v[0], _v[1] + rhs._v[1], _v[2] + rhs._v[2]); } /** Unary vector add. Slightly more efficient because no temporary * intermediate object. */ inline Vec3d& operator += (const Vec3d& rhs) { _v[0] += rhs._v[0]; _v[1] += rhs._v[1]; _v[2] += rhs._v[2]; return *this; } /** Binary vector subtract. */ inline const Vec3d operator - (const Vec3d& rhs) const { return Vec3d(_v[0] - rhs._v[0], _v[1] - rhs._v[1], _v[2] - rhs._v[2]); } /** Unary vector subtract. */ inline Vec3d& operator -= (const Vec3d& rhs) { _v[0] -= rhs._v[0]; _v[1] -= rhs._v[1]; _v[2] -= rhs._v[2]; return *this; } /** Negation operator. Returns the negative of the Vec3d. */ inline const Vec3d operator - () const { return Vec3d(-_v[0], -_v[1], -_v[2]); } /** Length of the vector = sqrt( vec . vec ) */ inline value_type length() const { return sqrt(_v[0] * _v[0] + _v[1] * _v[1] + _v[2] * _v[2]); } /** Length squared of the vector = vec . vec */ inline value_type length2() const { return _v[0] * _v[0] + _v[1] * _v[1] + _v[2] * _v[2]; } /** Normalize the vector so that it has length unity. * Returns the previous length of the vector. * If the vector is zero length, it is left unchanged and zero is returned. */ inline value_type normalize() { value_type norm = Vec3d::length(); if (norm > 0.0) { value_type inv = 1.0 / norm; _v[0] *= inv; _v[1] *= inv; _v[2] *= inv; } return(norm); } }; // end of class Vec3d /** multiply by vector components. */ inline Vec3d componentMultiply(const Vec3d& lhs, const Vec3d& rhs) { return Vec3d(lhs[0] * rhs[0], lhs[1] * rhs[1], lhs[2] * rhs[2]); } /** divide rhs components by rhs vector components. */ inline Vec3d componentDivide(const Vec3d& lhs, const Vec3d& rhs) { return Vec3d(lhs[0] / rhs[0], lhs[1] / rhs[1], lhs[2] / rhs[2]); } } // end of namespace osg #endif