Solving an ODE in power series

How do I find a series solution to an ODE? I do not mean taking the Taylor series of an exact solution; I want to solve nasty nonlinear differential equations locally via plug and chug. Surely, that is a built-in option in Mathematica.

-
Yes, there's a way. Do you have an example problem? –  Mark McClure May 17 '13 at 3:45
Somewhat surprisingly, the procedure you want isn't built-in, but one could certainly write a routine for this task. –  Ｊ. Ｍ. May 17 '13 at 3:46
@J.M. Not true! It's buried in the Holonomic context –  Mark McClure May 17 '13 at 3:50
@Mark, I see. I had something like this in mind... still, if memory serves, the functions there can only do linear ODEs, not things like $y^\prime=1+y^2$. –  Ｊ. Ｍ. May 17 '13 at 3:55
@J.M. Yes, that's definitely correct. Holonomic functions, in fact, satisfy linear ODEs by definition, which is exactly why I asked for an example, rather than post an answer. –  Mark McClure May 17 '13 at 4:04
show 1 more comment

This method of developing a truncated solution can be done as below. I illustrate with an example that DSolve does not seem much to like.

ode = x''[t] - t*x'[t] + Sin[t] == 0;
initconds = {x'[0] == 1, x[0] == 0};


We create a differentail operator to create this ode.

odeOperator = D[#, {t, 2}] - t*D[#, t] + Sin[t] &;


Now set up our Taylor series as a symbolic expansion using derivatives of x evaluated at the origin. I use an order of 15 but that is something one would probably make as an argument to a function, if automating all this.

xx = Series[x[t], {t, 0, 15}];


Next apply the differential operator and add the initial conditions. Then find a solution that makes all powers of t vanish.

soln = SolveAlways[Join[{odeOperator[xx] == 0}, initconds], t];


Let's look at the Taylor polynomial.

truncatedSol = Normal[xx /. soln[1]]

(* Out[500]= t + t^5/120 + t^7/1260 + (29 t^9)/362880 + ( 13 t^11)/1995840 + (2861 t^13)/6227020800 + (4649 t^15)/163459296000 *)

To assess how accurate it might be I will compare with NDSolve.

approxSol = NDSolve[Join[{ode}, initconds], x[t], {t, 0, 4}];

Plot[{truncatedSol, x[t] /. approxSol[[1]]}, {t, 0, 4}]
`

-
+1 This is a good start. Power series solutions, though, are frequently used to obtain recursion equations for the coefficients (of any solution that might be analytic within a neighborhood of the point of expansion). It would be nice, then, to have a function that outputs these equations (given a differential operator as input), rather than just obtaining an approximate solution with a limited radius of accuracy. In order to analyze singular points, it would also be useful to consider slightly more general series of the form $z^\alpha(a_0+a_1z+a_2z^2+\cdots)$ for non-integral $\alpha$. –  whuber May 17 '13 at 18:07
Thank you for your answer; it is very clear. My equation is very similar to the one you use; its nonlinear term is of the form sin(2f). I am more familiar with Maple which has the basic command dsolve(equation, series) and I was sure that MMA must have such a command also. –  rick May 18 '13 at 2:02
@rick, that's a bit unclear; is the thing inside the sine the dependent or independent variable? –  Ｊ. Ｍ. May 18 '13 at 15:03
My equation has the term sin (2y(x)) so it is far from linear. When looking at the documentation for NDSolve I did not see see where one might specify the integration method to use-I was searching for a Taylor series option. –  rick May 19 '13 at 15:25