#3173 coefmatrix problem

[Daniel Volinski]

The following set of comands breaks Maxima

load(eigen)$
A:matrix([1,1],[4,1]);
v:columnvector([x(t),y(t)]);
L:list_matrix_entries(A.v);
Eq1:'diff(x(t),t)=L[1];
Eq2:'diff(y(t),t)=L[2];
sol:desolve([Eq1,Eq2],[x(t),y(t)])$
define(x(t),rhs(sol[1]));
define(y(t),rhs(sol[2]));
[vals,vecs]:eigenvectors(A);
coefmatrix([x(t),y(t)],[exp(vals[1][1]*t),exp(vals[1][2]*t)]);
Message from maxima's stderr stream: 
*** - Program stack overflow. RESET
SERVER: Lost socket connection ...
Trying to restart Maxima.

build_info(version="5.38.0_dirty",timestamp="2016-04-03 10:09:05",
host="i686-w64-mingw32",lisp_name="CLISP",lisp_version=
"2.49 (2010-07-07) (built on maxima [192.168.43.2])")

[Daniel Volinski]

I forgot a very important line before the last line:

Thanks,
Daniel Volinski

[Daniel Volinski]

The line before the last line is: [vals,vecs]:eigenvectors(A);

[Mike]

According to the manual, the define command evaluates its second argument, which causes a problem here since x(0) and y(0) are found in sol. I suggest using

x(t):=rhs(sol[1]);

and

y(t):=rhs(sol[2]);

instead. Alternatively, the evaluation of x(0) and y(0) could be manually suppressed by defining sol as

sol:subst([x(0)='x(0),y(0)='y(0)],desolve([Eq1,Eq2],[x(t),y(t)]))$

in which case your define commands should work.

[Robert Dodier]

This is a good solution too. I think it's simpler than mine, so it's to be preferred.

[Daniel Volinski]

Hi Mike,

I'll take you suggestion.

Thanks,

Daniel

[Robert Dodier]

Enclose code in four tildes for formatting; also paste in missing line.

[Robert Dodier]

The problem is that you have defined x(t) and y(t) in terms of x(0) and y(0) but the definition doesn't say what to do if t=0. So x(t) calls x(0) and x(0) calls x(0) and there is a stack overflow. (Not a problem in coefmatrix.)

There are a couple of ways to handle it. Here's one way which is to check for t=0 and return a noun expression (i.e. prevent function evaluation) in that case.

(%i2) load(eigen)$
(%i3) A:matrix([1,1],[4,1]);
                     [ 1  1 ]
(%o3)                [      ]
                     [ 4  1 ]
(%i4) v:columnvector([x(t),y(t)]);
                     [ x(t) ]
(%o4)                [      ]
                     [ y(t) ]
(%i5) L:list_matrix_entries(A.v);
(%o5)      [y(t) + x(t), y(t) + 4 x(t)]
(%i6) Eq1:'diff(x(t),t)=L[1];
             d
(%o6)        -- (x(t)) = y(t) + x(t)
             dt
(%i7) Eq2:'diff(y(t),t)=L[2];
            d
(%o7)       -- (y(t)) = y(t) + 4 x(t)
            dt
(%i8) sol:desolve([Eq1,Eq2],[x(t),y(t)])$
(%i9) define(x(t), subst (body = rhs(sol[1]),
                                                   '(if t = 0 then 'x(0) else body)));
(%o9) x(t) := if t = 0 then 'x(0)
                        3 t
      (y(0) + 2 x(0)) %e
 else ---------------------
                4

                     - t
   (y(0) - 2 x(0)) %e
 - ---------------------
             4
(%i10) define(y(t), subst (body = rhs(sol[2]),
                                                      '(if t = 0 then 'y(0) else body)));
(%o10) y(t) := if t = 0 then 'y(0)
                        3 t
      (y(0) + 2 x(0)) %e
 else ---------------------
                2
                     - t
   (y(0) - 2 x(0)) %e
 + ---------------------
             2
(%i11) [vals,vecs]:eigenvectors(A);
(%o11) [[[3, - 1], [1, 1]], 
                           [[[1, 2]], [[1, - 2]]]]
(%i12) coefmatrix([x(t),y(t)],[exp(vals[1][1]*t),exp(vals[1][2]*t)]);
               [ y(0) + 2 x(0)    ]
               [ -------------  0 ]
               [       4          ]
(%o12)         [                  ]
               [ y(0) + 2 x(0)    ]
               [ -------------  0 ]
               [       2          ]

Note that the x(0) and y(0) you see here are noun expressions. You'll have to take that into account if you try to do some operations on them.

(%i13) noundisp:true $
(%i14) [x(t), y(t)];
                            3 t
        ('y(0) + 2 'x(0)) %e
(%o14) [-----------------------
                   4

                       - t
   ('y(0) - 2 'x(0)) %e
 - -----------------------, 
              4
                    3 t                       - t
('y(0) + 2 'x(0)) %e      ('y(0) - 2 'x(0)) %e
----------------------- + -----------------------]
           2                         2

Another way to handle the problem is to use simplification rules and simplify only if t # 0. I'll leave that aside for now.

[Daniel Volinski]

Thank you Mike and Robert for the responses.
Daniel