数学——Euler方法求解微分方程详解（python3）

算法的数学描述图解

实例

用Euler算法求解初值问题 $$ \frac{dy}{dx}=y+\frac{2x}{y^2}$$ 初始条件$y(0)=1$，自变量的取值范围$x \in [0, 2]$

算法Python3代码求解

``` # 导入包 import numpy as np import matplotlib.pyplot as plt # 定义求解函数 y_dot = y + 2*x/(y*y) def fx(y, x): return y + 2*x/(y*y) # 算法定义 def ode_euler(f, y0, tf, h): """ Solve and ODE using Euler method. Solve the ODE y_dot = f(y, t) Parameters ------------ :param f: function Function describing the ODE :param y0: array_like Initial conditions. :param tf: float Final time. :param h: float Time step :return: y : array_like Solution to the ODE. t : array_like Time vector. """

y0 = np.array(y0)
ts = np.arange(0, tf + h, h)
y = np.empty((ts.size, y0.size))
y[0, :] = y0
for t, i in zip(ts[1:], range(ts.size - 1)):
    y[i + 1, :] = y[i, :] + h * f(y[i, :], t)
return y, ts

实例应用案例

def newton_cooling_example():
print('Solving Newton Cooling ODE...')
y, ts = ode_euler(fx, 1, 2, 0.01)
print('Done.')
plt.figure()
plt.plot(ts, y)
plt.xlabel('time [s]')
plt.title('Solution to the Newton cooling equation')
plt.show()

<h2>代码中的部分函数理解</h2>
<h3>numpy.array</h3>
`numpy.array(object, dtype=None, copy=True, order='K', subok=False, ndmin=0)`
[参考numpy.array]( https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.array.html)
output：**创建一个array**，返回类型为ndarray
<strong>实例</strong>

np.array([1, 2, 3.0]) # array([1., 2., 3.])
np.array([[1, 2], [3, 4]]) # array([[1, 2], [3, 4]])
np.array([1, 2, 3], dtype=complex) # array([1.+0.j, 2.+0.j, 3.+0.j])

<h3>numpy.arange</h3>
[参考numpy.arange](https://docs.scipy.org/doc/numpy/reference/generated/numpy.arange.html)
`numpy.arange([start, ]stop, [step, ]dtype=None)`
**作用**：在给定间隔内返回均匀间隔的值。
值在半开区间[start, stop)内生成（换句话说，包括开始但不包括终止）。返回的是ndarray而不是列表。
np.arange()函数返回一个有终点和起点的固定步长的排列，如[1,2,3,4,5]，起点是1，终点是5，步长为1。 
参数个数情况： np.arange()函数分为一个参数，两个参数，三个参数三种情况 ：

    1. 一个参数时，参数值为终点，起点取默认值0，步长取默认值1。 
    2. 两个参数时，第一个参数为起点，第二个参数为终点，步长取默认值1。 
    3. 三个参数时，第一个参数为起点，第二个参数为终点，第三个参数为步长。其中步长支持小数。

**案例**

np.arange(3,7) # array([3, 4, 5, 6])
np.arange(3,7,2) # array([3, 5])

<h3>numpy.ma.size</h3>
`numpy.ma.size(obj, axis=None)`
[参考](https://docs.scipy.org/doc/numpy-1.10.1/reference/generated/numpy.ma.size.html#numpy.ma.size)
**案例**

a = np.array([[1,2,3],[4,5,6]])
np.size(a) # 6
np.size(a,1) # 3
np.size(a,0) # 2

<h3>numpy.empty</h3>
[参考](https://docs.scipy.org/doc/numpy/reference/generated/numpy.empty.html)
`numpy.empty(shape, dtype=float, order='C')`
**shape** : int or tuple of int  Shape of the empty array, e.g., (2, 3) or 2.
**out** : ndarray
**案例**

np.empty([2, 2])

结果

array([[ -9.74499359e+001, 6.69583040e-309],
[ 2.13182611e-314, 3.06959433e-309]]) #random
np.empty([2, 2], dtype=int)

结果

array([[-1073741821, -1067949133],
[ 496041986, 19249760]]) #random