2
$\begingroup$

I am trying to solve a recurrence relation using generating functions method: $$a_n=a_{n-1}+(n-1)a_{n-2}+(0.5n^2-1.5n+1)a_{n-3}$$

After some long calculations, I have arrived to this second-order differential equation: $$0.5 x^5 y''(x)+(2x^4+x^3)y'(x)+\left(x^3+x^2+x-1\right)y(x)+1=0$$

and these conditions: $y(0)=1, y'(0)=1$. $y(x)$ is the function that needs to be expanded as Taylor Series at $x=0$ to obtain the sequence from the coefficients. However, when I try to solve it both using DSolve and NDSolve, I have no luck. With DSolve it just returns the request itself:

$$\text{DSolve}\left[\left\{0.5 x^5 y''(x)+(2. x+1) x^3 y'(x)+\left(1. x^3+x^2+x-1\right)y(x)+1=0,y(0)=1,y'(0)=1\right\},y,x\right]$$

And with NDSolve I just receive errors and no equation:

Power::infy: Infinite expression 1/0.^5 encountered.
Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered.
NDSolve::ndnum: Encountered non-numerical value for a derivative at x == 0.`.

$$\text{NDSolve}\left[\left\{0.5 x^5 y''(x)+(2. x+1) x^3 y'(x)+\left(1. x^3+x^2+x-1\right)y(x)+1=0,y(0)=1,y'(0)=1\right\},y,\{x,0,1\}\right]$$

How could I resolve this problem?

Code:

Simplify[y[x] - (1 + x + 2 x^2)]
l = Expand[%]
Simplify[x (y[x] - (1 + x))]
r1 = Expand[%]
Simplify[x*D[x^2 (y[x] - 1), x] - x^2 (y[x] - 1)]
r2 = Expand[%]
Simplify[0.5 x*D[x*D[x^3*y[x], x], x] - 1.5 x*D[x^3*y[x], x] + 
  x^3*y[x]]
r3 = Expand[%]
eq = FullSimplify[r1 + r2 + r3 - l]
DSolve[{eq == 0, y[0] == 1}, y, x]
NDSolve[{eq == 0, y[0] == 1, y'[0] == 1}, y, {x, 0, 1}]
$\endgroup$
4
  • $\begingroup$ Please post code that can be readily copy-pasted by anyone who wants to try and answer your question. $\endgroup$ Oct 14 '20 at 17:44
  • $\begingroup$ NDSolve fails here, because the ODE is singular at x = 0. Try integrating from a slightly larger value of x, but even then expect some difficulties. $\endgroup$
    – bbgodfrey
    Oct 14 '20 at 18:47
  • $\begingroup$ RecurrenceTable[{a[n] == a[n - 1] + (n - 1) a[n - 2] + (n^2 - 3 n + 2)/2 a[n - 3], a[1] == 1, a[2] == 1, a[3] == 2}, a, {n, 11}] yields {1, 1, 2, 8, 22, 82, 334, 1370, 6338, 30692, 155722}. $\endgroup$
    – bbgodfrey
    Oct 15 '20 at 11:30
  • $\begingroup$ I'm starting from 0, a[0]==1, a[1]==1, a[2]==2 $\endgroup$ Oct 16 '20 at 11:33
5
$\begingroup$

Error Corrected.

Because neither DSolve nor AsymptoticDSolveValue produced a solution, I expanded the ODE itself and computed the coefficients of each term.

nmax = 10;
Normal@Series[x^5  y''[x]/2 + (2 x + 1) x^3 y'[x] + (x^3 + x^2 + x - 1) y[x] + 1, 
    {x, 0, nmax}]
CoefficientList[%, x] /. Derivative[n_][y][0] -> c[n] /. y[0] -> c[0]
Flatten@Solve[Thread[% == 0], Array[c, nmax + 1, 0]] // Values
%/(Range[0, nmax]!)

(* {1, 1, 2, 5, 14, 46, 166, 652, 2780, 12644, 61136} *)
$\endgroup$
0
3
$\begingroup$

When I rationalize eq, AsymptoticDSolveValue returns a solution for me:

eq = FullSimplify[r1 + r2 + r3 - l]
1 + (-1 + x + x^2 + x^3) y[x] + x^3 (1 + 2 x) y'[x] + 1/2 x^5 y''[x]
seq = CoefficientList[
 AsymptoticDSolveValue[{eq == 0, y[0] == 1}, y, {x, 0, 30}],
 x]
(*
{1, 1, 2, 5, 14, 46, 166, 652, 2780, 12644, 61136, 312676, 1680592, \
9467680, 55704104, 341185496, 2170853456, 14314313872, 97620050080, \
687418278544, 4989946902176, 37286121988256, 286432845428192, \
2259405263572480, 18280749571449664, 151561941235370176, \
1286402259593355776, 11168256342434121152, 99099358725069658880, \
898070590439513534464, 8306264068494786829696}
*)

And FindSequenceFunction returns the original recursion:

an = FindSequenceFunction[%]
(*
DifferenceRoot[
 Function[{\[FormalY], \[FormalN]}, {-\[FormalN] (1 + \[FormalN]) \
\[FormalY][\[FormalN]] + (-2 - 2 \[FormalN]) \[FormalY][
       1 + \[FormalN]] - 2 \[FormalY][2 + \[FormalN]] + 
     2 \[FormalY][3 + \[FormalN]] == 0, \[FormalY][1] == 
    1, \[FormalY][2] == 1, \[FormalY][3] == 2}]]
*)
Table[an[n+1], {n, 0, 30}]
(*  same as above: 1, 1, 2, 5, 14...  *)

I guess I'm not sure what exactly is being sought, though.

$\endgroup$
8
  • $\begingroup$ Nice soluton (+1). I tried AsymptoticDSolveValue[{x^5 y''[x]/2 + (2 x + 1) x^3 y'[ x] + (x^3 + x^2 + x - 1) y[x] + 1 == 0, y[0] == 1, y'[0] == 1}, y, {x, 0, 1}] instead, which returned unevaluated. Any idea why inclusion of y'[0] == 1 caused AsymptoticDSolveValue to fail? $\endgroup$
    – bbgodfrey
    Oct 14 '20 at 21:47
  • $\begingroup$ @bbgodfrey I guess it has to do with $x=0$ being an irregular singular point. I'm not well-acquainted irregular singular points. But try a higher order, such as AsymptoticDSolveValue[{eq == 0, y[0] == 1, y'[0] == 1}, y, {x, 0, 2}]. (I don't actually recall removing y'[0] == 1, and it may have been carelessness on my part in copying code. But it worked and I didn't catch it.) $\endgroup$
    – Michael E2
    Oct 15 '20 at 1:49
  • $\begingroup$ Amazing! {x, 0, 1} and x -> 0 don't work, but {x, 0, 2} and higher do. A bug? $\endgroup$
    – bbgodfrey
    Oct 15 '20 at 2:23
  • $\begingroup$ @bbgodfrey I don't know (whether a bug). Possibly a heuristic about how far to develop the series to solve for the order 1 leads to some problem; however, SeriesTermGoal does not help. $\endgroup$
    – Michael E2
    Oct 15 '20 at 2:48
  • 1
    $\begingroup$ @MichaelMykhaylov Did you "rationalize" the floating-point coefficients (i.e., change 0.5 to 1/2 and 1.5 to 3/2)? (bbgodfrey did this, too.) I gave the value of eq that I used in the first output line in my answer. You should use this and avoid floating point. Floating point cannot represent integers exactly above around 10^16, and Mathematica displays only the first six digits of FP numbers. You can see the complete number if you convert it to InputForm or FullForm. $\endgroup$
    – Michael E2
    Oct 15 '20 at 13:35

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.