Is there a numerical ODE solver facility in Maxima

Richard Fateman wrote:

> Using a computer algebra system there is another approach, which is to
> solve for the taylor series coefficients in f(x).
> There are various packages to do this, though I don't know of one that
> is set up in Maxima.
> 
> A trivial setup (but expensive to run) is to define f, using deftaylor,
> with symbolic coefficients a0, a1, ...
> and then expand f to some order, substitute in the DE, and solve for a0,
> a1, etc.
> 
> RJF

The colnew program solves the differential equation by approximating the
solution by splines with indeterminate coefficients, and then imposing the
differential equation at so-called "collocation points", which are similar
to Gauss points in numerical integration. This is proved to give better
convergence. Taylor expansions are most useful at singular points of a
solution, when they exist. In this case it is a criterion for
integrability, known as the Kowalewski-Painlev? analysis. It is trying to
plug such a singular expansion in the top equation that Sophie Kowalewski
found her celebrated integrability case. This sort of thing is very
appropriate for Maxima computations. For example i have done the exercise
here (please be indulgent i don't know all Maxima tricks by far, but at
least this worked for me).  

-- 
Michel Talon

niobe% cat kowa.max
/* The aim is to show that one can find Laurent developments for the
dynamical
 * variables p, q, r, etc. of the Kowalevski top, in the cases where the
 * determinant of the recursive system has sufficiently many integer roots
*/

tt: elapsed_real_time ();

/* The matrix of the recursive system */

Mat: matrix(
[(m-1)*A, -A1*r, -A1*q, 0, z, -y],
[-B1*r, (m-1)*B, -B1*p, -z, 0, x],
[-C1*q, -C1*p, (m-1)*C, y, -x, 0],
[0, h, -g, m-2, -r, q],
[-h, 0, f, r, m-2, -p],
[g, -f, 0, -q, p, m-2] );

Mat: subst([A1 = B-C, B1 = C-A, C1 = A-B], %)$

/* to get %i^2 = -1 */
algebraic: true$

/* the first Kowalevski solution */
Mat1: subst([B=A, A=2*C, y=0, z=0, h=0, q=%i*p, r=2*%i, f=-2*C/x,
                g=-2*%i*C/x] , Mat)$

Det1: factor(determinant(Mat1));


/* The second Kowalevski solution */
Mat2: subst([B=A, A=2*C, y=0, z=0, q=2*%i, r=0, f=-4*C/x,
                g=0, h=4*%i*C/x, p=0], Mat)$

Det2: factor(determinant(Mat2));

/* Now we need to solve the recursive equations for the Laurent series */

/* this is a column vector, the unknowns */
v: matrix([p[m]], [q[m]], [r[m]], [f[m]], [g[m]], [h[m]])$
/* note that v[1] = [p[m]] and part(v[1],1) = p[m] */

/* the right hand side of the recursion relations, starting with 0 */
w: [P[m], Q[m], R[m], F[m], G[m], H[m]]$
P[1]: 0$ Q[1]:0$ R[1]: 0$ F[1]: 0$ G[1]: 0$ H[1]: 0$

vars: makelist(part(v[j],1), j, 1, 6)$

print ( "The first Kowalevski case.")$

/* We solve the first Kowalevski case */

eqns: makelist(part((Mat1.v)[j],1) = w[j], j, 1, 6)$

/* solve for m=1 */
solve(subst(m=1, eqns), subst(m=1, vars))$
/* assign the solution to vars, :: evaluates right and left */
subst(m = 1, vars) :: map(rhs, %[1])$
print("First iteration. Vars:", %)$

/* now for m=2,3,4 */
for mm from 2 thru 4 do (
P[mm]: C*sum(q[j]*r[mm-j], j, 1, mm-1),
Q[mm]: -C*sum(p[j]*r[mm-j], j, 1, mm-1),
R[mm]: 0,
F[mm]: sum(r[j]*g[mm-j] - q[j]*h[mm-j], j, 1, mm-1),
G[mm]: sum(p[j]*h[mm-j] - r[j]*f[mm-j], j, 1, mm-1),
H[mm]: sum(q[j]*f[mm-j] - p[j]*g[mm-j], j, 1, mm-1),

solve(ev(subst(m=mm, eqns)), subst(m=mm, vars)),
/* in a block, % has to be replaced by %% */
subst(m=mm, vars) :: map(rhs, %%[1]) )$

/* Display the Laurent solution */

