2010夏合宿藤原 - P.S.P.T. 計算物理研究会

減衰のない強制振動

• 28期　藤原
• 2010年度夏合宿
• 使用言語 C言語
• 対象OS
• ライブラリ等

drawps2

• 添付資料
• 添付ソースコード

概要 結果 考察 感想 参考文献 使用方法 ソースコード

---

概要

抵抗力のない強制振動の運動方程式は,
質量 の質点とバネ定数 のばねからなる系に,周期的な外力 が加わっているとき,

となる。
両辺を で割り、 を 、 を とおくと,

となる。

のとき

外力はゆっくり変化し、強制振動は外力と同位相になる。

のとき

外力は早く変化し、強制振動は だけ位相が遅れて外力に対し逆位相で振動する。

のとき

共鳴する。

のとき

うなりが観測できる。

結果

考察

感想

参考文献

---

使用方法

ソースコード

#include "drawps2.h"
 
#define GN 160
#define dt 0.01
 
DRAWPS2 drps2, graph, phase;
 
// 変数     
double t,x,v;
double f = 1.7;
double omega = 1.2;
double m = 4.0;
double k = 4.0;
double omega0 = sqrt(k/m);
double T, U, E;
double graphX[GN], graphV[GN];
double F;
 
 
int count = 0;
 
// プロトタイプ宣言
void Init();
void Calculate();
void Display();
void subDisplay();
void subDisplay2();
 
// 連立微分方程式
double fx(double t,double x,double v)
{
	return v;
}
double fv(double t,double x,double v, double omega, double omega0, double f)
{
	return -omega0 * omega0 * x + f * cos(omega * t);
}
 
void psMain()
{
	drps2.window.Size(800, 300);
	drps2.window.Position(0, 50);
	drps2.event.Init(Init);
	drps2.event.Calculate(Calculate);
	drps2.event.Display(Display);
	drps2.psCreateWindow();
 
	graph.window.Size(430, 650);
	graph.window.Position(620, 75);
	graph.event.Display(subDisplay);
	graph.psCreateWindow();
 
	phase.window.Size(620, 500);
	phase.window.Position(0, 225);
	phase.event.Display(subDisplay2);
	phase.psCreateWindow();
}
 
 
void Init()
{
	int i;
 
	t = 0.0;
	x = 3.0;
	v = 0.0;
 
	for(i=0;i<GN;i++){
		graphX[i] = 0.0;
		graphV[i] = 0.0;
	}
}
 
void Calculate()
{
	int i;
	double kx[4],kv[4];
 
	kx[0]=dt*fx(t,x,v);
	kv[0]=dt*fv(t,x,v,omega,omega0,f);
 
	kx[1]=dt*fx(t+dt/2.0,x+kx[0]/2.0,v+kv[0]/2.0);
	kv[1]=dt*fv(t+dt/2.0,x+kx[0]/2.0,v+kv[0]/2.0,omega,omega0,f);
 
	kx[2]=dt*fx(t+dt/2.0,x+kx[1]/2.0,v+kv[1]/2.0);
	kv[2]=dt*fv(t+dt/2.0,x+kx[1]/2.0,v+kv[1]/2.0,omega,omega0,f);
 
	kx[3]=dt*fx(t+dt,x+kx[2],v+kv[2]);
	kv[3]=dt*fv(t+dt,x+kx[2],v+kv[2],omega,omega0,f);
 
	t += dt;
	x += (kx[0]+2.0*kx[1]+2.0*kx[2]+kx[3])/6.0;
	v += (kv[0]+2.0*kv[1]+2.0*kv[2]+kv[3])/6.0;
 
	T = 0.5*m*v*v;
	U = 0.5*k*x*x;
	E = T + U;
 
        F = f*cos(omega*t);
 
	count++;
	if(count%30 == 0){	
		for(i=0;i<GN-1;i++){
			graphX[i] = graphX[i+1];
			graphV[i] = graphV[i+1];
		}
		graphX[GN-1] = x;
		graphV[GN-1] = v;
	}
}               
 
 
void Display()//
{
	drps2.EraseOrbit();
 
	//外力
	drps2.SetPen(PS_SOLID, 4, RGB(250, 0, 0));
	drps2.DrawVector(
		320.0+x*10.0,
		100.0,
		F*40.0,
		0.0
	);
 
	drps2.SetPen(PS_SOLID, 1, RGB(0, 0, 0));
 
	drps2.DrawLine(//壁
		100.0,
		80.0,
		100.0,
		120.0
	);
	drps2.DrawLine(//床
		100.0,
		120.0,
		520.0,
		120.0
	);
	drps2.DrawLine(//バネ
		100.0,
		100.0,
		300.0+x*10.0,
		100.0
	);
	drps2.DrawRectangle(//物体
		300.0+x*10.0,
		120.0,
		340.0+x*10.0,
		80.0
	);
}
 
