1
$\begingroup$

I posted a similar question earlier today, but unfortunately I realized after several hours that I over-simplified the problem for this forum. There is a very good answer to the Computing a sequence by recurrence question. I just can't do the simplification step that turns the sums into $Q_{N-1}$'s. This is because the form of my problem is not as simple as the question I posed.

I am trying to create a list where each element of the list is calculated using the sum of all previous elements of the list and also uses the $N$th element of another list. The initial value of the list is a known quantity. $\alpha$, $\beta$, $\gamma$, $\delta$, $\phi$, and $\theta$ are all known constants $k_N$ is a known list of $N$ elements.

$\quad \quad Q_N=\alpha\frac{\beta-\sum_{i=0}^{N-1}{Q_i\delta}}{k_N}\sqrt{(\frac{k_N}{\gamma}\frac{\beta+\sum_{i=0}^{N-1}{Q_i\delta}}{\beta-\sum_{i=0}^{N-1}{Q_i\delta}})^{\phi}-(\frac{k_N}{\gamma}\frac{\beta+\sum_{i=0}^{N-1}{Q_i\delta}}{\beta-\sum_{i=0}^{N-1}{Q_i\delta}})^{\theta}},\,N>0$ $\quad \quad Q_{N=0}=Q_0$

I'm not quite sure how to code this. I think the Nest function is somewhat close to what I want, but I need to keep plugging in the sum of all previous values, not just the previous value.

FoldList and Nest seem to be close to what I'm trying to do, but I can't figure it out. It may be relevant that $N$ is going to be rather large. The context is a numeric solution to the time evolution of a system. So I want to have smaller time-steps, which leads to a large number of elements.

Any help would be greatly appreciated.

Note: I'm not asking the forum to figure out this rather complicated form necessarily; I just wanted to show the form to say that the simplification made to use FoldList in the other answer won't work here. I think the computation requires using all previous values, and not just the previous value. Answering the other question just without using the simplification may suffice.

Computing a sequence by recurrence

$\endgroup$
  • $\begingroup$ How long is your sequence, ie, how large of a value of $N$ are you looking for? Note that defining $R_N=\sum_{i=0}^{N-1}Q_i$, one has $R_{N+1}=R_N+Q_N$, and the expression for $Q_N$ only depends on $R_N$. So as long as $N$ is a couple million or less, you can just compute both $Q_N$ and $R_N$ at each stage, store the results, and then compute $Q_{N+1}$. $\endgroup$ – DumpsterDoofus Jan 9 '15 at 23:32
  • $\begingroup$ If you're computing as large as $N$ in the tens of millions or billions, then you may have to think more carefully, but the above method is efficient in the sense that each new term is computed in $O(1)$ time, and thus is probably asymptotically optimal. $\endgroup$ – DumpsterDoofus Jan 9 '15 at 23:33
  • $\begingroup$ I think an upper bound of 1 million for N. $\endgroup$ – matt Jan 10 '15 at 0:54
  • $\begingroup$ In that case, the method I used in my answer below should work fine. My example take 9 seconds for 1 million terms, so your slightly more complicated example should probably take 1 minute or less to compute. $\endgroup$ – DumpsterDoofus Jan 10 '15 at 1:05
4
$\begingroup$

This is the kind of problem which is most efficiently solved via memoization. As I mentioned in a comment, defining $R_N=\sum_{i=0}^{N-1}Q_i$, we have

$$\begin{align} Q_N&=f(R_N,k_n)\\ R_N&=R_{N-1}+Q_{N-1} \end{align}$$

where $f$ is some function which contains the specifics. In your case,

$$f(R_N,k_N)=\alpha\frac{\beta-\delta R_N}{k_N}\sqrt{\left(\frac{k_N}{\gamma}\frac{\beta+\delta R_N}{\beta-\delta R_N}\right)^{\phi}-\left(\frac{k_N}{\gamma}\frac{\beta+\delta R_N}{\beta-\delta R_N}\right)^{\theta}}$$

but I'll use the simpler example

$$\begin{align} f(R_N,k_N)& =\sin\left(R_N+k_N\right)\\ k_N& =N. \end{align}$$ The following implements the memoization:

Clear[Q, R]
Q[0] = R[0] = 0;
Q[n_] := Q[n] = Sin[1.0 R[n] + n]
R[n_] := R[n] = R[n - 1] + Q[n - 1]

The following computes the first 100,000 values of $Q_N$ and stores their timings in a list, which is then plotted (after smoothing to remove jitter):

times = Table[First@AbsoluteTiming@Q[i], {i, 100000}];
ListLinePlot[GaussianFilter[times, 1000], PlotRange -> All]

enter image description here

As you can see, the performance is $O(1)$ and is thus asymptotically optimal. On my computer, computing 100,000 values takes 0.9 seconds, so you can compute the first 1,000,000,000 values of $Q_N$ in 2.5 hours.

It is straightforward to modify the code to work on your example (you just have to change the definition of Q[n]). If you need more than a billion values of $Q_N$, you may want to consider using a different programming language, or ask a separate question on how to performance-tune memoization.

| improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.