[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
## Re: suite

**From**: |
Paul Kienzle |

**Subject**: |
Re: suite |

**Date**: |
Sat, 11 Dec 2004 21:00:14 -0500 |

Extending this just a little bit (after a minor correction),
U_n = a^n x + b \sum_{i=0}^{n-1}{ a^i }
= a^n x + b (a^n - 1)/(a-1)
= (a^n ( x (a-1) + b) - b) / (a-1)
which in Octave is:
(a.^n * ( x * (a-1) + b) - b ) / (a-1)
Some times:
octave> n=1000; x = 4; a = 0.3; b = 2;

` octave> tic; z = zeros(1,n+1); z(1)=x; for i=2:n+1; z(i)=a*z(i-1)+b;
``end; toc
` ans = 0.29167
octave:26> tic; z2 = (a.^[0:n] * ( x * (a-1) + b) - b ) / (a-1); toc
ans = 0.016500
octave> norm(z-z2)
ans = 2.0830e-15
Just getting U_1000 is not much cheaper:
octave:28> tic; U_1000 = (a^n * ( x * (a-1) + b) - b ) / (a-1); toc
ans = 0.014853
- Paul
On Dec 11, 2004, at 3:37 PM, Hall, Benjamin wrote:

If you write out some of the terms you get
U_0 = x
U_1 = a * x + b
U_2 = a ( a*x + b ) + b
U_3 = a ( a ( a*x + b ) + b ) + b = a^3*x + a^2b + ab + b
U_n = a^(n-1)*x + a^(n-2)*b + a^(n-3)*b + ... ab + b
so you could calculate U_n as follows
U_n = a^(n-1)*U_1 + sum( cumprod( repmat(a, [1 n-2] ) ) * b )

`It's harder to read than Mike's suggestion and you only get U
``evaulated for
``a single case plus its is only marginally faster than the for-loop
``below on
`my machine:
0.66 seconds for n=200 repeated 20 times for "for loop"
0.55 seconds for formula above with n=200 and repeated 20 times
So much for the fancy formula...
Ben
-----Original Message-----
From: Mike Miller [mailto:address@hidden
Sent: Saturday, December 11, 2004 1:12 PM
To: adel.essafi
Cc: help
Subject: Re: suite
On Sat, 11 Dec 2004, adel.essafi wrote:

hi list
please, how can I define a suite in octave
U(n+1)=3*U(n)+2 (for example)
and calculate u(1000)

A simple method:
n=1000;
U=zeros(1,n);
U0=4; #or whatever you want
U(1)=3*U0+2;
for i=2:n,
U(i)=3*U(i-1)+2;
end
That would work, but there are probably much nicer ways to do this and
other people will have more sophisticated ideas than mine!
Mike

-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.
Octave's home on the web: http://www.octave.org
How to fund new projects: http://www.octave.org/funding.html
Subscription information: http://www.octave.org/archive.html
-------------------------------------------------------------

**suite**, *adel\.essafi*, `2004/12/11`
**RE: suite**, *Hall, Benjamin*, `2004/12/11`
**Re: suite**,
*Paul Kienzle* **<=**