自治微分方程

自治微分方程指的是$$\frac{dy}{dt}=f(y(t))$$的右边是只包含$$y$$的函数的微分方程。例如前面章节出现的人口增长模型。

任何自治微分方程都是变量课分离的，但是积分可能会很困难。

给定一个自治微分方程的解，将其沿着水平方向横移，我们可以获得很多组不同的解。

例如之前也出现过的$$\frac{dy}{dt}=y(1-y)$$

对于自治微分方程，斜率场中的冗余信息太多，用相线（phase line）来简洁地表示这些信息。

微分方程$$\frac{dy}{dt}=y(1-y)$$的平衡解为$$y(t)=0$$和$$y(t)=1$$

1. 在$$y>1$$的区域，$$\frac{dy}{dt}<0$$，即形为$$y(t_0)=y_0, y_0="">1$$的初值问题，随着$$t$$的增加，函数值终会收敛于$$y(t)=1$$
2. 在$$y<0$$的区域，$$\frac{dy}{dt}<0$$，即形为$$y(t_0)=y_0, y_0 <0$$的初值问题，随着$$t$$的增加，函数值终会收敛于$$y(t)=-\infty$$
3. 在$$00$$,即形为$$y(t_0)=y_0, 0<y_0 <1$$的初值问题，，随着$$t$$的增加，函数值终会收敛于$$y(t)=1$$

用一根直线表示$$y(t)$$的值，用点表示平衡解，用箭头表示各个区域的初值问题，随着$$t$$的增加，$$y(t)$$值的走向：

    import sympy
    from sympy.abc import t
    from sympy import Function, Derivative, dsolve, Eq, solve

    def phaseLine(formula, points = [], domain = np.linspace(-1,1), domainFix = False):
        if points != []:
            solutions = points
        else:
            solutions = solve(Eq(formula,0),y(t))
        try:
            solutions.sort()
            ran = solutions[-1]-solutions[0]
            if len(solutions) == 1:
                intervals = [solutions[0]-0.5, solutions[0], solutions[0]+0.5]
            else:
                intervals = [float(solutions[0]-0.25*ran)]+solutions+[float(solutions[-1]+0.25*ran)]

            if domainFix:
                ydomain = domain
            elif len(solutions) > 1:
                ydomain = np.linspace(float(solutions[0]-0.25*ran), float(solutions[-1]+0.25*ran))
            elif len(solutions) == 1:
                ydomain = np.linspace(float(solutions[0]-1),float(solutions[0]+1))
        except:
            if domainFix:
                ydomain = domain
            else:
                ydomain = np.linspace(-1,1)

            if 0 in solutions:
                signVals = [ydomain[12],ydomain[38]]
            else:
                signVals = [0]
            solutions = []

        fig = plt.figure(num=1, figsize=(1,4))
        plt.plot([0 for dummy in ydomain],ydomain, color = 'black')
        plt.plot([0 for dummy in solutions], solutions, 'or')
        plt.axis(xmin = -0.5, xmax = 0.5)

        if solutions != []:
            for i in range(len(intervals)-1):
                midpoint = (intervals[i]+intervals[i+1])/2.0
                intervalSign = sign(formula.evalf(subs={'y(t)':midpoint}))
                if intervalSign == -1:
                    plt.text(0, midpoint, u'\u2193',fontname='STIXGeneral', color = 'blue', size=30, va='center', ha='center', clip_on=True)
                elif intervalSign == 1:
                    plt.text(0, midpoint, u'\u2191',fontname='STIXGeneral', color = 'blue', size=30, va='center', ha='center', clip_on=True)
        else:
            for val in signVals:
                intervalSign = sign(formula.evalf(subs={'y(t)':val}))
                if intervalSign == -1:
                    plt.text(0, val, u'\u2193',fontname='STIXGeneral', color = 'blue', size=30, va='center', ha='center', clip_on=True)
                elif intervalSign == 1:
                    plt.text(0, val, u'\u2191',fontname='STIXGeneral', color = 'blue', size=30, va='center', ha='center', clip_on=True)
        return fig

    y = Function('y')
    formula = y(t)*(1-y(t))
    fg = phaseLine(formula)
    fg.show()

练习：画出$$\frac{dy}{dt}=y^2-4y-12$$的相线

    formula = y(t)**2 - 4*y(t) - 12
    fg2 = phaseLine(formula)
    fg2.show()

练习：画出$$\frac{dy}{dt}=y^2cosy$$的相线

    formula = y(t)**2*sympy.cos(y(t))
    fg3 = phaseLine(formula, [-3/2.0*pi,-1.0/2*pi,0,1.0/2*pi,3/2.0*pi])
    fg3.show()

绘制$$\frac{dy}{dt}$$关于$$y$$的函数图，有助于绘制相线。