力学の入門 : 4 「ニュートンの運動方程式」

GRAM · 投稿記事 by **GRAM** » 13年前

さて、超初歩的な力学をやるのに最低限必要と思われる武器は大体そろったので、ニュートンの運動方程式（ふつう単なる運動方程式といった場合これを指す）に
ついて話したいと思います。

[ 1 ] 運動方程式とは

高校１年生で物理を習った場合はやったと思いますが、運動方程式とは次のような「超短い式」であらわされる方程式のことです
ma = F

この式が言っていることはごく単純です。
「物体のある時刻の加速度は、その時刻における力の向きと同じで、かつその大きさは、力の大きさを質量で割ったものに等しい」ということです。
これは直観的かと思います。

つまり物体は力を加えた方向に加速し、力が大きいほど早く加速します。　一方質量が大きいほどゆっくり加速し、軽いと加速は容易だということです。

このあまりにも当たり前な方程式は　「世の中の目に見えるほとんどの運動」を記述できてしまします。

[ 2 ]運動方程式の適用範囲

しかし運動方程式はどんなときにも使えるわけではありません。たとえば「加速する車に固定した座標系」などではそのまま使えません。
運動方程式が使えるのは

「慣性系」

と呼ばれる座標系に限ります。

慣性系とは「力を加えなければ物体は静止するか、等速度で運動するかのどちらか」という座標系です。

つまり自動車で加速するときなどは、後ろに引っ張られるような感覚とともにものが勝手に動いてしまいますので、運動方程式の適用範囲外といえます。

[ 3 ] プログラム

さて、加速度が力から求まってしまえば、運動をシミュレートすることは容易です。
つまり、加速度を力と質量から求め、　加速度から速度を、速度から位置を求めれば、もうあとは表示するだけです。
また、複数の力が加わるときは力を足していき最後に加速度を求めればいいので一旦力をメンバ変数にバッファしておいて、Updateで一気に加速度を求めるという
処理にしています。

このプログラムでは、マウスポインタの位置との距離に比例した力を質点に加え、マウスポインタの位置へ向かうようになっています。
それだけでは面白くないので、ついでに質点の速度に比例した抵抗を加えています。
抵抗を加えることでマウスポインタの位置まで来たら減速して止まる（後の減衰振動）ような運動をします。

ClearDrawScreen() をコメントアウトしてみるとわかりますが、質点は単純にマウスポインタの位置に一直線に向かうのではなく、
それまでの運動を継続しながらマウスポインタの位置に向かっているため、マウスの運動を追従しようとしているように見えます。

► スポイラーを表示

CODE:

#include 
#include 
#include 
#include 
using namespace Linear;

static const double dt = 1.0 / 60;
class Particle
{
public :
	//質点の作成
	static std::shared_ptr Create( double m, Vector s, Vector v )
	{
		return std::shared_ptr( new Particle(m,s,v) );
	}

	//質点に力を加える
	void AddForce( Vector f )
	{
		f_ += f;
	}

	//質点の位置を取得
	const Vector S() const
	{
		return s_;
	}

	//質点の速度を取得
	const Vector V() const
	{
		return v_;
	}

	//質点の更新
	void Update()
	{
		//運動方程式に従い加速度、速度、位置を更新
		a_ = f_ / m_;
		v_ += a_ * dt;
		s_ += v_ * dt;

		//力を0に
		f_ = Vector();
	}

	void Draw() const
	{
		DxLib::DrawCircle( s_.x, s_.y, 10, DxLib::GetColor( 0x2C, 0xCC, 0x7c ) );
	}

private :
	double			m_;	//質量		
	Vector	s_;	//位置
	Vector	v_;	//速度
	Vector	a_;	//加速度
	Vector	f_;	//力

private:
	Particle( double m, Vector s, Vector v )
		:m_(m),s_(s),v_(v)
	{
	}
	Particle( const Particle& );
	Particle& operator=( const Particle& );
};

//質点に毎フレームごとに力を与えるクラス
class Environment
{
public:
	static Environment& Instance()
	{
		static Environment e;
		return e;
	}

	//質点を設定
	void SetParticle( std::weak_ptr particle )
	{
		particle_ = particle;
	}

	//減衰を設定
	void SetH( double h )
	{
		h_ = h;
	}

