0477-1992 东大工 数学 WEB

常微分方程 差分方程 数值分析 稳定性分析

2階の線形常微分方程式

$$\frac{{\mathrm{d}}^{2}x}{\mathrm{d}{t}^2}+p^{2}x=0\text{\hspace{1em}(e)}$$
の解を数値的に求めるための差分（階差）方程式として以下のものを考える。
$$x_{n+1}-(2-{\alpha }^2)x_n+x_{n-1}=0\text{\hspace{1em}(f)}$$
ここで、$p$ は正の定数、$\alpha$ はパラメータ（$\beta \left(0\le \beta <\frac{1}{4}\right)$）に依存する正の数であり、
$${\alpha }^2=\frac{{\theta }^2}{1+\beta {\theta }^2},{\theta }^2=p^2{h}^2$$
で与えられる。また、$h$ は差分化する際の時間幅である（$t=nh$）。

(1) (f) 式の解として、

$$x_n=C_1{e}^{{\lambda }_{1}n}+C_2{e}^{{\lambda }_{2}n}\text{\hspace{1em}(g)}$$
の形のものを考える。ただし、$C_1,C_2$ は同時に $0$ ではない任意の定数、${\lambda }_1,{\lambda }_2$ は $\alpha$ によって定まる相異なる定数である。解 (g) が $n\to \pm \infty$ で発散しないためには、$\alpha$ はどのような範囲にあればよいか。また、$\alpha$ がこの範囲にあるとき、時間幅 $h$ と (e) 式を解析的に解いて得られる解の振動周期 $T$ との比 $\frac{h}{T}$ の値の範囲を $\beta$ を用いて表せ。

(2) $\alpha$ が (1) で与えられた範囲にあるとき、適当な初期条件のもとで (f) 式の解を数値的に求めよ。得られる解の周期 $\tilde{T}$ の相対誤差 $\frac{\tilde{T}-T}{T}$ を求めよ。

---

解答：
(1)
(f) 式に $x_n=e^{\lambda n}$ を代入すると、以下の特性方程式を得る。

$$e^{2\lambda }-(2-{\alpha }^2)e^{\lambda }+1=0$$
$z=e^{\lambda }$ とおくと、
$$z^2-(2-{\alpha }^2)z+1=0$$
この2次方程式の解を $z_1,z_2$ とする。解 $x_n=C_1{z}_1^n+C_2{z}_2^n$ が $n\to \pm \infty$ で発散しないためには、$|z_1|=|z_2|=1$ かつ $z_1\ne z_2$ である必要がある。
根と係数の関係より $z_1{z}_2=1$ であるため、$z_1,z_2$ は互いに共役な複素数とならなければならない。
したがって、判別式 $D=(2-{\alpha }^2{)}^2-4<0$ を満たす必要がある。
$${\alpha }^4-4{\alpha }^2<0$$
$\alpha >0$ であるから、
$$\boxed{0<\alpha <2}$$

(e) 式の解析解の振動周期は $T=\frac{2\pi }{p}$ である。$\theta =ph$ より、

$$\frac{h}{T}=\frac{hp}{2\pi }=\frac{\theta }{2\pi }$$
与えられた ${\alpha }^2=\frac{{\theta }^2}{1+\beta {\theta }^2}$ を ${\theta }^2$ について解くと、
$${\theta }^2=\frac{{\alpha }^2}{1-\beta {\alpha }^2}$$
$0\le \beta <\frac{1}{4}$ および $0<\alpha <2$ より、$1-\beta {\alpha }^2>1-\frac{1}{4}\cdot 4=0$ であるため、${\theta }^2$ は ${\alpha }^2\in (0,4)$ において単調増加する。
$$0<{\theta }^2<\frac{4}{1-4\beta }$$
$$0<\theta <\frac{2}{\sqrt{1-4\beta }}$$
したがって、$\frac{h}{T}$ の値の範囲は、
$$\boxed{0<\frac{h}{T}<\frac{1}{\pi \sqrt{1-4\beta }}}$$

(2)
$\alpha \in (0,2)$ のとき、特性方程式の解は共役複素数となり、

$$z=\frac{2-{\alpha }^2\pm i\alpha \sqrt{4-{\alpha }^2}}{2}$$
$z=e^{\pm i\varphi }$ とおくと、実部を比較して
$$\cos \varphi =1-\frac{{\alpha }^2}{2}(0<\varphi <\pi )$$
数値解は、初期条件で定まる任意定数を $A,B$ として、次のように表される。
$$x_n=A\cos (\varphi n)+B\sin (\varphi n)$$
この解のステップ数 $n$ に対する周期は $N=\frac{2\pi }{\varphi }$ であり、実際の時間に対する周期 $\tilde{T}$ は
$$\tilde{T}=Nh=\frac{2\pi h}{\varphi }$$
解析解の周期は $T=\frac{2\pi h}{\theta }$ であるから、相対誤差は
$$\frac{\tilde{T}-T}{T}=\frac{\tilde{T}}{T}-1=\frac{\frac{2\pi h}{\varphi }}{\frac{2\pi h}{\theta }}-1=\frac{\theta }{\varphi }-1$$
以上より、
$$\boxed{\begin{array}{rl}{x}_n& =A\cos (\varphi n)+B\sin (\varphi n)\left(A,B\text{ は任意定数、}\varphi =\arccos \left(1-\frac{{\alpha }^2}{2}\right)\right)\\ \frac{\tilde{T}-T}{T}& =\frac{\theta }{\arccos \left(1-\frac{{\alpha }^2}{2}\right)}-1\end{array}}$$

---

本题探讨了二阶线性常微分方程（简谐振子方程）的数值求解及其稳定性与误差分析。
第一问首先要求找出使得差分方程数值解不发散的条件。通过特征方程可以将差分方程转化为代数方程求解。由于物理上要求解不仅不发散，而且不衰减（即保持等幅振荡），因此特征根必须严格位于复平面的单位圆上。据此可以求出参数的取值范围。接着利用题目给定的参数关系式，解出单步时间步长对应的相角上限，从而得出步长与物理周期之比的允许范围。这是典型的时间积分算法（如Newmark-$\beta$方法）中关于条件稳定性的推导过程。
第二问要求写出数值解的形式并计算周期的相对误差。数值解的表达式可以直接通过将特征根写成复指数形式，利用欧拉公式展开得到实数域的三角函数解。数值解的周期与解析解的真实周期的偏差来源于离散化造成的频率畸变现象。通过对比离散相角与真实相角的比值，即可得出周期的相对误差，该误差指标通常用于评估数值积分算法在长时间模拟动力学问题中的相位保真能力。