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.

  • $\begingroup$ Yes, there's a way. Do you have an example problem? $\endgroup$ – Mark McClure May 17 '13 at 3:45
  • $\begingroup$ Somewhat surprisingly, the procedure you want isn't built-in, but one could certainly write a routine for this task. $\endgroup$ – J. M.'s torpor May 17 '13 at 3:46
  • 1
    $\begingroup$ @J.M. Not true! It's buried in the Holonomic context $\endgroup$ – Mark McClure May 17 '13 at 3:50
  • 1
    $\begingroup$ @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$. $\endgroup$ – J. M.'s torpor May 17 '13 at 3:55
  • 2
    $\begingroup$ @Mark, after playing with some of those functions for a while, I'm not sure if using them is any better than forming a DifferentialRoot[] and then using Series[] for expansion (altho I am sure those functions are used under the hood with what I'm proposing). $\endgroup$ – J. M.'s torpor May 17 '13 at 5:06

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}]

enter image description here

  • 3
    $\begingroup$ +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$. $\endgroup$ – whuber May 17 '13 at 18:07
  • 1
    $\begingroup$ 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. $\endgroup$ – rick May 18 '13 at 2:02
  • $\begingroup$ @rick, that's a bit unclear; is the thing inside the sine the dependent or independent variable? $\endgroup$ – J. M.'s torpor May 18 '13 at 15:03
  • $\begingroup$ 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. $\endgroup$ – rick May 19 '13 at 15:25
  • $\begingroup$ This's a great answer, but I don't know why we can't compute the following differential equation case:ode = x[t]*x''[t] - 1 - x'[t]^2 == 0; initconds = {x[0] == 1, x'[0] == 0};. $\endgroup$ – A little mouse on the pampas Aug 23 '20 at 3:25

J.M. in a comment above alluded to this possibility. Using Series and DifferentialRoot seems like a great way to expand linear differential equations in series. It seems worth putting down in an answer due to its simplicity:

lde = {y''[x] - x*y'[x] + Sin[x] == 0, y'[0] == 1, y[0] == 0};
Series[DifferentialRoot[Function @@ {{y, x}, lde}][x], {x, 0, 15}]

(* x + x^5/120 + x^7/1260 + (29 x^9)/362880 + (13 x^11)/1995840 +
    (2861 x^13)/6227020800 + (4649 x^15)/163459296000 + O[x]^16 *)

Here's a way to develop a series solution for a quasilinear ODE with analytic coefficients, $y^{(n)}=f(x,y,y',\dots,y^{(n-1)})$. The main idea to fold Solve over the series expansion of the ODE solve for the derivatives of $y$, something like this:

  Flatten[{#1, Solve[#2 == 0 /. #1]}] &,
  ics,                                               (* initial conditions *)
  Series[ode /. Equal -> Subtract, {x, x0, n}][[3]]  (* series coefficients *)

Getting this idea to work for others without their having tweaking it takes some work. I've handled most of the things I've run into in the code below, but there are probably some things I have not tested. A malformed ODE may fail, probably ungracefully. Suggestions and comments are welcome if you run into other problems.

seriesDSolve[ode_, y_, iter : {x_, x0_, n_}, ics_: {}] := 
 Module[{dorder, coeff},
  dorder = Max[0, Cases[ode, Derivative[m_][y][x] :> m, Infinity]];
  (coeff = Fold[
         Solve[#2 == 0 /. #1, 
          Union@Cases[#2 /. #1, Derivative[m_][y][x0] /; m >= dorder, 
            Infinity]]}] &,
       SeriesCoefficient[ode /. Equal -> Subtract, {x, x0, k}],
       {k, 0, Max[n - dorder, 0]}]
    Series[y[x], iter] /. coeff) /; dorder > 0

Example (constructed from a comment by the OP):

 y''[x] + y'[x] == Sin[2 y[x]],
 y, {x, 0, 3},
 {y[0] -> a, y'[0] -> b}         (* symbolic initial values *)

Mathematica graphics

(Looking at it, I suppose it would be nice if the initial condition could be specified just like in DSolve. That would take some parsing and this approach evolved according to in-house needs -- that is to say, mine.)

  • 1
    $\begingroup$ As we know, any Nth order ordinary differential equation can be recast as a coupled set of N first order ODEs, but not every set of N first order ODEs can necessarily be written as a neat Nth order ODE. So , if there are a system of coupled ODEs, how can we use seriesDSolve to obtain the series solution? $\endgroup$ – Mark_Phys Nov 9 '16 at 15:41
  • 1
    $\begingroup$ @user10709 You'd have to change the code to handle a vector of dependent variables instead of just a single (scalar) y. That includes handling their different derivatives by adapting the pattern Derivative[m_][y][x0]. The Solve command would have to solve a system for the next higher-order derivatives of the variables at each stage. That seems feasible, but easier if it's a first order system of the form $dx/dt = F(x, t)$ (for vector $x$). As you say, a system could be put in such form. $\endgroup$ – Michael E2 Nov 9 '16 at 18:55
  • $\begingroup$ @MichaelE2, This is very useful code. I used it to determine series solution for second order nonlinear differential equation. Can you extend it to four coupled second order differential equation? $\endgroup$ – Aschoolar Jan 30 at 0:45

On Mathematica 11.3 using a new function AsymptoticDSolveValue

Linear ODE:

AsymptoticDSolveValue[{y''[x] - x*y'[x] + Sin[x] == 0, y'[0] == 1, 
y[0] == 0}, y[x], {x, 0, 10}]

(* x + x^5/120 + x^7/1260 + (29 x^9)/362880  *)

Quasilinear ODE:

AsymptoticDSolveValue[{y''[x] + y'[x] == Sin[2*y[x]], y[0] == a, 
y'[0] == b}, y[x], {x, 0, 5}]

(* a + b x + 1/6 x^3 (b + 2 b Cos[2 a] - Sin[2 a]) + 
1/2 x^2 (-b + Sin[2 a]) + 
1/24 x^4 (-b - 4 b Cos[2 a] + Sin[2 a] - 4 b^2 Sin[2 a] + 
2 Cos[2 a] Sin[2 a]) + 
1/120 x^5 (b + 6 b Cos[2 a] - 8 b^3 Cos[2 a] + 4 b Cos[2 a]^2 - 
Sin[2 a] + 16 b^2 Sin[2 a] - 4 Cos[2 a] Sin[2 a] - 
12 b Sin[2 a]^2)*)

Nonlinear ODE:

AsymptoticDSolveValue[{y'[x] == 1 + y[x]^2, y[0] == 0}, y[x], {x, 0, 10}]

(* x + x^3/3 + (2 x^5)/15 + (17 x^7)/315 + (62 x^9)/2835 *)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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