用 python 解决 ode 得到错误的解决方案

Question

我想解决以下颂歌

K T + C T' = Q

给定的示例数据是我下面的代码

import numpy as np
import scipy as sp
# Solve the following ODE
# K*T + C*T' = Q
# T' = C^-1 ( Q - K * T )
T_start=sp.array([ 151.26,  132.18,  131.64,  146.55,  147.87,  137.87])

K = sp.array([[-0.01761969,  0.02704873,  0.00572222,  0.        ,  0.        ,
         0.        ],
       [ 0.02704873, -0.03546941,  0.        ,  0.        ,  0.00513177,
         0.        ],
       [ 0.00572222,  0.        ,  0.03001858, -0.04752982,  0.        ,
         0.02030505],
       [ 0.        ,  0.        , -0.04752982,  0.0444405 ,  0.00308932,
         0.        ],
       [ 0.        ,  0.00513177,  0.        ,  0.00308932,  0.02629577,
        -0.01793915],
       [ 0.        ,  0.        ,  0.02030505,  0.        , -0.01793915,
         0.00084506]])
Q = sp.array([ 1.66342077,  0.16187956,  0.65115035, 0.71274755,2.54614269,  0.13680399])

C_invers = sp.array([[ 3.44827586,  0.        ,  0.        ,  0.        ,  0.        ,
        -0.        ],
       [ 0.        ,  1.5625    ,  0.        ,  0.        ,  0.        ,
        -0.        ],
       [ 0.        ,  0.        ,  2.63157895,  0.        ,  0.        ,
        -0.        ],
       [ 0.        ,  0.        ,  0.        ,  2.17391304,  0.        ,
        -0.        ],
       [ 0.        ,  0.        ,  0.        ,  0.        ,  1.63934426,
        -0.        ],
       [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         2.38095238]])

time = np.linspace(0, 20, 10000)

#T_real = sp.array([[ 151.26,  132.18,  131.64,  146.55,  147.87,  137.87]])

def deriv(T, t):
    return sp.dot( C_invers,  Q - np.dot(K, T) )

T_sol = sp.integrate.odeint(deriv, T_start, time)

我知道结果是

sp.array([ 151.26,  132.18,  131.64,  146.55,  147.87,  137.87])

当且仅当我使用它作为 T_start 条件时，该解决方案才是“稳定的”

但如果我改变我的开始条件例如

T_start=sp.array([ 0,  0,  0,  0,  0,  0])

它不会收敛我得到以下结果：

我的错在哪里？负值对我的系统没有意义：/你能帮助我吗？谢谢 ;）

Answer 1

数组

array([ 151.26,  132.18,  131.64,  146.55,  147.87,  137.87])

是系统的平衡（大约）。您可以通过将方程组的右侧设置为 0 来找到这一点，这会导致Teq = inv(K)*Q：

In [9]: Teq = np.linalg.solve(K, Q)

In [10]: Teq
Out[10]: 
array([ 151.25960795,  132.17972469,  131.6402527 ,  146.55025359,
        147.87025015,  137.87029892])

这就是为什么当您使用这些值作为起点时您的解决方案看起来是稳定的。该解非常接近平衡，因此不会发生太大变化。

然而，从长远来看，解最终将偏离 Teq，因为平衡点不稳定。您的系统T' = inv(C)*(Q - K*T)是线性的T，因此您可以通过计算 系数矩阵的特征值来确定稳定性T。也就是说，写T = inv(C)*Q - inv(C)*K*T. 的系数矩阵T为-inv(C)*K。以下是找到该矩阵的特征值的方法：

In [11]: A = -C_invers.dot(K)

In [12]: np.linalg.eigvals(A)
Out[12]: 
array([-0.2089754 ,  0.12257481, -0.06349952, -0.01489581,  0.00146708,
        0.05878143])

系数矩阵 A 具有三个正特征值。这些对应于随时间呈指数增长的模式。也就是说，均衡是不稳定的，因此您看到的增长是可以预期的。