![[PukiWiki]](image/ROBO192.png "[PukiWiki]")
#author("2023-06-10T17:45:17+09:00","external:moriat","moriat") #author("2023-06-10T17:46:21+09:00","external:moriat","moriat") #topicpath ** 数値計算法 [#j01c23d1] *** お題 [#ke057a34] - シミューレートする方程式 $$ \begin{eqnarray} \frac{dx}{dt} & = & v \\ \frac{dv}{dt} & = & -k x \end{eqnarray} $$ *** Euler 法 [#vf8465e4] - 概要~ 普通に差分化して計算する方法はオイラー法と呼ばれている。 #code(python){{ def Euler(x) : xn = x[0] vn = x[1] dt = 0.01 k = 3 xnn = xn + vn * dt vnn = vn - k * dt * xn return [xnn, vnn] }} ところが、これを実行すると発散する。 #code(python){{ x0 = [2,3] for i in range(10000): x1 <- Euler(x0) x0 <- x1 print x0 }} - [[Google Colaboratory でのプログラムと実行結果>https://colab.research.google.com/drive/1PsdppFVgM6wJhwI3dU7zmJAwOxPqFyM2#scrollTo=63rA4El90-wk&line=1&uniqifier=1]] *** Leap Frog 法 [#d94c05bb] - 概要~ $$ \begin{eqnarray} x_{n+1} &=& x_n + v_{n+1/2} \times {\Delta} t + O(\Delta t^3) \end{eqnarray} $$ となり、$ \Delta t^2 $の項が消える。これは、 $$ \begin{eqnarray} x_{n+1} &=& x_{n+1/2} + x'_{n+1/2} \frac{\Delta t}{2} + \frac{1}{2} x''_{n+1/2} \left(\frac{\Delta t}{2}\right)^2 + \cdots \\ x_{n} &=& x_{n+1/2} - x'_{n+1/2} \frac{\Delta t}{2} + \frac{1}{2} x''_{n+1/2} \left(\frac{\Delta t}{2}\right)^2 + \cdots \end{eqnarray} $$ の2式の間で引き算することで得られる。より高次の項は、$ \Delta t^3 $ から始まる。~ &ref(IMG_3128.jpg,,30%); - コード~ この計算では発散しない。 #code(python){{ def LeapFrog(x) : xn = x[0] vn = x[1] dt = 0.01 k = 3 xnn = xn + vn * dt vnn = vn - k * dt * xnn return [xnn, vnn] x0 = [2,3] for i in range(100000): x1 = LeapFrog(x0) x0 = x1 print(x0) }} - [[Google Colaboratory でのプログラムと実行結果>https://colab.research.google.com/drive/1PsdppFVgM6wJhwI3dU7zmJAwOxPqFyM2#scrollTo=CEQHY9fO3TUl&line=1&uniqifier=1]] *** 4次の Runge-Kutta法 [#fd5f0bac] - 参考URL~ -- [[倭算数理研究所>https://wasan.hatenablog.com/entry/2013/11/02/201200]] - 概要~ $$ \begin{eqnarray} \frac{dx}{dt} &=& f(x) \\ x_{n+1} & = & x_{n} + \frac{\Delta t}{6} ( k_1+2k_2 + 2k_3 + k_4)\\ \ & & \\ k_1 & = & f(x_n ) \\ k_2 & = & f\left(x_n + \frac{k_1}{2} \Delta t \right) \\ k_3 & = & f\left(x_n + \frac{k_2}{2} \Delta t \right)\\ k_4 & = & f\left(x_n + {k_3} \Delta t \right) \end{eqnarray} $$ - コード~ #code(Python){{ def RungeKutta(x): xn = x[0] vn = x[1] dt = 0.01 k = 3 k1x = vn k1v = -k * xn k2x = vn + 0.5 * dt * k1v k2v = -k * (xn + 0.5 * dt * k1x) k3x = vn + 0.5 * dt * k2v k3v = -k * (xn + 0.5 * dt * k2x) k4x = vn + dt * k3v k4v = -k * (xn + dt * k3x) xnn = xn + (dt / 6.0) * (k1x + 2 * k2x + 2 * k3x + k4x) vnn = vn + (dt / 6.0) * (k1v + 2 * k2v + 2 * k3v + k4v) return [xnn, vnn] x0 = [2, 3] for i in range(10000): x1 = RungeKutta(x0) x0 = x1 print(x0) }} - [[Google Colaboratory でのプログラムと実行結果>https://colab.research.google.com/drive/1PsdppFVgM6wJhwI3dU7zmJAwOxPqFyM2#scrollTo=fFxyFkuu5e4z&line=22&uniqifier=1]]
