1

Higher Derivatives of NDSolve don't match the results from DSolve solution, can I somehow increase the accuracy of numerical method to match arbitrarily high derivatives?

Equation

I am using SeriesCoefficient to calculate probability distribution from the probability generating function $G(x)$ which should satisfy this parametric equation: $$ 4 x (x - 1)^2 G''(x)+((4 d - 4 L + 6) x^2 - (8 d - 4 L + 10)x + 4 (d + 1)) G'(x) +(L - d)(L - d - 1)(x - 1) G(x) = 0 $$ With boundary conditions $G(0)=0,G(1)=1$, where $G(x) = \sum_{n=d}^{L}P_{n}x^n \rightarrow P(n)=\frac{G^{(n)}(0)}{n!}$

Code

    sol = ParametricNDSolveValue[
    {4 x (x - 1)^2 f''[x] + 
    ((4 d - 4 L + 6) x^2 + 4 (d + 1) - (8 d - 4 L + 10) x) f'[x] + 
    ((L - d) (L - d - 1)) (x - 1) f[x] == 0, 
    {f[0] == 0, f[1] == 1}},f, {x, 0, 1}, {L, d}, 
    Method -> "FiniteElement"];

G[x_] = (x^d Hypergeometric2F1[d/2 - L/2, 1/2 + d/2 - L/2, 1 + d, x])/
Hypergeometric2F1[d/2 - L/2, 1/2 + d/2 - L/2, 1 + d, 1];

Manipulate[{
    Plot[{
        sol[l, k][x], G[x] /. {L -> l, d -> k}}, {x, 0, 1},
    PlotRange -> {0, 1}, PlotStyle -> {Thick, Dashed}],
    ListLinePlot[{
        Table[SeriesCoefficient[sol[l, k][x] /. {L -> l, d -> k}, {x, 0, n}], {n, 0, l}], 
        Table[SeriesCoefficient[G[x] /. {L -> l, d -> k}, {x, 0, n}], {n,0, l}]} ,
    PlotRange -> {0, 1}, PlotStyle -> {Thick, Dashed}]},
 {l, 0, 25, 1}, {k, 0, 25, 1}]

Problem

  1. I don't know if I can increase precision of the interpolated function of NDSolve to match the results of the analytical solution. Parameter $L$ defines the size of the system, and it can assume the values from $0$ to an arbitrarily large number.
    If $L=100$ then I need to calculate up to the 100th derivative $\frac{\partial^{100}}{\partial x^{100}}G(x)$. There are significant slowdowns for $L \rightarrow l >50$ already, and solutions also don't match when varying $d \rightarrow k$ parameter
  2. Alternatively I would like to know what other numerical methods I could use for this kind of ODE or PDE besides those available in Mathematica, because only Method -> "FiniteElement" gives some result when concerning BVP problem with singularities.

Related

highest-derivative-from-ndsolve
series-expansion-of-interpolatingfunction

Using InterpolationOrder->All or Method -> "ExplicitRungeKutta" and increasing "DifferenceOrder" don't seeem to work.

SvenDeanson
  • 11
  • 1
  • Assuming you've got a grip on the singularities, which I haven't, the series expansion will be limited by the interpolation order at the step at which the series is expanded, which for FEM is 2 and for most (all?) ODE methods is 2n+1 at the first step, where n is the differential order of the ODE. Reinterpolation seems to work here: https://mathematica.stackexchange.com/questions/133531/sudden-truncation-of-numerical-series-expansion -- I rather think getting the 100th derivative numerically is a tall order. – Michael E2 Sep 30 '21 at 11:20
  • BTW, 4 x (x - 1)^2 f''[x] + ((4 d - 4 L + 6) x^2 + 4 (d + 1) - (8 d - 4 L + 10) x) f'[x] + ((L - d) (L - d - 1)) (x - 1) f[x] == 0 /. f -> G /. {L -> 50, d -> 20, x -> 1/2} // Simplify yields False, suggesting that G is not a solution to the ODE. – Michael E2 Sep 30 '21 at 13:22

0 Answers0