Closed form formula instead of recursive sequence

Question

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!

Answer 1

$$ 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")