help-octave
[Top][All Lists]

## 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-----
Sent: Saturday, December 11, 2004 1:12 PM
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
-------------------------------------------------------------