Discretized normal distribution

Question

The normal distribution $e^{-x^2/2}$ satisfies the ODE $y'=-xy$ on the real line.

Is there a non-trivial, non-negative solution to the discretized version of this ODE $$u(n+1)-u(n-1)=-2nu(n)?$$

Answer 1

First, any solution $u$ satisfies $u(-n)=u(n)$ for all $n$ — because both $u(n)$ and $u(-n)$ satisfy the same recurrence, $u(0)=u(-0)$ and $u(1)=u(-1)\impliedby u(1)-u(-1)=-2\cdot 0\cdot u(0).$

Next, the function $$U(x)=\sum_{n=0}^{\infty}u(n)\frac{x^n}{n!},\qquad |x|<\frac12,$$ satisfies $(1+2x)U''(x)+2U'(x)-U(x)=0$ which is solved by $$U(x)=C_1 I_0(\sqrt{1+2x})+C_2 K_0(\sqrt{1+2x})$$ where $I_0$ and $K_0$ are modified Bessel functions (and $C_1, C_2$ are constants).

This gives general solutions of the recurrence in question. Further, as we have $$I_0(\sqrt{1+2x})=\sum_{k=0}^{\infty}\frac{1}{k!^2}\Big(\frac{1+2x}{4}\Big)^k\\=\sum_{k=0}^{\infty}\frac{2^{-2k}}{k!^2}\sum_{n=0}^{k}\frac{k!(2x)^n}{n!(k-n)!}\\=\sum_{n=0}^{\infty}\frac{(2x)^n}{n!}\sum_{k=n}^{\infty}\frac{2^{-2k}}{k!(k-n)!}\\=\sum_{n=0}^{\infty}\frac{1}{n!}\Big(\frac{x}{2}\Big)^n\sum_{k=0}^{\infty}\frac{2^{-2k}}{k!(n+k)!},$$ the choice $C_1=1, C_2=0$ gives a nonnegative solution: $$u(n)=\sum_{k=0}^{\infty}\frac{2^{-|n|-2k}}{k!(|n|+k)!}.$$ This is just $I_n(1)$; actually, the general solution can be seen as $$u(n)=C_1 I_n(1) + C_2 (-1)^n K_n(1).$$ Note that $I_n(1)$ decays (and $K_n(1)$ grows) when $n\to\infty$.

Answer 2

Edit : in fact, the sequence one gets through the given second order recurrence relationship is, whatever $u_0$ and $u_1$, most often oscillating and rapidly divergent.

I propose to take instead the following discretization with a non-symmetrical derivative (with a step $\Delta x=1$) :

$$-n u_n=u_n-u_{n-1} \ \ \ \text{with} \ u_0=1\tag{1}$$

giving the first order (instead of second order) recurrence relationship :

$$u_n=\tfrac{1}{1+n}u_{n-1}\tag{2}$$

One gets in this way a decreasing sequence tending to $0$ from above, already "in the same spirit" as Laplace-Gauss.

But a nice surprise comes, when one takes a much smaller step, for example $\Delta x = 0.01$, instead of step $\Delta x = 1$, i.e. by replacing (2) by :

$$u_n=\tfrac{1}{1+n/100}u_{n-1}\tag{3}$$

When one plots $u_n$ as a function of variable $n$, one obtains the following "curve", very close to the (right part of) Laplace-Gauss curve :

Fig. 1 : The 21 first terms of recurrent sequence defined by (3) $u_0=1, u_1=\tfrac{100}{101}, \cdots$ represented as a curve.

Remark : One should make a study of the "most adequate" $\Delta x$. We will not do it because we are already far enough from the initial question of the OP.

---

Previous "answer", merely a comment now.

For sequences, It is often advisable to take a look at the OEIS Encyclopedia of sequences. But in this case, I have only found results for the cousin sequence A001053 : https://oeis.org/A001053.

$a(n+1) = n*a(n) + a(n-1)$ with $a(0)=1, a(1)=0,$

mentionning solutions with modified Bessels, etc. (as some of those found by @Metamorphy and explained in his very interesting answer).

Answer 3

Clearly $u(-n)=u(n)$ by symmetry, and for any solution multiplying all the $u(n)$ by a constant $k$ also gives a solution. So in a sense, having a correct value for the ratio of $u(0)$ and $u(1)$ determines all the other values from the recursion.

