2020-12-17

复杂数学系统 Seminar 8 的解释说明：

---

 首先我们附上源代码：

import matplotlib.pyplot as plt
import numpy as np


def seminar8(scenario=1):
    # seminar8 Exploration of integration methods
    # Here we compare four integration schemes: euler, RK2, midpoint quadrature
    # and RK4.
    # 比较四种积分方案：欧拉、RK2、中点求积、RK4

    if scenario == 1:
        myfunc = lambda x, t: f1(x, t)
        tmax = 400  # why???
        h = 2  # 0.1  1  2
        t = np.arange(0, tmax + h, h)
        x = np.zeros((len(t), 1))  # pre-allocate x, x0=0 预先分配
        x[0] = 0
    else:
        myfunc = lambda x, t: f2(x, t)
        tmax = 10  # why???
        h = 2  # h = 0.1  1  2
        t = np.arange(0, tmax + h, h)
        x = np.zeros((len(t), 1))  # pre-allocate x, x0=0;
        x[0] = 0.5

    # Euler integration 欧拉积分
    for i in np.arange(len(t) - 1):
        x[i + 1] = x[i] + h * myfunc(x[i], t[i])

    # Plot result
    plt.figure()
    plt.plot(t, x, 'k-')

    # Now use RK2 integration
    for i in np.arange(len(t) - 1):
        xmidpoint = x[i] + h / 2. * myfunc(x[i], t[i])  # why we use a float here?and
        x[i + 1] = x[i] + h * myfunc(xmidpoint, t[i] + h / 2.)
    plt.plot(t, x, 'b-')

    # Now use mid-point quadrature formula, that is:
    # x(t+h)=x(t)+h/2 f(x(t),t) +h/2 f(x(t+h/2),t+h/2)
    for i in np.arange(len(t) - 1):
        xmidpoint = x[i] + h / 2. * myfunc(x[i], t[i])
        x[i + 1] = x[i] + h / 2. * myfunc(x[i], t[i] + h / 2) + h / 2. * myfunc(xmidpoint, t[i] + h / 2.)
    plt.plot(t, x, 'r-')

    # Now use RK4 integration
    for i in np.arange(len(t) - 1):
        k1 = h * myfunc(x[i], t[i])
        k2 = h * myfunc(x[i] + k1 / 2., t[i] + h / 2.)
        k3 = h * myfunc(x[i] + k2 / 2., t[i] + h / 2.)
        k4 = h * myfunc(x[i] + k3, t[i] + h)
        x[i + 1] = x[i] + (k1 + 2. * k2 + 2. * k3 + k4) / 6.
    plt.plot(t, x, 'g-')

    if scenario == 2:  # Here analytical solution is possible write it in the plot instruction!
        # 解析写到清洁说明里面
        plt.plot(t, (t + 1) ** 2 - 0.5 * np.exp(t), 'k.')  # ???

    plt.show()


# Function 1: exponential decay + periodic input (which stops after t=200)
# 函数1:指数衰减+周期性输入(t=200后停止)
def f1(x, t):
    tau = 10
    period = 30
    k = 1
    if t < 200:
        dxdt = -1. / tau * x + np.sin(2 * np.pi * t / period);
    else:
        dxdt = -1. / (k * tau) * x;
    return dxdt


# Some function! 
def f2(x, t):
    dxdt = x - t ** 2 + 1
    return dxdt


if __name__ == "__main__":
    seminar8(1)

 

第一个场景使用函数1，它大致对应于一个由周期性输入供电的CTRNN神经元。周期输入的频率可以通过修改。指数衰减的强度可以通过修改。在时间200处，周期性输入停止，并且可以使用进一步调谐衰变速度。

第二个场景使用函数2。这个函数可以解析积分。考虑了四种积分格式：Euler，Runge-Kutta二阶，中点求积，Runge-Kutta四阶

我们先来研究函数1：

先将main（）函数中的入口seminar8()的参数填写为1：

分别取不同的h值（0.1、1、2）来进行函数运算，导出结果和曲线图：

 

 

 

 

 

 

综上：指数衰减+周期性输入（t=200之后开始停止）可以看到随着h的值不断增长，可以方法之间的偏差也越来越大，当h=2时，四种方法之间的偏差已经非常明显。在准确性方面，RK4应该是最优秀的。

现在我们再来继续研究函数2：

此时只需要将最后一行代码的输入改成2，即可。

我们分别对h取0.1、1、2，运行代码，绘出图像：

 

 

综上：可以看到随着h的值不断增长，各种方法之间的偏差也越来越大。当很轴向右取到5时，各个曲线之间开发发生分离。但无论是那一种h值情况下，均是RK4方法

的变化率最大，可见其优越性。

总结：函数2的阶数更高，即scenario2更为复杂，所以RK4与其他方法的偏差就会越大。

h的值取得越大，那么在固定长度的区间上所取的点就会越多，步长也就越小。有利于数值的准确度。当h很小时，步长变大，就会导致函数在区间内的数值变化很快，出现很大的误差。

标签：plot,myfunc,12,函数,17,arange,plt,2020,np
来源： https://blog.csdn.net/cheng1215peng/article/details/111310385