Mathematica 6 and system of matrix and scalar ODEs in NDSolve: NDSolve::nlnum

*To*: mathgroup at smc.vnet.net*Subject*: [mg83667] Mathematica 6 and system of matrix and scalar ODEs in NDSolve: NDSolve::nlnum*From*: Michael Schmidt <michaelschmidt.hd at googlemail.com>*Date*: Wed, 28 Nov 2007 05:24:10 -0500 (EST)

Hi you all, I encountered a strange behaviour of NDSolve in Mathematica 6.0. In the mathematica package REAP which we developed for Mathematica 5.0, a system of ordinary differential equations has to be solved. Part of it are matrices, part of it are ordinary scalars. Mathematica 6.0, however, will show the error message NDSolve::nlnum , if the variable name of the matrix equation is "Yd". In the case, that it is "Hd", everything is fine. One might argue, that Yd is used as an internal variable in NDSolve, but then it is even more strange that the same happens for "REAP`RGEMSSM`Private`Yd". Everything works fine in Mathematica 5.x . Also, this problem will only show up, if there is a system of ODEs where there are matrices and scalars. Has anyone of you encountered the same error? Does anyone of you know, why the result of NDSolve depends on the name of the parameter? Here is some sample code: NDSolve[{D[y[t], t] == -((3 y[t]^3)/(16 \[Pi]^2)), D[Yd[t], t] == 1/(16 \[Pi]^2) (Yd[t]. Transpose[Conjugate[Yd[t]]].Yd[t] - (16 y[t]^2)/3 Yd[t]), y[Log[20000000000000000]] == 0.6984, Yd[Log[20000000000000000]] == {{0.5, 0, 0}, {0, 0.6, 0}, {0, 0, 0.7417}}}, {y, Yd}, {t, Log[100], Log[2 10^16]}]; results in NDSolve::nlnum: "The function value \ {-0.00647232,{{-0.00744611,0.,0.},{0.,-0.0085173,0.},{0.,0.,-0.\ 00963568}}} is not a list of numbers with dimensions {10} at \ {t,y[t],Yd[t]} = \ {37.5307,0.698424,{{0.500028,0.,0.},{0.,0.600032,0.},{0.,0.,0.741736}}\ }." NDSolve[{D[y[t], t] == -((3 y[t]^3)/(16 \[Pi]^2)), D[Hd[t], t] == 1/(16 \[Pi]^2) (Hd[t]. Transpose[Conjugate[Hd[t]]].Hd[t] - ( 16 y[t]^2)/3 Hd[t]), y[Log[20000000000000000]] == 0.6984, Hd[Log[20000000000000000]] == {{0.5, 0, 0}, {0, 0.6, 0}, {0, 0, 0.7417}}}, {y, Hd}, {t, Log[100], Log[2 10^16]}]; NDSolve[{D[y[t], t] == -((3 y[t]^3)/(16 \[Pi]^2)), D[REAP`RGEMSSM`Private`Yd[t], t] == 1/(16 \[Pi]^2) (REAP`RGEMSSM`Private`Yd[t]. Transpose[ Conjugate[ REAP`RGEMSSM`Private`Yd[t]]].REAP`RGEMSSM`Private`Yd[t] - ( 16 y[t]^2)/3 REAP`RGEMSSM`Private`Yd[t]), y[Log[20000000000000000]] == 0.6984, REAP`RGEMSSM`Private`Yd[ Log[20000000000000000]] == {{0.5, 0, 0}, {0, 0.6, 0}, {0, 0, 0.7417}}}, {y, REAP`RGEMSSM`Private`Yd}, {t, Log[100], Log[2 10^16]}]; Best regards, Michael