Incorrect solution of 3 by 3 linear ODE systems

[normalhuman]

dsolve gives wrong solutions of the following 3 by 3 homogeneous linear system. Specifically, the first equation is the one not being satisfied by the solution. If the equations are reordered, the solution changes so that again, 2nd and 3rd equations are satisfied but the 1st is not. From Stack Overflow

k2, k3, t = symbols('k2 k3 t')
x, y, z = symbols('x y z', cls=Function)
eq1 = (Eq(Derivative(x(t),t), k3*y(t)), Eq(Derivative(y(t),t), -(k3+k2)*y(t)), Eq(Derivative(z(t),t), k2*y(t)))
sol1 = dsolve(eq1)
print([eqn.subs({e.lhs: e.rhs for e in sol1}).doit().simplify() for eqn in eq1])
eq2 = (Eq(Derivative(z(t),t), k2*y(t)), Eq(Derivative(x(t),t), k3*y(t)), Eq(Derivative(y(t),t), -(k3+k2)*y(t)))
sol2 = dsolve(eq2)
print([eqn.subs({e.lhs: e.rhs for e in sol2}).doit().simplify() for eqn in eq2])

Expected output: [True, True, True] for both solutions. Actual output from current master:

[Eq(-C3*k3*(k2 + k3)*exp(-t*(k2 + k3))/k2, (C2*k2*exp(t*(k2 + k3)) - C3*k3*(k2 + k3))*exp(-t*(k2 + k3))/k2), True, True]

and

[Eq(C3*k2*exp(-t*(k2 + k3)), C2 + C3*k2*exp(-k2*t - k3*t)), True, True]

[vishalg2235]

I have debugged it and found that in this routine M and evector of both eq are following:

# for first order of equation
(Pdb) M
Matrix([
[0,       k3, 0],
[0, -k2 - k3, 0],
[0,       k2, 0]])
(Pdb) evector
[(0, 2, [Matrix([
[1],
[0],
[0]]), Matrix([
[0],
[0],
[1]])]), (-k2 - k3, 1, [Matrix([
[        k3/k2],
[-(k2 + k3)/k2],
[            1]])])]

# for second order of equation
(Pdb) M
Matrix([
[0, 0,       k2],
[0, 0,       k3],
[0, 0, -k2 - k3]])
(Pdb) evector
[(0, 2, [Matrix([
[1],
[0],
[0]]), Matrix([
[0],
[1],
[0]])]), (-k2 - k3, 1, [Matrix([
[-k2/(k2 + k3)],
[-k3/(k2 + k3)],
[            1]])])]

M vector in both case are different which makes their eigenvectors different and I think it is root cause of wrong result.

@normalhuman which algorithm is used here to solve ODE?

[normalhuman]

A wrong one. The method _linear_neq_order1_type1 assumes that when an eigenvalue is repeated, there is one Jordan block for it (i.e., there is one eigenvector, and the rest are generalized eigenvectors). So it returns incorrect solutions for any matrix that has an eigenspace of dimension greater than 1.

I don't see a way to correct it either, while keeping the overall structure of the method. When an eigenspace is two-dimensional, one can't simply grab an eigenvector v and solve (A-lambda)*w = v to find a generalized eigenvector w.

I think the method needs to be replaced. Two options I considered are

1. Implementing Putzer's algorithm
2. Using the already-implemented matrix exponential exp(t*A) where A is the coefficient matrix.

I don't know which one is better. Putzer does not rely on our ability to compute Jordan blocks; it only needs a flat list of n eigenvalues of A (possibly repeated). We need PR #14323 just to be able to get that list. On the other hand, there is repeated integration (formula 7 in the paper).

Matrix exponential should be much easier to implement, but at present exp(t*A) has complex numbers even when A is real, so the solutions we get appear to be complex. I don't know if this is something to fix in the ODE module, or in the matrix exponential itself. A simple example:

>>> exp(Matrix([[1, -1], [1, 1]]))
Matrix([
[     exp(1 - I)/2 + exp(1 + I)/2, -I*exp(1 - I)/2 + I*exp(1 + I)/2],
[-I*exp(1 + I)/2 + I*exp(1 - I)/2,      exp(1 - I)/2 + exp(1 + I)/2]])

All the entries are actually real, because exp(A) is real when A is real. But this isn't obvious. What exp(A) should be is:

Matrix([
[E*cos(1), -E*sin(1)],
[E*sin(1),  E*cos(1)]]

[normalhuman]

How about this: at the end of the MatrixBase.exp method, check if all entries of self are real, and if they are, return re(ret) instead of ret. Will this break anything?

This change should be made (if it should be made) separately from the ODE issue, because it's not about the ODE, but about better form of exp(Matrix).

[vishalg2235]

I will try to implement putzer's algorithm for solving 'type6' ODE's.

I'am opening new issue for matrix exponential.

[vishalg2235]

What exp(A) should be is:

Matrix([
[E*cos(1), -E*sin(1)],
[E*sin(1),  E*cos(1)]]

What is E here?

[normalhuman]

In SymPy, E denotes Euler's number.

[vishalg2235]

@normalhuman A doubt after computing exp(A*t). As x = exp(A*t), x = Matrix([x(t), y(t), z(t)]) and exp(A*t) will be n*n matrix. Then how to get x(t), y(t) and z(t). Please correct me if I'am wrong somewhere.

[normalhuman]

A search for "matrix exponential, ODE system" should bring up some sources to learn from, such as this and this.

[vishalg2235]

Ok! I should have searched it before.
I think docstring also should be changed ?

[normalhuman]

If the method is rewritten, the docstring should be rewritten too, to match the new method.

[jksuom]

Matrix exponential should be much easier to implement, but at present exp(t*A) has complex numbers even when A is real, so the solutions we get appear to be complex.

Now that was fixed. Maybe this method should also be implemented.

[skirpichev]

Long time ago this was fixed by @cbm755, see #9235.