	void Update()
	{
		if( auto p = particle_.lock() )
		{
			int x, y;
			//マウスの位置を取得
			DxLib::GetMousePoint( &x, &y );
			Vector mouse( x, y, 0);

			//マウスポインタと質点との相対位置ベクトルを取得
			Vector dist = mouse - p->S();

			//その距離に従って質点に力を加える
			p->AddForce( dist );

			//減衰　速度に比例した抵抗を加える
			p->AddForce( -(p->V()) * h_);
			
		}
	}

	void Draw()
	{
		if( auto p = particle_.lock() )
		{
			int x, y;
			//マウスの位置を取得
			DxLib::GetMousePoint( &x, &y );
			Vector mouse( x, y, 0);


			//マウスポインタと質点との相対位置ベクトルを取得
			Vector s =  p->S();
			Vector dist = mouse - s;

			DxLib::DrawLine( mouse.x, mouse.y, s.x, s.y, GetColor( 0xFF,0,0 ) );
			
		}
	}

private:
	std::weak_ptr particle_;
	double h_;							//減衰定数
private:
	Environment():h_(1){}
	Environment( const Environment& );
	Environment& operator=( const Environment& );
};


  
 
int WINAPI WinMain( HINSTANCE hInstance, HINSTANCE hPrevInstance, LPSTR lpCmdLine, int nCmdShow )
{
    ChangeWindowMode(true);
    if( DxLib_Init() ( 300, 200, 0 ), Vector() );

	Environment::Instance().SetParticle( particle );
	//数値を変えると、質点への抵抗が変化します。（大きいと液体の中を動いているような感じになる。0で抵抗なし）
	Environment::Instance().SetH( 0.75 );

    while( ProcessMessage() == 0 && CheckHitKey( KEY_INPUT_ESCAPE ) == 0 )
    {
        ClearDrawScreen() ;
        //メインループ
        {
			Environment::Instance().Update();
			particle->Update();

			Environment::Instance().Draw();
			particle->Draw();
        }
        ScreenFlip() ;
    }
 
    DxLib_End() ;
    return 0 ;
}

Vector.hの中身

► スポイラーを表示

CODE:

#pragma once
#include 
//#include 
#include 
#include 
//#include 

namespace Linear{

	

	static const double VTOL = 0.000001;

	template
	struct Vector{
	public:
		NType x;
		NType y;
		NType z;

		inline Vector():x(0),y(0),z(0){}

		template 
		inline Vector(OtherType init_x, OtherType init_y, OtherType init_z)
			:x(init_x),y(init_y),z(init_z){
		}
		
		
		template 
		inline Vector(const Vector& other)
			:x(other.x),y(other.y),z(other.z){
		}

		template 
		inline Vector(const OtherType& other)
			:x(other.x),y(other.y),z(other.z){
		}

		template 
		inline Vector& operator=(const Vector& other){
			x=other.x;
			y=other.y;
			z=other.z;
			return *this;
		}
		inline ~Vector(){}



		void Reset(){
			x = NType(0);
			y = NType(0);
			z = NType(0);
		}

		template 
		void Reset(OtherType init_x, OtherType init_y, OtherType init_z){
			x=init_x;
			y=init_y;
			z=init_z;
		}

		inline operator bool() const {return x||y||z;}
		inline bool operator==(const Vector& other) const{
			return abs(other.x-x)& other) const{
			return abs(other.x-x)>NType(VTOL) || abs(other.y-y)>NType(VTOL) || abs(other.z-z)>NType(VTOL);
		}

		//以下固有メンバ関数
		inline NType Length() const{ //ベクトルの大きさ
			return sqrt(x*x+y*y+z*z);
		}

