解微分方程

In [1]:

  1. %pylab inline
  1. Populating the interactive namespace from numpy and matplotlib

积分求解

简单的例子

\frac{dy}{dt} = sin(t)

In [2]:

  1. def dy_dt(y, t):
  2.     return np.sin(t)

积分求解：

In [3]:

  1. from scipy.integrate import odeint
  2.  
  3. t = np.linspace(0, 2*pi, 100)
  4.  
  5. result = odeint(dy_dt, 0, t)

In [4]:

  1. fig = figure(figsize=(12,4))
  2. p = plot(t, result, "rx", label=r"$\int_{0}^{x}sin(t) dt $")
  3. p = plot(t, -cos(t) + cos(0), label=r"$cos(0) - cos(t)$")
  4. p = plot(t, dy_dt(0, t), "g-", label=r"$\frac{dy}{dt}(t)$")
  5. l = legend(loc="upper right")
  6. xl = xlabel("t")

04.07 解微分方程 - 图1

高阶微分方程

抛物运动（竖直方向）：

\frac{d^2x}{dt^2} = g - \frac{D}{m}\frac{dx}{dt}

改写成如下形式：

y = \left[x, \frac{dx}{dt}\right]

\begin{aligned}\frac{dy_0}{dt} &= y_1 \\frac{dy_1}{dt} &= -g - \frac{D}{m} y_1 \\end{aligned}

In [5]:

  1. def dy_dt(y, t):
  2.     """Governing equations for projectile motion with drag.
  3.     y[0] = position
  4.     y[1] = velocity
  5.     g = gravity (m/s2)
  6.     D = drag (1/s) = force/velocity
  7.     m = mass (kg)
  8.     """
  9.     g = -9.8
  10.     D = 0.1
  11.     m = 0.15
  12.     dy1 = g - (D/m) * y[1]
  13.     dy0 = y[1] if y[0] >= 0 else 0.
  14.     return [dy0, dy1]

In [6]:

  1. position_0 = 0.
  2. velocity_0 = 100
  3. t = linspace(0, 12, 100)
  4. y = odeint(dy_dt, [position_0, velocity_0], t)

In [7]:

  1. p = plot(t, y[:,0])
  2. yl = ylabel("Height (m)")
  3. xl = xlabel("Time (s)")

04.07 解微分方程 - 图2

In [8]:

  1. y, infodict = odeint(dy_dt, [position_0, velocity_0], t, full_output=True, printmessg=True, )
  2. print sorted(infodict.keys())
  3. print "cumulative number of function evaluations at each calculated point:", infodict['nfe']
  4. print "cumulative number of time steps", infodict['nst']
  1. Integration successful.
  2. ['hu', 'imxer', 'leniw', 'lenrw', 'message', 'mused', 'nfe', 'nje', 'nqu', 'nst', 'tcur', 'tolsf', 'tsw']
  3. cumulative number of function evaluations at each calculated point: [ 45  49  51  53  55  59  61  61  63  65  67  67  69  71  73  73  75  77
  4.   77  79  79  81  81  83  85  85  87  87  89  89  91  91  93  95  95  97
  5.   97  99  99 101 101 103 103 105 107 107 109 109 111 111 113 113 115 115
  6.  117 117 119 119 121 121 123 123 123 125 125 127 127 129 129 131 131 131
  7.  133 133 135 135 135 137 137 139 139 139 141 141 143 143 143 145 145 147
  8.  147 149 149 149 154 158 274 280 280]
  9. cumulative number of time steps [ 20  22  23  24  25  27  28  28  29  30  31  31  32  33  34  34  35  36
  10.   36  37  37  38  38  39  40  40  41  41  42  42  43  43  44  45  45  46
  11.   46  47  47  48  48  49  49  50  51  51  52  52  53  53  54  54  55  55
  12.   56  56  57  57  58  58  59  59  59  60  60  61  61  62  62  63  63  63
  13.   64  64  65  65  65  66  66  67  67  67  68  68  69  69  69  70  70  71
  14.   71  72  72  72  73  75 130 133 133]

原文: https://nbviewer.jupyter.org/github/lijin-THU/notes-python/blob/master/04-scipy/04.07-ODEs.ipynb