One asks for the probability $r$ starting at $1$ to eventually reach $0$. The dynamics is invariant by translations hence $r$ is also the probability starting at $2$ to eventually reach $1$. Consider the first step of a random walk starting at $1$. Either the first step is to $0$ then one hits $0$ eventually since one is already at $0$. Or the first step is to $2$ then to hit $0$, one must first hit $1$ starting from $2$ and after that, hit $0$ starting from $1$. This yields the equation $r=\frac13+\frac23r^2$, whose solutions are $r=1$ and $r=\frac12$.
If $r=1$, let us assume the random walk continues with the same dynamics after its first return to $0$. The new portion of the walk is distributed as before hence one returns to $0$ a second time. And so on, hence, calling $X_n$ the position at time $n$, one sees that $X_n=0$ for infinitely many times $n$.
Consider now the homogeneous random walk $(Y_n)$ on the whole integer line whose steps are $+1$ and $-1$ with probabilities $\frac23$ and $\frac13$ respectively. Then $Y_n=Z_1+\cdots+Z_n$ where $(Z_n)$ is i.i.d. and $\mathrm E(Z_n)=\frac13(-1)+\frac23(+1)=\frac13\gt0$. By the strong law of large numbers, $\frac1nY_n\to\frac13$. One can recover $(X_n)$ from $(Y_n)$ through the change of time $\tau_{n+1}=\min\{k\gt\tau_n\mid Y_k\geqslant0,\,Y_k\ne Y_{\tau_n}\}$ for every $n\geqslant0$, and $\tau_0=0$. Then $X_n=Y_{\tau_n}$ and $\tau_n\geqslant n$ hence $\frac1nX_n=\frac1nY_{\tau_n}\geqslant\frac1{\tau_n}Y_{\tau_n}$. One sees that $\liminf\limits_{n\to\infty}\frac1nX_n\geqslant\lim\limits_{n\to\infty}\frac1nY_n=\frac13$.
This is impossible if $X_n=0$ infinitely often, hence $r=\frac12$.
For a random walk whose steps are $+1$ and $-1$ with probability $p\gt\frac12$ and $1-p$ respectively, the same argument yields $r=\frac{1-p}p$.
For a random walk whose steps are $+2$ and $-1$ with probability $p$ and $1-p$ respectively, the crucial argument that to reach $0$ from $n\gt0$, one must reach $n-1$ from $n$, then reach $n-2$ from $n-1$, and so on, is still valid. Hence $r=pr^3+1-p$. If $r\gt\frac13$ the drift $p(+2)+(1-p)(-1)=3p-1$ is positive and one knows that $r\lt1$ hence $r$ is the positive root of the equation $p(r^2+r)=1-p$, that is, $r=\frac1{2p}\left(\sqrt{4p-3p^2}-p\right)$.
The same trick can be applied to any random walk whose steps are almost surely $\geqslant-1$.
First let's consider a closely related enumeration problem.
Define the range of a cycle $\{x_0, x_1, x_2, \ldots, x_{k-1}, x_k=x_0\}$ to be $\max_i \left|x_{i+1}-x_i\right|$. How many cycles on $\mathbb{Z}$ are there with $k$ edges and range $\le r$? (Two cycles are equivalent if they differ by a translation or a rotation.)
If $N_k^{(r)}$ is the number of cycles with $k$ edges and range no greater than $r$, then the probability of returning safely to $0$ by making random jumps of length $\le r$ is given by
$$
P^{(r)}=\sum_{k=2}^{\infty}\frac{kN_{k}^{(r)}}{(2r)^{k}}.
$$
The factor of $k$ is needed because the starting point ($0$) can be identified with any of the nodes in a given cycle; and the factor of $(2r)^{-k}$ is the probability of taking any specific sequence of $k$ jumps. For $r=1$, there is only one legal cycle (a jump to the right and then a jump to the left), and it has two edges, so
$$
P^{(1)}=\frac{2 N_2^{(1)}}{(2r)^{2}}=\frac{2\cdot 1}{(2\cdot 1)^2}=\frac{1}{2}.
$$
For $r=2$, there are two cycles for each $k\ge 2$ (depending on whether the first rightward jump is $+1$ or $+2$), so
$$
P^{(2)}=\sum_{k=2}^{\infty}\frac{kN_k^{(2)}}{(2r)^k}=2 \sum_{k=2}^{\infty}k 4^{-k}=-2x\frac{d}{dx}\sum_{k=2}^{\infty}x^{-k}\bigg\vert_{x=4}=-2x\frac{d}{dx}\left(\frac{1}{x(x-1)}\right)\bigg\vert_{x=4}=\frac{2(2x-1)}{x(x-1)^2}\bigg\vert_{x=4}=\frac{2\cdot 7}{4\cdot 9}=\frac{7}{18}.
$$
This recovers the two known results. For $r\ge 3$, the desired enumeration can be done using a transfer matrix method. I'll make this answer "community wiki" in case someone (maybe me) has the time and energy to fill in the details for $r=3$.
The gist of the transfer matrix approach is that we can list the "states" that a particular node in the cycle can have, figure out which state-to-state transitions are allowed moving from left to right on the number line, and express the number of ways to get from the left end of a cycle to the right end (i.e., the number of different cycles) in terms of a product of many identical local matrices. In this case the "state" consists of the set of jumps to, from, and over the node, along with the range remaining to each jump. Because we're counting cycles, each node will have exactly one incoming edge and one outgoing edge. The states for $r=2$ are simple:
A blue (red) line represents a jump to the right (left), and the number next to a line indicates its remaining range. A legal transition consists of these steps:
- Decrement all numbers on the right edge by the same value (at least $1$, but no more than the smallest number on the tile).
- Select a new tile with the same number of blue and red lines on its left edge as the current tile has on its right edge.
- Place the new tile to the right of the current tile and connect lines of the same color. Make sure that all lines that pass through the new tile maintain their (decremented) values.
Legal transitions for $r=2$ are: $A$ to $(B,C,D)$ (and $A$ to $D$ can happen two ways), $B$ to $(C,D)$, and $C$ to $(B,D)$. The corresponding transfer matrix is
$$
\left(\begin{matrix}0 & 0 & 0 & 0 \\
1 & 0 & 1 & 0 \\
1 & 1 & 0 & 0 \\
2 & 1 & 1 & 0
\end{matrix}\right).
$$
Now, the number of cycles with $k$ edges will be given by
$$
N^{(2)}_k = \left(\begin{matrix}0 & 0 & 0 & 1 \end{matrix}\right)
\left(\begin{matrix}0 & 0 & 0 & 0 \\
1 & 0 & 1 & 0 \\
1 & 1 & 0 & 0 \\
2 & 1 & 1 & 0
\end{matrix}\right)^{k-1}
\left(\begin{matrix}1 \\ 0 \\ 0 \\ 0
\end{matrix}\right)
=
\left(\begin{matrix}0 & 0 & 1 \end{matrix}\right)
\left(\begin{matrix}
0 & 1 & 0 \\
1 & 0 & 0 \\
1 & 1 & 0
\end{matrix}\right)^{k-2}
\left(\begin{matrix}
1 \\ 1 \\ 2
\end{matrix}\right),
$$
where we can eliminate the first row and column of the full transfer matrix because the initial state ($A$) never recurs (i.e., the first row is all zeroes). Note that we are counting the number of ways to go from the initial state to the final state, $D$, in $k-1$ steps. In this case, the vector on the right-hand side is an eigenvector of the (reduced) transfer matrix with eigenvalue $1$, so
$$
N^{(2)}_k=\left(\begin{matrix}0 & 0 & 1 \end{matrix}\right)
\left(\begin{matrix}
1 \\ 1 \\ 2
\end{matrix}\right) = 2
$$
for any $k \ge 2$. In general, though, we will have $N_k=\sum_i c_i \lambda_i^{k-2}$, where the $\lambda_i$ are eigenvalues of the reduced transfer matrix.
The $r=2$ case can be made even simpler by noting that the set of states is symmetric under exchange of red and blue ($A$ and $D$ are themselves symmetric, and $B$ and $C$ are symmetric images of each other), so we can treat mirror pairs as single states. The transitions are then: $A$ to $D$ (two ways), $A$ to $B/C$ (two ways), $B/C$ to $B/C$, and $B/C$ to $D$. So
$$
N^{(2)}_k = \left(\begin{matrix}0 & 0 & 1 \end{matrix}\right)
\left(\begin{matrix}0 & 0 & 0 \\
2 & 1 & 0 \\
2 & 1 & 0
\end{matrix}\right)^{k-1}
\left(\begin{matrix}1 \\ 0 \\ 0
\end{matrix}\right) = \left(\begin{matrix}0 & 1 \end{matrix}\right)
\left(\begin{matrix}
1 & 0 \\
1 & 0
\end{matrix}\right)^{k-2}
\left(\begin{matrix}
2 \\ 2
\end{matrix}\right) \\ = \left(\begin{matrix}0 & 1 \end{matrix}\right)\left(\begin{matrix}
2 \\ 2
\end{matrix}\right) = 2.
$$
This additional simplification, though hardly necessary for $r=2$, will come in handy for $r=3$.
For a less trivial application of the transfer matrix approach (but still not the full $r=3$ case), I considered the asymmetric model where each jump is uniformly drawn from $\{-2,-1,1,2,3\}$. The states in this case are:
The legal transitions are: $A$ to $(B,C,D,E,H)$ (and $A$ to $H$ can happen two ways), $B$ to $(C,D,H)$ (two way from $B$ to $H$), $C$ to $(D,H)$, $D$ to $(B,H)$, $E$ to $(F,G)$, $F$ to $G$, and $G$ to $H$. So the number of cycles with $k$ edges is given by
$$
\left(\begin{matrix}0 & 0 & 0 & 0 & 0 & 0 & 1 \end{matrix}\right)\left(\begin{matrix}
0 & 0 & 1 & 0 & 0 & 0 & 0\\
1 & 0 & 0 & 0 & 0 & 0 & 0\\
1 & 1 & 0 & 0 & 0 & 0 & 0\\
0 & 0 & 0 & 0 & 0 & 0 & 0\\
0 & 0 & 0 & 1 & 0 & 0 & 0\\
0 & 0 & 0 & 1 & 1 & 0 & 0\\
2 & 1 & 1 & 0 & 0 & 1 & 0\\
\end{matrix}\right)^{k-2}
\left(\begin{matrix}
1 \\ 1 \\ 1 \\ 1 \\ 0 \\ 0 \\ 2
\end{matrix}\right).
$$
Using the Jordan decomposition of the reduced transfer matrix, this is
$$
\left(\begin{matrix}1 & 0 & 0 & 0 & 1 & 1 & 1 \end{matrix}\right)\left(\begin{matrix}
0 & 1 & 0 & 0 & 0 & 0 & 0\\
0 & 0 & 1 & 0 & 0 & 0 & 0\\
0 & 0 & 0 & 1 & 0 & 0 & 0\\
0 & 0 & 0 & 0 & 0 & 0 & 0\\
0 & 0 & 0 & 0 & \nu_1 & 0 & 0\\
0 & 0 & 0 & 0 & 0 & \nu_2 & 0\\
0 & 0 & 0 & 0 & 0 & 0 & \nu_2^{*}
\end{matrix}\right)^{k-2}\left(\begin{matrix}
-1 \\ 0 \\ 1 \\ 1 \\ c_1 \\ c_2 \\ c_2^{*}
\end{matrix}\right) \\ =c_1 \nu_1^{k-2} + c_2 \nu_2^{k-2} + c_2^{*} (\nu_2^{*})^{k-2} - \delta_{k,2} + \delta_{k,4} + \delta_{k,5},
$$
where $\nu_1=1.32472$, $\nu_2=-0.662359 - 0.56228 i$, $c_1=2.945975$, and $c_2=0.0270124 + 0.118443 i$. Then the probability of safe return to $0$ is given by
$$
P^{(3,2)}=\sum_{k=2}^{\infty}\frac{kN_k^{(3,2)}}{5^k}=-\frac{2}{5^2}+\frac{4}{5^4}+\frac{5}{5^5}+\sum_{i} c_i \sum_{k=2}^{\infty}\frac{k \nu_i^{k-2}}{5^k} \\ = -\frac{9}{125} + \sum_i \frac{c_i(10-\nu_i)}{5(5-\nu_i)^2} = -0.072 + 0.378409 + 2\cdot 0.00289355 \\ = 0.3121961.
$$
For comparison, in a simple simulation of $10^8$ trials (running each trial for at most $1000$ jumps), $31218647$ asymmetric jumpers returned safely to the origin. This gives an estimate of $0.31219 \pm 0.00006$ for the probability; i.e., the simulation agrees perfectly with the matrix calculation above.
Finally, for $r=3$ we have the following set of states:
Several pairs of states are related by red-blue exchange symmetry (e.g., $B$ and $B'$). Using this symmetry, we can write the transfer matrix in terms of symmetric states only (e.g., $B/B'$), as described earlier for the $r=2$ case. The transitions out of unprimed states are:
- $A$ to $(B,B',C,C',D,H(\times 3))$
- $B$ to $(B',C,C',E,H(\times 2))$
- $C$ to $(B',H)$
- $D$ to $(F,F',G,G')$
- $E$ to $(F',G')$
- $F$ to $G$
- $G$ to $(C,H)$.
So
$$
N_k^{(3)} = \left(\begin{matrix}0 & 0 & 0 & 0 & 0 & 0 & 1 \end{matrix}\right)\left(\begin{matrix}
1 & 1 & 0 & 0 & 0 & 0 & 0\\
2 & 0 & 0 & 0 & 0 & 1 & 0\\
0 & 0 & 0 & 0 & 0 & 0 & 0\\
1 & 0 & 0 & 0 & 0 & 0 & 0\\
0 & 0 & 2 & 1 & 0 & 0 & 0\\
0 & 0 & 2 & 1 & 1 & 0 & 0\\
2 & 1 & 0 & 0 & 0 & 1 & 0\\
\end{matrix}\right)^{k-2}
\left(\begin{matrix}
2 \\ 2 \\ 1 \\ 0 \\ 0 \\ 0 \\ 3
\end{matrix}\right).
$$
This sequence begins with $3,6,14,30,62,130,\ldots$. After diagonalizing the matrix, we arrive at
$$
N_k^{(3)}=\delta_{k,2}+\left(\begin{matrix}1 & 1 & 1 & 1 \end{matrix}\right)\left(\begin{matrix}
\nu_1 & 0 & 0 & 0 \\
0 & \nu_2 & 0 & 0 \\
0 & 0 & \nu_3 & 0\\
0 & 0 & 0 & \nu_3^{*}
\end{matrix}\right)^{k-2}\left(\begin{matrix}
c_1 \\ c_2 \\ c_3 \\ c_3^{*}
\end{matrix}\right)=\delta_{k,2} + \sum_i c_i \nu_i^{k-2},
$$
where
$$
\nu_1 = -0.716673 \\
\nu_2 = 2.10692 \\
\nu_3 = 0.304877 - 0.754529 i \\
c_1 = -0.188335 \\
c_2 = 3.136179 \\
c_3 = -0.4739206 - 0.3006352 i.
$$
Finally, the desired probability is
$$
P^{(3)} = \sum_{k=2}^{\infty}\frac{kN_k^{(3)}}{6^k}=\frac{1}{18}+\sum_i \frac{c_i(12-\nu_i)}{6(6-\nu_i)^2}=0.325873.
$$
Checking by simulation for $10^8$ trials yields an estimate of $0.32587 \pm 0.00006$, wholly consistent with this result.
Automating the process
One can do all of the following by making a number of formalizations.
It is advantageous to mechanize the process before proceeding to larger cases; in particular, we will represent a state as a triple $(X,Y,P)$ where $X$ is the set of remaining ranges of rightward (red) edges, $Y$ is the set of leftward (blue) edges, and $P$ is a partition of $X\cup Y$ representing which edges have been connected (noting that if $x=y$ then the edges represented by $x\in X$ and $y\in Y$ originated together and are clearly connected). We will say $(X',Y',P')$ is an $n$-translate of $(X,Y,P)$ if we might see the former after moving our "window" $n$ steps to the right. For convenience, define
$$m=\min(X\cup Y)$$
$$S-n=\{s-n:s\in S\}$$
and, because the relations are not conducive to standard notation, I will leave to the reader the problem of relating $P'$ and $P$. Then, we can have the following situations.
Birth: In the state where the node in question has two rightward edges passing from it. We require:
$$X'=X-n\cup \{r\}$$
$$Y'=Y-n\cup \{r\}$$
We will consider $r$ to be its own connectedness class in $P'$, and otherwise just shift everything.
Right-Continuation: This is the state where the node receives a rightward edges and also outputs one, we need that for some $c\in X-n$ we have:
$$X'=X-n\cup\{r\}\setminus\{c\}$$
$$Y'=Y-n$$
The partition will essentially shift and add $r$ to the equivalence class of $c$.
Left-Continuation: This is where a node receives and produces a leftward edge. Flop $X$ and $Y$ in the above definition.
Death: This is where a node is a local maximum, receiving two edges. We need that
$$X'=X-n\setminus\{c_x\}$$
$$Y'=Y-n\setminus\{c_y\}$$
where $c_x$ and $c_y$ are appropriate members of $X-n$ and $Y-n$ respectively. Neither $X'$ nor $Y'$ may contain $0$. The equivalence classes representing $c_x$ and $c_y$ will be merged. If $c_x$ and $c_y$ were the last remaining representatives of their equivalence class, then the remaining edges were not connected to them; this transition is invalid unless $X'=Y'=\{\}$.
Then, in the matrix, the number of ways to transition from one state to another is the number of $n$'s such that the latter is an $n$-translate of the former. My code works by generating every conceivable state - i.e. every pair $(X,Y,P$) where $|X|=|Y|$ and $P$ is a partition of $X\cup Y$. Some of these are inaccessible from the starting state or cannot access the ending state, but I did not write code to prune such states away. I also did not take advantage of the symmetry $(X,Y,P)\rightarrow(Y,X,P)$. Either would likely lead to big increases in efficiency, allowing the computation of more values $P^{(r)}$.
It is worthy of note that
$$\sum_{k=2}^{\infty}a\cdot \frac{kM^{k-1}}{2k}\cdot b=a\cdot\left(\sum_{k=2}^{\infty}\frac{kM^{k-1}}{2r}\right)\cdot b$$
by linearity, meaning we can first use the identity that
$$\sum_{k=2}^{\infty}\frac{kM^{k-1}}{(2r)^k}=\frac{1}{4r^2}(2I-\frac{1}{2r}M)\cdot M\cdot(\frac{1}{2r}M-I)^{-2}$$
which is easily derived by noting that $M$ acts, for these purposes, just like a real number with $|M|<2r$ as $M$ must have no eigenvalues of absolute value $2r$ or more (given that there are only $(2r)^k$ paths of length $k$, whereas there are about $k\lambda^k$ cycles of length $k$ with identified starting point where $\lambda$ is the largest eigenvalue of $M$). Then, we can take the transfer matrix, directly compute the infinite sum as above, and extract the desired coefficient. It should be noted that this will always yield a rational answer. In particular, we receive:
$$P^{(1)}=\frac{1}2$$
$$P^{(2)}=\frac{7}{18}$$
$$P^{(3)}=\frac{4368595}{13405842}\approx 0.325872$$
$$P^{(4)}=\frac{33373525946827013707062414432976}{117177232862732965149684560993569}\approx 0.284812$$
Mathematica code may be found in this revision. It may be noted that the methods may be modified easily to handle asymmetric cases or cases with non-uniform probability - essentially, instead of using our transfer matrix as a count, we can use it to measure the probability of a cycle occurring given that we start on a given point of it.
Best Answer
The abstract general recipe looks like this. Given a domain $D$, two sets $A,B \subset D$ and a Markov process with generator $L$, the probability to hit $A$ before hitting $B$ starting from $x$ solves the following problem
\begin{aligned} (Lu)(x) = 0 & \quad x \not \in A \cup B \\ u(x) = 1, & \quad x \in A \\ u(x) = 0, & \quad x \in B. \end{aligned}
As is this is not the problem that you want to solve, because your problem has no $B$. To get around this, you can solve a sequence of problems where you "move $B$ to infinity". In the case of the symmetric random walk on $\mathbb{Z}$ where you are interested in hitting zero, you can take $A=\{ 0 \},L=P-I$ where $P$ is the transition probability matrix, and $B=\{ n \}$, where $n$ is a large say positive integer. (Since the walk is not symmetric, the positivity assumption changes matters: the situation is different depending on whether the initial jump is to the left or to the right.) The result is a linear system of equations, which in the positive case looks something like
$$u_0 = 1 \\ (1-p)u_0 - u_1 + pu_2 = 0 \\ (1-p)u_1 - u_2 + pu_3 = 0 \\ \dots \\ (1-p)u_{n-2} - u_{n-1} + pu_n = 0 \\ u_n = 0.$$
This can be solved by recurrence relation techniques; for $p \neq 1/2$ the general solution is a combination of exponentials, to which you can adjoin the boundary conditions. In this case when you send $n \to \infty$, you will find that for all $k>0$, $u_k$ converges to zero. On the other hand, if you start on the left and keep the process biased towards moving to the right, then $u_k$ will converge to $1$ instead.
Note that a closely related recipe works for finding the expected time to hit $0$ starting from a point $x$.