		inline Vector& Normalize(){	//ベクトルの正規化
			NType m = sqrt( x*x + y*y + z*z );
			if(m 
		inline	Vector& operator+=(const Vector& v){
			x += v.x;
			y += v.y;
			z += v.z;
			return *this;
		}

		template 
		inline	Vector& operator-=(const Vector& v){
			x -= v.x;
			y -= v.y;
			z -= v.z;
			return *this;
		}

		template 
		inline Vector& operator*=(const OtherType s){
			x *= s;
			y *= s;
			z *= s;
			return *this;
		}

		template 
		inline Vector& operator/=(const OtherType s){
			x /= s;
			y /= s;
			z /= s;
			return *this;
		}


		inline Vector operator-() const{
			return Vector(-x, -y, -z);
		}


		NType& operator[](int i){
			assert( 1 
		struct RCMul{
			typedef typename std::remove_const::type type;
		};
		template 
		struct RCDev{
			typedef typename std::remove_const::type type;
		};
		
		template 
		struct RCMul, RHS>{
		};
 
		template 
		struct RCMul>{
		}; 

		/*template 
		struct RCMul, RHS>{
		};
 
		template 
		struct RCMul>{
		}; */
	}

	//-------------------------------------------------------------------------------------------------------------------------------------

	template
	Vector GetValue(const Vector& lhs){
		return Vector(lhs.x.GetValue(), lhs.y.GetValue(), lhs.z.GetValue());
	}



	//ベクトルの和、差
	template 
	inline	Vector operator+(const Vector& lhs, const Vector& rhs){
		return Vector(lhs.x+rhs.x, lhs.y+rhs.y, lhs.z+rhs.z);
	}

	template 
	inline	Vector operator-(const Vector& lhs, const Vector& rhs){
		return Vector(lhs.x-rhs.x, lhs.y-rhs.y, lhs.z-rhs.z);
	}
	
	 //内積
    template 
    decltype( NType1(0) * NType2(0) )
    operator* ( const Vector& lhs, const Vector& rhs ){

        return lhs.x*rhs.x + lhs.y*rhs.y + lhs.z*rhs.z;
    
	}


	//外積
	template 
	inline const Vector::type >
	operator^(const Vector& lhs, const Vector& rhs){
		return Vector::type >
			(
				lhs.y*rhs.z - lhs.z*rhs.y,
				lhs.z*rhs.x - lhs.x*rhs.z,
				lhs.x*rhs.y - lhs.y*rhs.x
			);
	}

	//スカラーとの乗除	
	
	template 
	const Vector::type >
	operator*(const Vector& v, const NType2& s){
		return Vector::type >
			(v.x*s, v.y*s, v.z*s);
	}


	template 
	const Vector::type >
	operator*(const NType1& s, const Vector& v){
		return Vector::type >
			(s*v.x, s*v.y, s*v.z);
	}


	template 
	const Vector::type >
	operator/(const Vector& v, NType2 s){
		return Vector::type >
			(v.x/s, v.y/s, v.z/s);
	}


	//ベクトルの三重積
	template 

	inline decltype(NType1(0) * NType2(0) * NType3(0) )
	TripleScalar(const Vector& u,const Vector& v,const Vector& w){
		return 
			(u.x * (v.y*w.z-v.z*w.y)) +
			(u.y * (v.z*w.x-v.x*w.z)) +
			(u.z * (v.x*w.y-v.y*w.x));
	}







	//マイナー関数
	//----------------------------------------------------------------------------------------------------------


	//ベクトルの方向成分を得る
	template 
	Vector GetComponent(const Vector& lhs, Vector comp){
			comp.Normalize();
			auto d = (comp / ( NType2(1)*NType2(1) )) * lhs;
			return comp*d;
	}


	template 
	Vector VAbs( const Vector& v){
		return Vector(abs(v.x), abs(v.y), abs(v.z) );
	}


	template 
	//法線ベクトル(3Dは外積でいいから作らない)
	inline Vector GetNormalVector2D(Vector& v){
		return Vector(-v.y,v.x,0.0);
	}


	//ベクトル間のコサイン
	template 
	inline decltype( NType(0) / NType(0) )
	GetCos(Vector& lhs, Vector& rhs){
		NType mag1 = lhs.Magnitude();
		NType mag2 = rhs.Magnitude();
		if(abs(mag1)
	inline const Vector VRotate2D(const Vector& v, NType2 sin, NType2 cos){
		return Vector(v.x*cos-v.y*sin, v.x*sin+v.y*cos, 0.0);
	}

	template 
	inline decltype( NType(0) * NType(0) )
	VDistanceSq(const Vector& v1, const Vector& v2){
		return (v1-v2)*(v1-v2);
	}
	
	template 
	inline T VEuclideanDistance(const Vector& v1, const Vector& v2){
		return (v1-v2).Length();
	}
	
	template 
	inline T VManhattanDistance(const Vector& v1, const Vector& v2){
		T x = v1.x-v2.x;
		T y = v1.y-v2.y;
		T z = v1.z-v2.z;
		if(x<T(0))x=-x;
		if(y<T(0))y=-y;
		if(z<T(0))y=-z;
		return x+y+z;
	}
}

実行結果

次回はばねの運動かな？

		6月 2025
日	月	火	水	木	金	土
1	2	3	4	5	6	7
8	9	10	11	12	13	14
15	16	17	18	19	20	21
22	23	24	25	26	27	28
29	30