Python で解く常微分方程式とグラフ

    numpy と matplotlib.pylot の import を行っておく。
    >>> import numpy
    >>> import matplotlib.pyplot as plt
  
□ 関数の定義と作図

    >>> def f(x):
            return(x**2 - 4*x + 6)
    >>> f(2)
    >>> f(numpy.array([2, 3]))

    >>> t = numpy.linspace(1, 3, 101)
    >>> plt.plot(t, f(t))
    >>> plt.grid('on')
    >>> plt.show()

    >>> t = numpy.logspace(-1, 2, 100)
    >>> plt.semilogx(t, f(t))
    >>> plt.grid('on')
    >>> plt.show()
    >>> plt.semilogy(t, f(t))
    >>> plt.grid('on')
    >>> plt.show()
    >>> plt.loglog(t, f(t))
    >>> plt.grid('on')
    >>> plt.show()

□ 常微分方程式の解

   x'' + 4*x' + 13*x = 0, x'(0) = 1, x(0) = 1
   の解を求める．一階の常微分方程式に書き改めると次のようになる．

   x1' = x2
   x2' = -4*x2 - 13*x1
    
   微分方程式の数値階の計算には scipy.integrate.odeint を用いるため、
   import しておく。

    >>> from scipy.integrate import odeint
    >>> def xdot(x, t):
            return (x[1], -4*x[1] - 13*x[0])
    >>> x0 = numpy.array([1, 1])
    >>> t = numpy.linspace(0, 4, 1000)
    >>> y = odeint(xdot, x0, t)

    >>> plt.plot(t, y)
    >>> plt.grid('on')
    >>> plt.show()
    >>> plt.plot(t, y[:, 0])
    >>> plt.grid('on')
    >>> plt.show()
    >>> plt.plot(t, y[:, 1])
    >>> plt.grid('on')
    >>> plt.show()
    >>> plt.plot(y[:, 0], y[:, 1])
    >>> plt.grid('on')
    >>> plt.show()

   強制振動系
   x'' + 2*x' + (4*π)^2*x = sin(4πt), x'(0) = 0, x(0) = 0
   の解を求める．一階の常微分方程式に書き改めると次のようになる．

   x1' = x2
   x2' = -2*x2 - (4*π)^2*x1 + sin(4πt)

    >>> def xdot(x, t):
            return (x[1], \
                    -2*x[1] - (4*numpy.pi)**2*x[0] + numpy.sin(4*numpy.pi*t))
    >>> x0 = numpy.array([0, 0])
    >>> t = numpy.linspace(0, 4, 1000)
    >>> y = odeint(xdot, x0, t)

    >>> plt.plot(t, y)
    >>> plt.grid('on')
    >>> plt.show()
    >>> plt.plot(y[:, 0], y[:, 1])
    >>> plt.grid('on')
    >>> plt.show()

□ ３次元プロット

    mpl_toolkits.mplot3d.Axes3 を使うので、import する。

    >>> from mpl_toolkits.mplot3d import Axes3D

   0 ≦ x ≦ 4 and 0 ≦ y ≦ 4 での sin(x)*sin(y) の図

    >>> x = numpy.linspace(0, 4, 40)
    >>> y = x
    >>> X, Y = numpy.meshgrid(x, y)
    >>> Z = numpy.sin(X)*numpy.cos(Y)
    >>> fig = plt.figure(figsize = (8,8))
    >>> ax = fig.add_subplot(111, projection='3d')
    >>> ax.plot_surface(X, Y, Z)
    >>> plt.show()

   -π ≦ x ≦ π and -π/2 ≦ y ≦ π/2 での sin(x*y) の図

    >>> x = numpy.linspace(-numpy.pi, numpy.pi, 40)
    >>> y = numpy.linspace(-numpy.pi/2, numpy.pi/2, 40)
    >>> X, Y = numpy.meshgrid(x, y)
    >>> Z = numpy.sin(X * Y)
    >>> fig = plt.figure(figsize = (8,8))
    >>> ax = fig.add_subplot(111, projection='3d')
    >>> ax.plot_surface(X, Y, Z)
    >>> plt.show()