void subDisplay()//x-t,v-tグラフ
{                                     
	int i;
	char str[64];
 
	graph.EraseOrbit();
 
	sprintf(str, "位置x : %f", x);
	graph.DrawText(2, 2, str, RGB(0, 0, 255));
	sprintf(str, "速度v : %f", v);
	graph.DrawText(2, 22, str, RGB(0, 255, 0));
 
	graph.SetPen(PS_SOLID, 2, RGB(0, 0, 0));
	//x_tグラフ
	sprintf(str, "x");
	graph.DrawText(17.5, 100.0, str, RGB(0, 0, 0));
	graph.DrawVector( //x
		20.0,
		270.0,
		0.0,
		-150.0
	);
	sprintf(str, "t");
	graph.DrawText(395.0, 190.0, str, RGB(0, 0, 0));
	graph.DrawVector( //x_t
		2.0,
		200.0,
		380.0,
		0.0
	);
 
	//v_tグラフ
	sprintf(str, "v");
	graph.DrawText(17.5, 300.0, str, RGB(0, 0, 0));
	graph.DrawVector( //v
		20.0,
		470.0,
		0.0,
		-150.0
	);
	sprintf(str, "t");
	graph.DrawText(395.0, 390.0, str, RGB(0, 0, 0));
	graph.DrawVector( //v_t
		2.0,
		400.0,
		380.0,
		0.0
	);
 
	//x
	graph.SetPen(PS_SOLID, 1, RGB(0, 0, 255));
	for(i=0;i<GN-1;i++){
		graph.DrawLine(
			(double)(i*2),
			200.0-graphX[i]*5.0,
			(double)((i+1)*2),
			200.0-graphX[i+1]*5.0
		);
	}
	//v
	graph.SetPen(PS_SOLID, 1, RGB(0, 255, 0));
	for(i=0;i<GN-1;i++){
		graph.DrawLine(
			(double)(i*2),
			400.0-graphV[i]*5.0,
			(double)((i+1)*2),
			400.0-graphV[i+1]*5.0
		);
	}
}
 
void subDisplay2()//位相空間描写
{
	char str[64];
 
	/*エネルギー*/
	sprintf(str, "T : %f", T);	//運動エネルギー
	phase.DrawText(2, 2, str, RGB(0, 0, 0));
	sprintf(str, "U : %f", U);      //ポテンシャル・エネルギー
	phase.DrawText(2, 22, str, RGB(0, 0, 0));
	sprintf(str, "E : %f", E);	//全エネルギー
	phase.DrawText(2, 42, str, RGB(0, 0, 0));
 
	/*横軸ｘ*/
      	phase.SetPen(PS_SOLID, 1, RGB(0, 0, 0));
	phase.DrawVector(
		20.0,
		220.0,
		550.0,
		0.0
	);
	sprintf(str, "x");
	phase.DrawText(580.0, 210.0, str, RGB(0, 0, 0));
 
        /*縦軸ｖ*/
	phase.DrawVector(
		320.0,
		400.0,
		0.0,
		-350.0
	);
	sprintf(str, "v");
	phase.DrawText(317.0, 32.0, str, RGB(0, 0, 0));
 
	/*軌跡描画*/
	phase.SetPen(PS_SOLID, 1, RGB(0, 0, 250));
	phase.DrawEllipse(
		320.0+x*13.0-2.2,
		220.0+v*13.0-2.2,
		320.0+x*13.0+2.2,
		220.0+v*13.0+2.2
		);
}