equation solving problems

Question

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_] := 
Normal[Series[
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

Answer 1

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]  ;
 Off[Solve::ratnz];
 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 *)

Answer 2

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.

LeafCount[vtab]

(* 29998 *)

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*
    w^3), 
   -((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))} *)

Answer 3

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

Infinitycs = Module[{n = ninfphase, c}, 
Clear[c];
veqexp = 
CoefficientList[
 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 (
      v^\[Prime]\[Prime])[r])/
   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?