As an empirical approach, it is clear that the problem is that once $u(n)$ is negative while $u(n-1)$ is positive, you then get explosive oscillation. If this is supposed to be like a normal distribution then values should not be negative. So the aim should be to find a solution which prevents (or at least delays) a negative $u(n)$.

If you try $u(0)=1,u(1)=1$ then you get $u(2)=-1$, which is unacceptable. If you try $u(0)=1,u(1)=\frac12$ then you get $u(2)=0,u(3)=\frac12, u(4)=-3$ which is better. Then try $u(0)=1,u(1)=\frac14$ and you get $u(2)=\frac12,u(3)=-\frac74$ which goes too far, but suggests that you want $\frac14 \lt \frac{u(1)}{u(0)}\lt \frac12$.

You can then use a binary search algorithm to improve the calculation of $\frac{u(1)}{u(0)}$ to arbitrary precision, adjusting the estimate based on when the first negative value occurs from what I have elsewhere called failed recursions. Since you only want the ratio $\frac{u(1)}{u(0)}$, you can keep the calculation in integers to avoid rounding issues, expressing each step as:

• double the previous value used for $u(0)$
• double the previous value used for $u(1)$ and either adjust it by $+1$ if in the previous step the negative values occur for $u(n)$ with $n$ odd, or adjust it by $-1$ if in the previous step the negative values occur for $u(n)$ with $n$ even
• inspect the values for later $u(n)$ to find which are negative for the adjustment in the next step

You should get something like this

       u(0)    u(1)   u(2)    u(3)   u(4)   u(5)     u(6)     u(7)   next adj.

         1        1     -1       5    -31    253    -2561    30985      -1
         2        1      0       1     -6     49     -496     6001      -1
         4        1      2      -7     44   -359     3634   -43967      +1
         8        3      2      -5     32   -261     2642   -31965      +1
        16        7      2      -1      8    -65      658    -7961      +1
        32       15      2       7    -40    327    -3310    40047      -1
        64       29      6       5    -24    197    -1994    24125      -1
       128       57     14       1      8    -63      638    -7719      +1
       256      115     26      11    -40    331    -3350    40531      -1
       512      229     54      13    -24    205    -2074    25093      -1
      1024      457    110      17      8    -47      478    -5783      +1
      2048      915    218      43    -40    363    -3670    44403      -1
      4096     1829    438      77    -24    269    -2714    32837      -1
      8192     3657    878     145      8     81     -802     9705      -1
     16384     7313   1758     281     72   -295     3022   -36559      +1
     32768    14627   3514     571     88   -133     1418   -17149      +1
     65536    29255   7026    1151    120    191    -1790    21671      -1
    131072    58509  14054    2293    296    -75     1046   -12627      +1
    262144   117019  28106    4595    536    307    -2534    30715      -1
    524288   234037  56214    9181   1128    157     -442     5461      -1
   1048576   468073 112430   18353   2312   -143     3742   -45047      +1
   2097152   936147 224858   36715   4568    171     2858   -34125      +1
   4194304  1872295 449714   73439   9080    799     1090   -12281      +1
   8388608  3744591 899426  146887  18104   2055    -2446    31407      -1

which suggests $\frac{u(1)}{u(0)}$ should be marginally less than $0.44639$, and by taking more steps you can find this ratio to as high a precision as you wish. metamorphy's analysis shows this is the ratio of modified Bessel functions $\frac{I_1(1)}{I_0(1)}=0.4463899658965345\ldots$.

If you want to turn this into a probability mass function, i.e. summing to $1$, you then get $u(0) \approx 0.466$, $u(1)=u(-1) \approx 0.208$, $u(2)=u(-2) \approx 0.050$, $u(3)=u(-3) \approx 0.008$, $u(4)=u(-4) \approx 0.001$, and other $u(n)$ much smaller. These correspond to the so-called exponentially scaled Bessel functions $u(n)=\frac1e I_n(1)$, which can be derived from metamorphy's analysis.

Comparing this probability mass function in blue with the red density of a standard normal distribution gives the following chart. Despite appearances, the variances of the two distributions are both $1$, though the discrete distribution is leptokurtic with kurtosis $4$ rather than the normal distribution's $3$ (i.e. excess kurtosis $1$ rather than $0$).