소개

파이썬을 사용하여 상미분 방정식을 푸는 방법을 설명합니다.

OS/Python module version

사용한 OS는 Windows 10입니다.
파이썬 모듈 버전은 다음과 같습니다.

콘솔

> python --version
Python 3.4.3
> pip freeze
matplotlib==1.4.3
numpy==1.13.1+mkl
pandas==0.16.2
pyparsing==2.0.3
python-dateutil==2.4.2
pytz==2015.4
scipy==0.19.1
six==1.9.0

샘플 코드

다음 미분 방정식을 푸십시오.

y'' = -y \\

변수 z를 도입합니다.
변수를 늘림으로써 1개의 2층 미분 방정식을 2개의 1층 미분 방정식의 연립 방정식으로 생각합니다.

y' = z \\
z' = -y \\

초기 조건을 $t=0, y_0 = 1, z_0 = 0$로 합니다.
솔루션은 $cos(t)$입니다. 올바른 해가 용이하게 요구되기 때문에 오차의 평가를 할 수 있습니다.

미분 방정식을 푸는 Python 프로그램은 다음과 같습니다.
scipy.integrate.ode를 사용하여 해결합니다.

test.py

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import ode

def f(t, v): # t, y:v[0], y'=z:v[1]
    return [v[1], -v[0]] # y':return[0] y''=z':return[1]
v0 = [1.0, 0.0] # y0, y'0=z0

solver = ode(f)
solver.set_integrator(name="dop853")
solver.set_initial_value(v0)
tw = 10.0*2.0*np.pi
dt = tw / 1000;
t = 0.0
ts = []
ys = []
while solver.t < tw:
    solver.integrate(solver.t+dt)
    ts += [solver.t]
    ys += [solver.y[0]]
plt.figure(0)
plt.plot(ts, ys)
plt.figure(1)
plt.plot(ts, np.cos(ts)-ys)
plt.show()
print(np.max(np.cos(ts)-ys))

실행합니다.

콘솔

> python .\test.py

실행하면 그래프가 표시됩니다.

다음 그래프는 미분 방정식의 해입니다. cos(t)의 그래프가 10주기분 표시되어 있습니다.


다음 그래프는 오차를 플로팅합니다. t가 증가함에 따라 오차도 증가하고 있습니다.
오차는 $10^{-15}$의 스케일이 됩니다.