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, 2020 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, 2020 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, 2020 at 11:30
  • $\begingroup$ I'm starting from 0, a[0]==1, a[1]==1, a[2]==2 $\endgroup$ Oct 16, 2020 at 11:33

2 Answers 2

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, 2020 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, 2020 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, 2020 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, 2020 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, 2020 at 13:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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