Solving ode with python getting wrong solution

Question

i want to solve the following ode

KT + CT' = Q

to given example Data is my code below

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)

i know that the result is

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

the solution is "stable" if and only if i use this as the T_start condition

but if i change my start condition for example to

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

it won't converge im getting the following result:

where is my fault? Negative values make no sense for my system :/ Can you help me? thanks ;)

score 5 · Answer 1 · answered Feb 02 '18 at 03:08

The array

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

is the equilibrium of your system (approximately). You can find this by setting the right-hand side of your system of equations to 0, which leads to 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])

That's why your solution appears to be stable when you use these values for the starting point. The solution is very close to the equilibrium, so it doesn't change much.

Long term, however, the solution will eventually diverge away from Teq, because that equilibrium point is unstable. Your system, T' = inv(C)*(Q - K*T), is linear in T, so you can determine the stability by computing the eigenvalues of the coefficient matrix of T. That is, write T = inv(C)*Q - inv(C)*K*T. The coefficient matrix of T is -inv(C)*K. Here's how you can find the eigenvalues of that matrix:

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])

The coefficent matrix A has three positive eigenvalues. Those correspond to modes that will grow exponentially in time. That is, the equilibrium is unstable, so the growth that you see is to be expected.

So that means there is no way to determine an unstable equilibrium? Thank you :) — Search898, Feb 02 '18 at 03:23
Is it possible to use additional Information to determine this? For example using an apriori information for the result that the solution should be between 130 and 160? — Search898, Feb 02 '18 at 03:29
For my Problem is `C_invers` arbitrary choosable, but it should be a diagonal matrix. So maybe there is a way to determin `C_inverse` so that `A=-C_invers.dot(K)` has only negative eigenvalues? Thanks guys :) — Search898, Feb 02 '18 at 09:45

Solving ode with python getting wrong solution

1 Answers1