I need to obtain a closed form of $M(t)$, satisfying the following recurrence equation:
$$M(t+1)=a+bM(t)+\frac{c}{t+1}\sum_{t'=0}^tM(t')+df(t)$$
Where $f(t)$ is a known function and $a$, $b$, $c$ and $d$ are known constants. And with the initial condition $M(t=0)=0$.
I've never solved an equation of this kind, so I'm asking for a known analythic method for obtaining a closed form of the function $M(t)$. Maybe a reference book may help, or just the name of one method of solution.
On
I have an explicit solution, which I obtained by using the method of generating functions. First, compute what I called $P_n(1)$ from a recurrence relation. Then use those to compute explicit values of $M(k)$. $P_n(1)$ is the sum of the $M(k)$ values, for $0 \le k \le n$. You don't even need to compute all $P_n(1)$ values ahead of time. You can compute the $P$'s and the $M$'s in alternation (one $P$, one $M$, one $P$, one $M$, etc).
$ \begin{array}{rcl} P_n(1) &\!\!\!=&\!\!\! \sum_{k\,=\,1}^{n} b^k \sum_{j\,=\,0}^{k-1} \frac{\displaystyle (j+1)\,[\,a + d\,f(j)\,] + c\,P_j(1)}{\displaystyle (j+1)\,b^{j+1}} \qquad (n \ge 1) \\[0.2in] % M(k) &\!\!\!=&\!\!\! b^k \sum_{j\,=\,0}^{k-1} \frac{\displaystyle (j+1)\,[\,a + d\,f(j)\,] + c\,P_j(1)}{\displaystyle (j+1)\,b^{j+1}} \qquad (k>0)\,. \end{array} $
Now, the proof of this result took me 4 pages of LaTeX, with lots of macros of my own that are not supported by the LaTeX engine here so I'm not sure what the best way to share that proof is. For now, I'm going to post screen shots.
There is, I fear, not much to be done for lower values of $t$. But the asymptotics for $t\to\infty$ could reasonably be estimated, depending on the values of the parameters. For instance, if $b>0$, write $h(t)=tb^{-t}M(t)$ and consider the functional equation $$h(t+1)=\frac{t+1}th(t)+(t+1)b^{-t-1}\left(a+df(t)\right)+\frac cb \frac {t+1}{t}\int_0^tb^{u-t}h(u)\frac{\mathrm du}u$$ which is equivalent to your recurrence equation. As $t\to\infty$, one can expand $h(t+1)$ to the first order and neglect the second order terms in $t$ (i.e. using $t+1\simeq t$) and get a more addressable equation $$h'(t)=tb^{-t}\left(a+df(t)\right)+\frac cb\int_0^tb^{u-t}h(u)\frac{\mathrm du}u.$$
From there, everything depends on the function $t\mapsto tb^{-t}\left(a+df(t)\right)$. It is impossible to say something general and you will probably have to distinguish cases, that is find which term is leading the asymptotic behaviour.
However, in the case the function $t\mapsto tb^{-t}(a+df(t))$ is bounded you can use Grönwall's inequality to bound the solution. Let us say that $tb^{-t}(a+df(t))$ is between $A>0$ and $B>0$ for sufficiently large $t$. Since the antiderivative of $u\mapsto b^u/u$ is $u\mapsto\mathrm{Ei}\left(u\log b\right)$, you get the asymptotic behaviour of $h$ as roughly $h(t)\sim \exp\left(\frac cb \mathrm {Ei}(t \log b)\right)$ which gives $$M(t)\sim \frac{b^t}t\exp\left[\frac cb\mathrm{Ei}\left(t\log b\right)\right].$$ The leading constant would be in that case a number between $A$ and $B$.