print ( "Laurent solution depending on 5 parameters,
p, %r1, %r2, %r3, %r4")$

[p[0], q[0], r[0], f[0], g[0], h[0]]: [p, %i*p, 2*%i, -2*C/x, -2*%i*C/x, 0]$

for j in [p, q, r] do block([z],
z: fullmap(ev, j[0]/t + j[1] + t*j[2] + t^2*j[3] + t^3*j[4]),
print(j, "(t) = ", expandwrt_factored(z, t, p)),
define(ev(j(t)), z) )$ /* this is for the check below */

for j in [f, g, h] do block([z],
z: fullmap(ev, j[0]/t^2 + j[1]/t + j[2] + t*j[3] + t^2*j[4]),
print(j, "(t) = ", expandwrt_factored(z, t, p)),
define(ev(j(t)), z) )$

print ( "Checking integrals of motion.")$

/* Check that weight has constant magnitude, this is central for Poisson  */
taylor(f(t)^2 + g(t)^2 + h(t)^2, t, 0, 0);
/* Notice no terms in 1/t^4, 1/t^3, 1/t^2, 1/t !!! */

/* Check another integral of motion, in the centre of Poisson algebra  */
taylor(2*C*(p(t)*f(t)+q(t)*g(t))+C*r(t)*h(t), t, 0, 1);
/* Notice no terms in 1/t^3, 1/t^2, 1/t, t !!! */


/* The phase space is thus or real dimension 4 and not 6. In the integrable
 * case there are 2 integrals of motion
*/
/* Check the energy integral of motion */
taylor( C*r(t)^2+2*C*((p(t))^2+(q(t))^2) -2*x*f(t), t, 0, 2);
/* Notice no terms in 1/t^2, 1/t, t, t^2 !!! */

/* Check the new Kowalevski integral of motion */
taylor(
((p(t)+%i*q(t))^2+x*(f(t)+%i*g(t))/C)*((p(t)-%i*q(t))^2+x*(f(t)-%i*g(t))/C),
t, 0, 0);
/* Notice no terms in  1/t^4, 1/t^3, 1/t^2, 1/t !!! */

ttt: elapsed_real_time ()$ ttt -tt;

print ( "Now the second Kowalewski case." )$

/* We solve the second Kowalevski case */

/* Reset variables */

kill(allbut (Mat2, w, ttt))$
v: matrix([p[m]], [q[m]], [r[m]], [f[m]], [g[m]], [h[m]])$
w: [P[m], Q[m], R[m], F[m], G[m], H[m]]$
P[1]: 0$ Q[1]:0$ R[1]: 0$ F[1]: 0$ G[1]: 0$ H[1]: 0$


vars: makelist(part(v[j],1), j, 1, 6)$

eqns: makelist(part((Mat2.v)[j],1) = w[j], j, 1, 6)$

/* solve for m=1 */
solve(subst(m=1, eqns), subst(m=1, vars))$
/* assign the solution to vars, :: evaluates right and left */
subst(m = 1, vars) :: map(rhs, %[1])$
print("First iteration. Vars:", %)$

/* now for m=2,3,4 . Note the determinant vanishes to second order at m=2,
 * better find 2 solutions here.
*/
[mm]: -C*sum(p[j]*r[mm-j], j, 1, mm-1),
R[mm]: 0,
F[mm]: sum(r[j]*g[mm-j] - q[j]*h[mm-j], j, 1, mm-1),
G[mm]: sum(p[j]*h[mm-j] - r[j]*f[mm-j], j, 1, mm-1),
H[mm]: sum(q[j]*f[mm-j] - p[j]*g[mm-j], j, 1, mm-1),
solve(ev(subst(m=mm, eqns)), subst(m=mm, vars)),
subst(m=mm, vars) :: map(rhs, %%[1]) )$



/* Display the Laurent solution */

print ( "Laurent solution depending on 4 parameters,  %r5, %r6, %r7, %r8")$

[p[0], q[0], r[0], f[0], g[0], h[0]]: [0, 2*%i, 0, -4*C/x, 0, 4*%i*C/x]$

for j in [p, q, r] do block([z],
z: fullmap(ev, j[0]/t + j[1] + t*j[2] + t^2*j[3] + t^3*j[4]),
print(j, "(t) = ", expandwrt_factored(z, t, p)),
define(ev(j(t)), z) )$ /* this is for the check below */

for j in [f, g, h] do block([z],
z: fullmap(ev, j[0]/t^2 + j[1]/t + j[2] + t*j[3] + t^2*j[4]),
print(j, "(t) = ", expandwrt_factored(z, t, p)),
define(ev(j(t)), z) )$

print ( "Checking integrals of motion.")$

taylor(f(t)^2 + g(t)^2 + h(t)^2, t, 0, 0);

taylor(2*C*(p(t)*f(t)+q(t)*g(t))+C*r(t)*h(t), t, 0, 1);

taylor( C*r(t)^2+2*C*((p(t))^2+(q(t))^2) -2*x*f(t), t, 0, 2);

taylor(
((p(t)+%i*q(t))^2+x*(f(t)+%i*g(t))/C)*((p(t)-%i*q(t))^2+x*(f(t)-%i*g(t))/C),
t, 0, 0);

elapsed_real_time () - ttt;