Re: suite

[Paul Kienzle]

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
-------------------------------------------------------------