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$
$\endgroup$
6
4 Answers
$\begingroup$
$\endgroup$
5
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}]
answered May 17, 2013 at 15:52
-
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$– whuberCommented May 17, 2013 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$– rickCommented May 18, 2013 at 2:02
-
$\begingroup$ @rick, that's a bit unclear; is the thing inside the sine the dependent or independent variable? $\endgroup$ Commented May 18, 2013 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$– rickCommented May 19, 2013 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$ Commented Aug 23, 2020 at 3:25
$\begingroup$
$\endgroup$
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 *)
answered Sep 18, 2014 at 1:01
$\begingroup$
$\endgroup$
3
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:
Fold[
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.
Clear[seriesDSolve];
seriesDSolve[ode_, y_, iter : {x_, x0_, n_}, ics_: {}] :=
Module[{dorder, coeff},
dorder = Max[0, Cases[ode, Derivative[m_][y][x] :> m, Infinity]];
(coeff = Fold[
Flatten[{#1,
Solve[#2 == 0 /. #1,
Union@Cases[#2 /. #1, Derivative[m_][y][x0] /; m >= dorder,
Infinity]]}] &,
ics,
Table[
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):
seriesDSolve[
y''[x] + y'[x] == Sin[2 y[x]],
y, {x, 0, 3},
{y[0] -> a, y'[0] -> b} (* symbolic initial values *)
]
(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.)
answered Oct 10, 2015 at 19:51
-
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
seriesDSolveto obtain the series solution? $\endgroup$ Commented Nov 9, 2016 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 patternDerivative[m_][y][x0]. TheSolvecommand 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$ Commented Nov 9, 2016 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$ Commented Jan 30, 2021 at 0:45
$\begingroup$
$\endgroup$
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 *)
answered Mar 25, 2018 at 15:37
Holonomiccontext $\endgroup$DifferentialRoot[]and then usingSeries[]for expansion (altho I am sure those functions are used under the hood with what I'm proposing). $\endgroup$