Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I have some equation:

$$ veq=-2-lr-l^2r+2(r+ir^3\omega) v' + (-2+r)r^2v'^2 + (-2+r) r^2 v''==0 $$

or in Mathematica form:

-2 - l r - l^2 r + 
2 (r + I r^3 \[Omega]) Derivative[1][v][r] + (-2 + r) r^2 Derivative[
1][v][r]^2 + (-2 + r) r^2 (v^\[Prime]\[Prime])[r]

then generate veqexp with the code

veqexp[n_] := 
veq /. {v[r_] :> Sum[c[i]/r^i, {i, 1, n}], 
  v'[r_] :> Sum[-i c[i]/r^(i + 1), {i, 1, n}], 
  v''[r_] :> 
   Sum[i (i + 1) c[i]/r^(i + 2), {i, 1, n}]}, {r, \[Infinity], 
 n - 2}]];

I sub into this an ansatz for solution $v=\sum^n c_i r^i$. Calling this new equation in terms of the $c_i$ 'veqexp' (i.e. veq expanded). The following chunk of code can solve this for my coefficients: (note I want to solve around infinity)

vcoeffs[nn_] := 
Block[{}, Clear[c]; 
Do[c[i] = c[i] /. Solve[
     Limit[veqexp[nn] r^(i - 2), 
       r -> \[Infinity]] == 0, c[i]][[1]];, {i, 1, nn}]] ;

and out correctly pop the $c_i$ coefficients in terms of $(\omega,\ell)$ parameters of my 'veq', so far so good. The problem is that these functions of $(\omega,\ell)$rapidly grow in size until by about the 30th Mathematica's memory gives out and dies. However if I try to calculate for a given $(\omega,\ell)=(0.1,1)$ so that each c[i] is just a number, I also hit problems to do with recursion that I have no idea about. For example

vcoeffs[nn_] := 
Block[{}, Clear[c]; 
Do[c[i] = 
  c[i] /. Solve[
     Limit[(veqexp[nn] /. {\[Omega] -> 0.1, l -> 1}) r^(i - 2), 
       r -> \[Infinity]] == 0, c[i]][[1]];, {i, 1, nn}]] ;

generates errors that read

 $RecursionLimit::reclim: Recursion depth of 256 exceeded. >>

when I would have expected just a series of $nn$ complex coefficients which would take no memory for Mathematica to store.

Can anyone suggest a way around this? thanks

share|improve this question
Is ir in your veq equation a single variable, or i*r and, in the latter case, is i then Mathematica's I? I'd prefer a non-ambiguous and above all copyable Mathematica definition instead of a pretty printing LaTeX one in this case. – Sjoerd C. de Vries Sep 27 '12 at 17:52
If veqexp is an equation as you say, what does it do in ` Limit[veqexp[nn] r^(i - 2), r -> [Infinity]]`? How do you take the limit of an equation? – Sjoerd C. de Vries Sep 27 '12 at 18:30
$i$ is the Mathematica $I$ yes. $veqexp[nn]$ gets hit by the $r^{i-2}$ then it turns out that the only terms involving c[i] don't have negative powers of $r$ as factors, so in the limit $r\to \infty$ we get a single variable equation for c[i] to feed into Solve[...]. As I say all this works as I want, but if I try to make it specific parameters $(\omega,\ell)$ then recursion issues. Alternatively the coefficients c[i] get too big to store in memory when $i>20$. – fpghost Sep 27 '12 at 18:46
Please add the full Mathematica definitions of veq and veqexp. – Sjoerd C. de Vries Sep 27 '12 at 18:49
I'll second Sjoerd's request. It makes no sense to have every respondent rewriting the HTML form of veq into Mathematica. – Daniel Lichtblau Sep 27 '12 at 23:24
up vote 1 down vote accepted

I understand that you want to use your ansatz, then choose the coefficients in order to kill the highest n powers so that your solution is well behaved at r=Infinity. The main problem I found is that I don't know how to use CoefficientList[poly] when the polynomial has negative powers. The following is my approach : I use Together because it will put the result of the substitution in the form 1/r^m \[Times] (polynomial with positive powers only).

This is your equation :

expr[v_] := -2 - l r - l^2 r + 2 (r + I r^3 \[Omega]) D[v, r] + (-2 + r) r^2
 D[v, r]^2 + (-2 + r) r^2 D[v, {r, 2}]

and this is my solution. I use numQ = 1 for an exact solution, numQ = 1. to use inexact coefficients.

