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).
When I ask, say, for the second homotopy derivative of y^2, I should get y[1]^2 + 2y[0]y[2] (remember, it's the derivative at q=0, so all other terms vanish). The problem is that the only way I can think of to get these kind of results is to predefine a finite series, which isn't very dynamic in the long run. You see, these calculations can become monstrous within just a few terms (after all, every previous solution contributes to the next one). It would be ideal if I could find the first few terms, see if that gets me close enough, and then, if necessary, just add the next few terms without having to recalculate the whole process. The problem is that when I begin with a truncated series, I may find that I need more terms than I can obtain accurately, which means I'd have to recalculate the whole problem.
So how about those sums? Is there a way to express an unknown infinite sum in an algebraically useful way in Maxima?
Maybe think in these terms: Let y and z be power series in q.
D(m, y*z) := block(
step one: calculate the resultant series out to order m
step two: choose the coefficient of q^m and return it (simplified)
);
Damn, that's actually a pretty good start...any thoughts at all would be greatly appreciated.
-- Casey