第八次作业

作者：敬雷

L1

3.4题

---

1

不想放代码 :(

待解方程：

$$\frac{d^2x}{dt^2} = -kx^{\alpha}$$
可化为：
$$\frac{dx}{dt} = v, \,\,\frac{dv}{dt} = -kx^{\alpha}$$
$\mathrm{Euler}$-$\mathrm{Cromer}$ 方法：

• $v(t + dt) = v(t) - kx^{\alpha}dt$
• $x(t + dt) = x(t) + v(t + dt)dt$

---

2

$\alpha = 1$，$k = 1$，即解方程 $d^2x / dt^2 = -x$。
初始条件设为 $x = A, dx/dt = 0$。其解析解为 $x = A\cos(t)$。
取$dt = 0.05$，用 $\mathrm{Euler}$-$\mathrm{Cromer}$ 方法得到的解的图像如下。

可以看出， $\mathrm{Euler}$-$\mathrm{Cromer}$ 得到的解非常稳定，振幅几乎没有变化。作为对比，红色虚线表示的 $\mathrm{Euler}$ 方法得到的解振幅一直在增加。

为了进一步说明 $\mathrm{Euler}$-$\mathrm{Cromer}$ 方法的稳定性，设振子质量为 1，计算其机械能：

可以看出，机械能在12.2至12.8之间波动，在振子的一个振动周期内波动两次。并且即使总时间增加到了30000秒，机械能仍然没有明显的变化。
有趣的是，$dt$取得很大时，比如 1.1s，振幅也不会一直增加或减小，而是周期性变化（下图）。这说明 $\mathrm{Euler}$-$\mathrm{Cromer}$ 方法确实吊爆了。

---

3

$\alpha = 1$，$k = 1$，即解方程 $d^2x / dt^2 = -x$。
初始条件设为 $x = A, dx/dt = 0$。其解析解为 $x = A\cos(t)$。
取$dt = 0.05$，对不同的振幅 $A$，用 $\mathrm{Euler}$-$\mathrm{Cromer}$ 方法得到的解的周期如下。

可以看出，周期恒为 6.28s，确实不随振幅变化。

---

4

$\alpha = 3$，$k = 1$，即解方程 $d^2x / dt^2 = -x^3$。
初始条件设为 $x = A, dx/dt = 0$。
取$dt = 0.05$，用 $\mathrm{Euler}$-$\mathrm{Cromer}$ 方法得到的解的图像和周期随振幅变化的图像。

可以看出，这个解和简谐振动挺像的，都是周期运动。而振动周期随振幅增加而减小。
一个直观的解释是，振幅很小时，其受力比简谐运动小很多，于是运动得慢，周期很大；振幅很大时，受力比简谐运动大很多，于是运动的快，周期很小。

---

L2

3.5题

---

无代码可放 :)

方程 $\frac{d^2x}{dt^2} = -kx^{\alpha}$ 的周期为：

$$T = K(\alpha)\sqrt{\frac{2(\alpha+1)}{k}} A^{\frac{1-\alpha}{2}}$$

其中$K(\alpha) = \frac{\sqrt{\pi}\Gamma(1 + \beta)} {\Gamma(1/2 + \beta)} + F(1/2, \beta, 1 + \beta, e^{i\pi/\beta})$，$\beta = \frac{1}{a+\alpha}$，$F$ 为超几何级数，$A$ 为振幅。取 $k = 1$，计算可得当 $\alpha = 1$ 时，$T = 2\pi$，与振幅$A$无关；而 $\alpha = 3$ 时，$T \approx 7.42/A$ 。程序算出来的和理论值一致。
.

.
.
.
.
..

.
.
.
.
.

.
.。
.
.
..
.
.
.
.
.

完