用python解决ode问题

Question

我想解决以下问题

KT + CT' = 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èu开始条件

但是如果我把我的开始条件改成

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

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

我的错在哪里？负值对我的系统没有意义：/Can you help me？谢谢；）

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有三个正特征值。它们对应的模式将在时间上呈指数增长。也就是说，平衡是不稳定的，所以你看到的增长是可以预期的。你知道吗