【数値解析】最終レポート1.オイラー法

数値解析のレポートをやります。これが終われば春休みになったも同然です。はい。
ただ、授業聞いてなかったので、まずはInputからです。はい。
問題はこれ

問.調和振動子の初期値問題: $$ { \displaystyle \frac{d}{dt} \left( \begin{array}{c} p \\\ q \end{array} \right) =J \left( \begin{array}{c} p \\\ q \end{array} \right), \hspace{1cm} J = \left( \begin{array} {cc} 0 & -1 \\\ 1 & 0 \end{array} \right) } $$ を数値的に解くことを考える。(ただし、 $t>0$ ) (1) この問題を陽的Euler法、陰的Euler法、台形則、4次Runge-Kutta法などで解き、解軌道を(p,q)平面状にかけ(時間の刻み幅は、各種法の解軌道の特徴が現れるように適切に選ぶこと。)

(2)上記の結果は、講義で説明した安定性条件からどのように理解できるか議論せよ。(数値実験できない場合は理論的に何が予想されるか議論せよ)

(3)もあったけど、問題が結構変わるので、今は割愛しました。
今回はNumpyを使ってやっていくことに。
そのためにまず、Numpyの知識を仕入れる。
(この過程もこのページに書いてると長くなりそうなので、違うページに書きました)

ってなわけで、次は、プログラムを組んで行きますよ〜
Pythonによるオイラー法の参考にしたのはこちら

特に、後者にはめっちゃお世話になりました...
また、厳密解が $p^{2}+q^{2}=1$ だったので、円を描写するべく以下のページも参考に。

Matplotlibで円や長方形などの図形を描画

道具は揃ったじゃコード組みますか！！！！

実際組み上がったコード

がこちら

import matplotlib.pyplot as plt
import matplotlib.patches as patches #楕円を書くため
import numpy as np

J = np.array([[0.00,-1.00],[1.00,0.00]])
#微分方程式右辺
def dxdt(x):
    return J.dot(x)

def EulerOrbit(x0,MAX_T,dt):
    x = x0.copy() #vector #ndarray
    t = 0.00
    orbit = []
    while t < MAX_T:
        orbit.append(list(np.round(x,2)))
        x += dxdt(x)*dt #ndayyay
        t += dt
    return orbit


dt = 0.1 #刻み幅
T = 100.00 #試行上限
x0 = np.array([1.00,0.00]) # t=0でのx

orbit = EulerOrbit(x0,T,dt)
nporbit = np.array(orbit) #listが、ndarray型に

p = nporbit[:, 0]
q = nporbit[:, 1]
plt.plot(p,q)

# 厳密解のプロット
#t = np.linspace(0, MAX_T, 1/DELTA_T) #等差数列の生成
#円を描くための準備
ax = plt.axes()
# fc = face color, ec = edge color
x = patches.Circle(xy=(0, 0), radius=1.0, fill=False, ec='r')
ax.add_patch(x) #描写

plt.xlabel($p$,fontsize=14)
plt.ylabel($q$,fontsize=14)
plt.show()

このコード生むのにめちゃくちゃ時間かかった...　
途中で

File "EulerMethod.py", line 17
    x += dxdt(x)*dt #ndayyay
    ^
SyntaxError: invalid syntax

というErrorにハマって...、for文の中身を1行ずつやったら通ったり、通らなかったり...
「マジでなんでだろう...」と思ってたら、この前の行のorbit.append(list(np.round(x,2)))の最後の)が抜けてるっていう、あるあるなしょうもないミスでした。はい。(キレそう)

よっしゃオイラー法終わり！と思ったら、陽的オイラー法が終わっただけでした!!!(キレ気味)
陰的は...わからん！調べます！！

雑談

ちなみに2こめの参考記事であるPDFを読んでる時に、

vecx = list(map(lambda v: v(t, *x), vectorfield))

という１行が出てきて「！？！？」となったPython弱者僕は、map関数なるものの存在を初めて知りました。
その時読んだ記事がこれ

高階関数...便利そう！！

雑談その２

今回の記事をかくに当たって、そもそも、はてなブログで数式描く術を知らなかったので、色々調べました(笑)
その時参考に下のがこれらの記事

これもあって時間かかりました(笑)

次のステージ（）へ

台形則やりますか！！！(空元気)

ナルミヤの備忘録(仮)

ナルミヤが学んだことなどを書き記していくブログ(方向性模索中。)

【数値解析】最終レポート1.オイラー法

実際組み上がったコード

雑談

雑談その２

次のステージ（）へ