2020-12-17

2023-10-11 13:14:22

复杂数学系统 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）来进行函数运算，导出结果和曲线图：