I have the following discrete time dynamical system
$$ y(t+1) = y(t) + \frac{1}{1+\exp(z+ u f y(t))} ,\quad y(0)=0,$$
where $z$ is a real number $f$ and $u$ are non-negative reals. I know I have little hope of obtaining a closed form solution for this process. But, actually, for my application, a "better" solution involves finding (making the dependence of $y$ on $f$ explicit by writing $y(t)$ as $y(t,f)$): 
$$\lim\_{n \to \infty} \frac{1}{n}\cdot y\left(nt,\frac{f}{n}\right).$$
**Update:** Initial simulations suggested that this converged to $$y(t) = \frac{t}{1+\exp(z)}.$$ But that was wrong. In conclusion, the differential equation is indeed a good approximation to the above limit. Sorry about the confusion.


How can I rigorously show this? Also appreciated are references to texts that discuss similar problems.


Thanks



Limit of a discrete time dynamical system
The continuous model of the problem suggests that the limit does depend on $f$ (and $u$). More precisely, it depends on how fast the parameter $f$ is suppressed in the expression whose limit you are taking; behavior of $\lim y(nt, f\_n)/n$ will depend on the limit of $nf\_n$ as $n \to \infty$. The answer will be a function of the limit of $nuf\_n$. Only when this limit is zero does one get the proposed formula.


The associated differential equation is $Y' = 1/(1+Ae^{BY})$ where $A = e^z$ and $B=uf$. Its solution vanishing at 0 is $Y(t) = H^{-1}(t)$ where $H(t) = t + (A/B)(e^{Bt} - 1)$. It does look like this matches the asymptotic behavior of your sequence for $y$ both in the large and small range. For small $t$, the expansion $Y(t) = t/(1+A) + O(t^2)$ corresponds to your formula, but I think the answer is not quite that simple: you have to establish whether the expression whose limit is calculated belongs to the small regime where $Y(t)$ is approximately linear, or the large regime where $Y(t)$ is logarithmic, $Y(t)=O(\log(t))$. The limit uses $n$ iterations so we want to know, as a function of $B \sim 1/n$, whether the transition between regimes happens at a point much larger than $1/B$. However, it's easy to calculate that the ratio $H(t)/t$ moves away from 1 (the difference is larger than some constant independent of $B$) as soon as $Bt$ is of order 1 (i.e., bounded below by a given positive constant) and this would spoil the limit if the differential equation is a good model of the difference equation. 


(ADDED: for comparison of $Y$ predictions with $y$ simulations, in the phase transition where $Bt$ is of order 1, $Y(t) \sim t/C$ and $H(t) \sim Ct$, with $C = 1 + A(e^q - 1)/q$, and $q = Bt =uft$. That is, $Y$ stays approximately linear but the coefficient goes to zero, consistent with the idea that it's turning into a logarithmic function. Let $t=nt\_0, \quad f=f\_0/n$, for some constant $f\_0$ and with $u$ and $t\_0$ also held constant while $f$ varies with $n$, so that the phase transition parameter is $q=uf\_0 t\_0$ and the predicted value of the limit, if $Y$ is a good approximation for $y$, is $L\_{pred} = \lim Y(nt\_0,f\_0/n)/n = \lim nt\_0/nC = t\_0/C = t\_0(q/(q + Ae^q - A))$. In the original notation of the question, $L = t/(1+{e^z}F(uft))$ where $F(x)=(e^x-1)/x$. Does this match the simulations?)


To see the small-$y$ behavior directly in the difference equation, it can be expanded in powers of $y$. 


$y(t+1) - y(t) = 1/(1+A) - (AB)/(1+A)^2)y + O(y^2)$


Your formula proposes that when $B \sim 1/n$, the effect of the $y^{\geq 1}$ terms is of order smaller than $n$ for $t \in [0,n]$. The sum of the first $t$ values of the $y^1$ term will be of order $t^2$, so one expects these corrections to be suppressed only on a short interval, $t << n^{1/2}$. The calculation with the differential equation suggests that $f\_n = f/n$ is too large a parameter ; this calculation with the truncated difference equation can be used to prove that $f\_n = f/n^k$ is small enough for any $k > 2$. Adding higher degree terms to the approximate difference equation would, I suppose, only get closer to the picture suggested by the differential equation. 


To prove rigorously the predictions from the differential equation you could try to control $y$ by trapping the sequence $y(n)$ between two trajectories of the ODE. If simulations are consistent with a heuristically "wrong" formula it would be very interesting to sort out what the truth is.