mysol[n_, numQ_] := First@Module[{cVec, rVec, exprAnsatz, coefficientRules},
  cVec = Subscript[c, #] & /@ Range[1, n];
  rVec = r^-# & /@ Range[1, n];
  exprAnsatz = Together[(expr[Dot[cVec, rVec]] // Expand )] ;
  coefficientRules = CoefficientRules[Together[(expr[Dot[cVec, rVec]] // Expand // Simplify)] [[2]], 
 r]  ;
 Solve[Thread[numQ coefficientRules[[1 ;; n, 2]] == 0], cVec]

This is keeping \[Omega], l symbolic (it's faster with numbers) :

mysol[40, 1]; // AbsoluteTiming
mysol[40, 1.]; // AbsoluteTiming   
(* {15.207711, Null} *)

Before checking against your solution I redefined your function

myvcoeffs[nn_] := Block[{}, 
 Clear[c];Table[c[i] = c[i] /. Solve[Limit[veqexp[nn] r^(i - 2), r -> \[Infinity]] == 0, c[i]][[1]], {i, 1, nn}]]

(myvcoeffs[10] == mysol[10][[All, 2]]) // Simplify
(* True *)
share|improve this answer
thanks for the reply. Could you possibly talk me through what is going on in your mysol[..], it seems to work nicely but I don't understand most of commands unfortunately.. – fpghost Sep 29 '12 at 10:08
Please have a look at the documentation, then ask specific questions. You can run the code inside Module one line a time and see what it does (though I tried to explain how it works). – b.gatessucks Sep 29 '12 at 10:12

I'd form a Laurent series as an unknown "function" v[] for which we need to find derivatives. Isolate powers of r and set each to zero.

n = 30;
v = Series[vv[r], {r, inf, n}] /. inf -> Infinity;
vpoly = -2 - l*r - l^2*r + 
   2*(r + I*r^3*w)*D[v, r] + (-2 + r)*r^2*D[v, r]^2 + (-2 + r)*r^2*
    D[v, {r, 2}];
vpolyn = Normal[vpoly];
vpolyn = r^-Exponent[vpolyn, r, Min]*vpolyn;

coeffs = Reverse[CoefficientList[vpolyn, r]];
vars = Select[Variables[coeffs], ! FreeQ[#, vv] &];

We now solve for this system. There are two key observations.

(1) The jth coefficient in the list is linear in the jth derivative, and only uses other derivatives of lower order. That is, the system is in some sense triangulated. As Solve will not necessarily grasp this we may do better to solve iteratively.

(2) If, after each derivative is found, we use Together to simplify it prior to back substitution, we might keep coefficient swell to a reasonable level.

The simple loop below avails itself of these observations.

Timing[vtab =
   Table[var = Derivative[j][vv][Infinity];
    soln = Solve[coeffs[[j]] == 0, var][[1]];
    soln = soln /. Rule[aa_, bb_] :> Rule[aa, Together[bb]];
    coeffs = coeffs /. soln;
    var /. soln
    , {j, Length[coeffs]}];]

(* {3.38, Null} *)

Total size is now manageable.


(* 29998 *)


Here are the solutions for the first few derivatives.

vtab[[1 ;; 5]]

(* {(I*(l + l^2))/(2*w), (l + l^2 + 2*I*w)/(2*w^2), 
   (-6*I*l - 5*I*l^2 + 2*I*l^3 + I*l^4 + 12*w - 12*l*w - 12*l^2*w)/(4*
   -((3*(6*l + 4*l^2 - 4*l^3 - 2*l^4 + 12*I*w - 22*I*l*w - 
        21*I*l^2*w + 
              2*I*l^3*w + I*l^4*w + 16*w^2))/(2*w^4)), 
   (1/(2*w^5))*(3*(60*I*l + 32*I*l^2 - 55*I*l^3 - 25*I*l^4 + 3*I*l^5 + 
           I*l^6 - 120*w + 300*l*w + 260*l^2*w - 80*l^3*w - 
      40*l^4*w + 288*I*w^2 - 
           136*I*l*w^2 - 136*I*l^2*w^2))} *)
share|improve this answer

I also now know of another way to find the coefficients:

Infinitycs = Module[{n = ninfphase, c}, 
veqexp = 
 Series[(-2 - l r - l^2 r + 
    2 (r + I r^3 \[Omega]) Derivative[1][v][
      r] + (-2 + r) r^2 Derivative[1][v][r]^2 + (-2 + r) r^2 (
   r /. {v'[r_] :> Sum[-i c[i]/r^(i + 1), {i, 1, n}], 
    v''[r_] :> 
     Sum[i (i + 1) c[i]/r^(i + 2), {i, 1, n}]}, {r, \[Infinity], 
   n - 1}], r^-1];
Do[c[i] = c[i] /. Simplify[Solve[veqexp[[i]] == 0, c[i]][[1]]]; 
Print, {i, 1, n}] ;
Table[c[i], {i, 1, n}]];

What is it about the functioning of 'Module' that means these coefficients $c_i$ now don't eat the memory? I thought Module just localized the variable, but would still be having to store these $n$ large coefficients in memory?

share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.