# Why earlier terms generated from AsymptoticDSolveValue change when increasing the order?

I was comparing my hand solution with Mathematica's. I noticed the particular solution terms generated for a Frobenius series solution change depending on the order. This should not happen. As when increasesing the order, more terms will show up, but the previous terms should not change their earlier values.

My question is why does this happen? And how could this be explained? Since when solving Frobenius series by hand this should not happen. (I can show my hand solution if needed).

Here is an example. The homogeneous solution to this ode is correct and increasing the order do not cause the earlier terms to change. It is only the particular solution which have this behavior.

ode = 2 x*y''[x] + (x + 1) y'[x] + 3 y[x] == x;
sol = AsymptoticDSolveValue[ode, y[x], {x, 0, 2}] // Expand


Here is screen shot, showing only the particular solution and how the earlier term change as more terms are added. I am not showing the $$y_h$$ as that is correct and not relevant.

I then compare this with Maple 2022.2 by changing the order and asking for series solution, and noticed there the terms do not change as the order increases

Order:=5:
ode:=2*x*diff(y(x),x$2)+(x+1)*diff(y(x),x)+3*y(x)=x: expand(convert(dsolve(ode,y(x),'series'),polynom))  You see that the terms do not change, but more terms are added as the order increases. This is what is expected. Why does the terms in particular solution change with AsymptoticDSolveValue? From help, it says it supports Frobenius series method. V 13.2 on windows 10. • Looks like a bug to me ! Commented Feb 27, 2023 at 20:46 • All the terms up to the order requested seem correct, though. You could add O[x]^(n+1) for order n to truncate the spurious terms (followed by // Normal // Expand if desired). Commented Feb 27, 2023 at 21:24 • @MichaelE2 do you happen to know/guess why this issue with the spurious terms affects only the particular solution and not the homog. solution? i.e. homog. solution terms do not change as order is increased. Commented Feb 27, 2023 at 22:50 • I don't know. It looks like there's quadratic term in the ode (the 2n in the 2n+1 terms you're getting), but I don't see how that makes sense. Or maybe it's miscounting the number of terms to compute because of the half-powers in the Frobenius series. I assume the homog. sol. is computed separately and correctly. It just messes up on the particular sol. Commented Feb 28, 2023 at 2:14 • The IVP {ode, y[0] == 0, y'[0] == 0} does not produce the spurious terms. Odd. Commented Feb 28, 2023 at 13:50 ## 3 Answers Here's what seems to be going on: variation of parameters is used to derive a particular solution: ode = 2 x*y''[x] + (x + 1) y'[x] + 3 y[x] == x; sol = AsymptoticDSolveValue[ode, y[x], {x, 0, 3}] // Collect[#, {C[1], C[2]}, Expand] &  {y1, y2} = Coefficient[sol, {C[1], C[2]}]; fx = First@SolveValues[ode /. {y[x] -> 0, y'[x] -> 0}, y''[x]]; -y1 Integrate[ Normal@Series[fx y2/Wronskian[{y1, y2}, x], {x, 0, 3}], x] + y2 Integrate[ Normal@Series[fx y1/Wronskian[{y1, y2}, x], {x, 0, 3}], x] % // Expand  It seems WRI merely neglected to truncate the result to the specified order. The terms beyond the order requested are not accurate terms of the series solution because the factors in the variation of parameters formula do not have accurate terms beyond the requested order. $Version

(* "13.2.1 for Mac OS X x86 (64-bit) (January 27, 2023)" *)

ode = 2 x*y''[x] + (x + 1) y'[x] + 3 y[x] == x;


The exact solution with the arbitrary constants set to zero is

sol = DSolveValue[ode, y[x], x] /. {C[1] -> 0, C[2] -> 0} // Simplify

(* 1/96 x (15 + x) -
1/96 E^(-x/2) Sqrt[\[Pi]/2] Sqrt[x] (15 - 10 x + x^2) Erfi[Sqrt[x]/Sqrt[2]] *)


This satisfies the ode

(ode /. y -> Function[{x}, Evaluate@sol]) // Simplify

(* True *)


If you knew the exact solution then the series expansion would be

series = Series[sol, {x, 0, 13}] // Normal

(* x^2/6 - x^3/18 + x^4/84 - x^5/540 + x^6/4455 - x^7/45045 + x^8/540540 - \
x^9/7518420 + x^10/119041650 - x^11/2115278550 + x^12/41701205700 - \
x^13/903526123500 *)


You should only expect the coefficients to be the same as the series up to the requested number of terms. Worrying about anything beyond the requested order is like worrying about the digits of a number beyond the specified precision.

Grid[
Join[
{{ReplacePart[Take[series, 12], 12 -> \[Ellipsis]],
SpanFromLeft}}, ({#, Total[{#[[1]], Style[#[[-1]], Red]}] &@
(Total /@
Partition[
List @@ Expand[
AsymptoticDSolveValue[
ode, y[x], {x, 0, #}] /. {C[1] -> 0, C[2] -> 0}],
UpTo[#]])}) & /@ Range[6]],
Alignment -> {{Center, Left}, Center},
Frame -> All]


• You should only expect the coefficients to be the same as the series up to the requested number of terms. exactly. Then why does it return terms not up to the order requested in this case? Asking the user to count manually the terms in order to figure when to stop counting is not correct behavior. Maple only returns terms up requested order and no more. Commented Feb 27, 2023 at 22:42
• But do you happen to know why this affects only the particular solution and not the homog. solution? i.e. homog. solution terms do not change as order is increased. Commented Feb 27, 2023 at 22:49
• If you don't want the other terms use Take[ AsymptoticDSolveValue[ode, y[x], {x, 0, #}] // Expand, #] & /@ Range[6] Commented Feb 27, 2023 at 22:59

Inspired by the trick shown by Michael in the comment that when using using zero IC, spurious terms are gone, this is a function which will do this automatically. It only works with second order ode's that are meant to be solved with NO IC's.

Now all terms in $$y_p$$ show up are at the correct order as $$y_h$$ and the user does not need to figure which ones in $$y_p$$ are correct and which are not.

mySeriesSolution[odeIn_Equal, y_, x_, order_Integer?(Positive), at_ : 0] :=
Module[{ode = First@odeIn - Last@odeIn, lhs, rhs, z, solyp, ysol},
rhs = Total@Cases[ode, z_ /; FreeQ[z, y]];
lhs = Total@Cases[ode, z_ /; Not@FreeQ[z, y]];
ysol = AsymptoticDSolveValue[lhs == 0, y[x], {x, at, order}];
If[rhs =!= 0,
solyp = AsymptoticDSolveValue[{lhs == rhs, y[0] == 0, y'[0] == 0}, y[x], {x, at, order}];
ysol += solyp
];
ysol
]


Using the above on the example in the question gives

ode = 2 x*y''[x] + (x + 1) y'[x] + 3 y[x] == x;
mySeriesSolution[ode, y, x, 2]
mySeriesSolution[ode, y, x, 3]
mySeriesSolution[ode, y, x, 4]
mySeriesSolution[ode, y, x, 5]


Compare to direct call to AsymptoticDSolveValue

ode = 2 x*y''[x] + (x + 1) y'[x] + 3 y[x] == x;
Collect[Expand@AsymptoticDSolveValue[ode, y[x], {x, 0, 2}], {C[1], C[2]}]
Collect[Expand@AsymptoticDSolveValue[ode, y[x], {x, 0, 3}], {C[1],  C[2]}]
Collect[Expand@AsymptoticDSolveValue[ode, y[x], {x, 0, 4}], {C[1],  C[2]}]
Collect[Expand@AsymptoticDSolveValue[ode, y[x], {x, 0, 5}], {C[1],  C[2]}]


Where all terms in red should not be there.