Closed form formula instead of recursive sequence I'm trying to create a computational model for a neuroscience project, but the computation times are too long for it to be useful. In particular, there is an iterative recursive step that is too slow (because I go through it millions of times). It would greatly help me if I have a closed form expression to compute this step directly. However, I cannot figure out how to obtain it.
The iterative procedure is as follows:
at $t=1$, $V_t$ equals $S+P$,
at $t>1$, $V_t$ equals $S+P-S^2(V_{t-1} + S)^{-1}$,
where,
$S$ and $P$ are positive constants, and $V_{t-1}$ equals $V_t$ at the previous timepoint.
Timepoints t are integers $\geq 1$.
Is there a closed form formula with which I can directly compute $V_t$ for any given $t$? An approximation, fairly precise for $t<100$, would also be very useful!
 A: $$
V_t = S + P - \frac{S^2}{V_{t-1} + S}\\
= \frac{(S+P)V_{t-1} + (S+P)S - S^2}{V_{t-1}+S}\\
= \frac{(S+P)V_{t-1} + PS}{V_{t-1}+S}\\
V_{t+1} = \frac{(S+P)V_t + PS}{V_t+S}\\
\alpha = \frac{2S+P}{1} = 2S+P\\
\beta = \frac{(S+P)S-PS}{1^2} = S^2\\
x_{t+2} - (2S+P) x_{t+1} + (S^2) x_t = 0\\
x_t \equiv C r^t\\
r^2 - (2S+P) r + S^2 = 0\\
r_{\pm} = \frac{(2S+P) \pm \sqrt{(2S+P)^2 - 4 S^2}}{2}\\
= \frac{(2S+P) \pm \sqrt{P^2 + 4SP}}{2}\\
x_t = C_+ r_+^t + C_- r_-^t\\
\frac{x_2}{x_1} = y_1 = V_1 + S = 2S+P\\
$$
WLOG take $x_1=1$. So we have two equations to solve to get $C_\pm$.
$$
C_+ r_+ + C_- r_- = x_1 = 1\\
C_1 r_+^2 + C_- r_-^2 = x_2 = 2S+P\\
C_1 r_+^2 + C_- r_+ r_- = r_+\\
C_- (r_-^2 - r_+ r_-) = 2S+P-r_+\\
C_- = \frac{2S+P-r_+}{r_-^2 - r_+ r_-}\\
C_+ = \frac{1-C_- r_-}{r_+}\\
$$
Now from that to $V_t$ we say
$$
V_t = y_t - S = \frac{x_{t+1}}{x_t} - S\\
= \frac{C_+ r_+^{t+1} + C_- r_-^{t+1}}{C_+ r_+^t + C_- r_-^t} - S
$$
from math import sqrt
import typing
S=5
P=1
def Vt(t:int,S=S,P=P) -> float:
    negb = (2*S+P)
    sqrt_term = sqrt(P**2+4*S*P)
    rplus = (negb + sqrt_term)/2
    rminus = (negb - sqrt_term)/2
    Cminus = (negb - rplus)/(rminus**2-rplus*rminus)
    Cplus = (1-Cminus*rminus)/rplus
    xt = Cminus*rminus**t + Cplus*rplus**t
    xtplus1 = Cminus*rminus**(t+1) + Cplus*rplus**(t+1)
    return xtplus1/xt - S
print(f"V_1 = {Vt(1)}, should be {S+P}")
actual = Vt(1)
for i in range(2,10):
    expected = S+P-S**2/(actual+S)
    actual=Vt(i)
    print(f"V_{i} = {actual}, should be {expected} and that is {actual-expected} off")

