Casey Coleman wrote:
> I can't work with nth-term formulas because I don't know the series
beforehand. I have to solve for each of the coefficients of q^n, which in this
case are each functions of a single variable x. I do this (in theory) by
setting up an infinite chain of linear recurrence relations. This infinite
chain is obtained by taking the homotopy derivative(s) of the so called
zeroth-order deformation equation:
>
> (1-q)L[y-u[0]] = qhN[y]
>
> Which has the following properties:
> q is the homotopy parameter ranging from 0 to 1
> L is an auxiliary linear operator (independent of q) which I choose
> u[0] is my initial guess solution
> N[y] is the nonlinear operator I'm interested in solving
> y is a homotopy series whose coefficients are to be determined (in such
a way that it solves N[y]=0.
> h is the convergence control parameter, which is basically just an
iteration factor.
>
> The principle task at hand is to compute D[m](N[y]) for any m. Doing so is
easy in principle since only three operations are needed:
> 1.) D[m](y) = y[m], the coefficient of q^m. This is true for any power
series in q.
> 2.) Taking derivatives of quantities like e^y (for strongly nonlinear
problems?)
> 3.) determining the coefficients of series defined by ordinary algebraic
expressions like y^2.
>
> So here is what I'm talking about:
> y = y[0] + y[1]q + y[2]q^2 + y[3]q^3 +...
> y[0] is a guess, so I know what it is. Using the deformation equations,
which are governed primarily by the linear operator I choose and my initial
guess, I will find the y[n], for any n that I choose, recursively (the y[n]'s
are functions of x in this case
If this may be of any help, here is an example of a similar problem, in which
the Kowalevski top equations are "solved" by plugging a Laurent series (the
so-called Kowalevsky-Painlev? criterion) and recursively solving them.
/* 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 */
/* Can be removed by setting: globalsolve:true */
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.
*/
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)),
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;
--
Michel Talon