引言

微分方程是一种方程中带有微分项的特殊方程,在不借助计算机的情况下,计算微分方程,人们通常会首先将方程进行分类,确定属于哪一类微分方程,如齐次线性微分方程或者其他。当确定微分方程的种类之后,人们会根据这一类方程的特有解法进行求解微分方程。将会得到y关于x的函数,然后带入x之后,可以得到y的值。以下为百度百科中一阶齐次线性微分方程的通解:

在某些特殊领域,比如机器人控制中,当我们需要求解导纳控制的期望位置xd时,就会需要用到微分方程,这时肯定不可能让你手算结果,因此如何用计算机得到微分方程的解,是控制中必要的一环,也是本节讲述的内容。

通解

在解微分方程的时候,我们首先要定义原本的公式,这里推荐使用dsolve函数:

    # 初始化
    t = symbols('t')
    M = symbols('M')
    B = symbols('B')
    K = symbols('K')
    fext = symbols('fext')
    theta = symbols('f', cls=Function)
    # 表达式
    expr= Eq(M*theta(t).diff(t, t) + B*theta(t).diff(t) + K*theta(t), fext)
    # 求解微分方程
    r1 = dsolve(expr)

如上所示,以导纳控制为例。我们首先定义微分方程的公式,然后使用求解器求得微分的解:

得到:

C1*exp(t*(-B - sqrt(B**2 - 4*K*M))/(2*M)) + C2*exp(t*(-B + sqrt(B**2 - 4*K*M))/(2*M)) + fext/K

强制输出为latex格式后,得到:

可以看到,这时我们得到的解是一个通解。而作为微分方程来说,我们总可以根据实际情况定义初始值,在这里我们定义初始值为:

f0=0  f0'=0

有两种方式达到求解C的目的,即如何求特解。

特解

1、在特解的求解时,我们可以根据上面定义的初始值,利用solve函数,通过求解方程组的方式,得到C1,C2,在得到C时,也可以将参数指定。

    M,B,K = symbols('M,B,K')
    fext = symbols('fext')
    c1,c2 = symbols('c1,c2')
    solution = solve([(c1+c2+fext/K).subs({K:1,fext:10}),(c1*(-B - sqrt(B**2 - 4*K*M))/(2*M)+c2*(-B + sqrt(B**2 - 4*K*M))/(2*M)).subs({M:2,B:1,K:1})],[c1,c2])
    print(solution[c1],solution[c2])

2、当前还有另一种求解的方式:

那就是先通过f0=0求解C1,得到C1关于C2的方程;

然后将C1带入f0’=0求出的C2的公式中,得到最终的C2的值。

    C_eq_1 = Eq(r1.rhs.subs(t,0),0)
    C_sol_1 = solve(C_eq_1)
    # print(C_sol_1[0])
    C_eq_2 = Eq(0,diff(r1.rhs,t).subs(t,0).subs(C_sol_1[0]))
    C_sol_2 = solve(C_eq_2)

计算

在计算时,可以直接将C1,和C2的值带入求解:

r1.subs(C_sol_1[0]).subs(C_sol_2[0]).subs({M:1,B:3,K:1,fext:10}).rhs.subs(t,0.5).evalf()

May the force be with you!