单摆运动数值模拟

一、单摆运动

单摆是能够产生往复摆动的一种装置，将无重细杆或不可伸长的细柔绳一端悬于重力场内一定点，另一端固结一个重小球，就构成单摆。若小球只限于铅直平面内摆动，则为平面单摆，若小球摆动不限于铅直平面，则为球面单摆。
重力对单摆的力矩为：
$$m{l}^2\theta +\gamma l\theta +mgsin(\theta )=Fcos({\omega }_{D}t)$$

令${\omega }_D=\sqrt{\frac{g}{l}}$

$$\theta +\gamma \frac{\theta }{m}+{\omega }_0^{2}sin(\theta )=\frac{F}{ml}cos({\omega }_{D}t)$$

小角度近似得
$$\theta +2\beta \theta +{\omega }_0^{2}sin()\theta =fcos({\omega }_{D}t)$$

数值处理
令
$$\frac{d\theta }{dt}=\omega$$
则
$$\frac{d\omega }{dt}+\frac{g}{l}sin(\theta )=0$$

根据有限差分法可得
$$\frac{{\omega }_{n+1}-{\omega }_n}{\delta t}+\frac{g}{l}sin({\theta }_n)=0$$

二、代码转换

#include <stdio.h>
#include <math.h>
#include <stdlib.h>

void put(FILE *outstream, double data1, double data2, double data3)
{
	fprintf(outstream, "%f,%f,%f\n,", data1, data2, data3);
} 

int main(int argc, char *argv[])
{
	double theta = 0.02;//0.02
	double w = 0;
	double g=9.8;
	double l = 1;
	double beta = 0.05;//0.05
	double w_d = 3.0;
	double f = 4.25;//7.3185
	float t = 0.01;
	int i;
	
	FILE *outstream;
	if((outstream=fopen("hahaha.dat", "ab")) == NULL)
	{
		printf("Cannot open file, press any key to exit!\n");
		exit(1);
	}
	
	for(i=0; i<10000; i++)
	{
		w = w +(-g/l*sin(theta)-2*beta*w+f*cos(w_d*t*i))*t;
		theta = theta + w*t;
		put(outstream,w,theta,t*(i));
	}
	fclose(outstream);
	printf("%f\n",sqrt(g/l));
	return 0;
}

三、 结果分析

（一）无阻尼无驱动

（二）附加阻尼

（三）附加阻尼和驱动力

f=1

f=6

f=5

f=7

四、总结

总体的思路非常简单，对于一阶常微分方程，差分法使用起来非常方便，对于二阶，通常可将二阶降为一阶进行计算，但是降阶后对于代码的顺序需要调整，即把先要用到的结果放在前面，总的来说算是对差分法的回顾吧。