% A geometric series is a sum in which each
% term is obtained from the previous term by
% multiplying by the same value.
%
% For example,
%
%              2    3    4
%     1 + z + z  + z  + z  + ...
%
% is a geometric series with common ratio z.
%
% There is no requirement that the common ratio
% be real. The value z may be any complex
% number.


% You are welcome to change this value of 'z'.

z = 0.2 + 0.5j

fprintf('\nThe common ratio z:\n')
z

fprintf('\nThe magnitude of z:\n')
abs(z)


% First, consider the finite geometric sum
%
%          N
%        -----
%         \     n
%          )   z
%         /
%        -----
%        n = 0
%
% which means
%
%                    2              N
%     1 + z + z  + ... + z .
%
% Let this sum be S:
%
%                         2              N
%     S = 1 + z + z  + ... + z .
%
% Multiplying both sides by z gives
%
%                              2              N+1
%     z*S = z + z  + ... + z  + z.
%
% Subtracting the second equation from the first
% causes almost every term to cancel:
%
%                         N+1
%     S - z*S = 1 - z.
%
% Factoring the left side gives
%
%                         N+1
%     S*(1 - z) = 1 - z.
%
% Therefore, provided z is not equal to 1,
%
%             N+1
%        1 - z
%     S = --------.
%         1 - z
%
% Thus,
%
%          N
%        -----            N+1
%         \     n    1 - z
%          )   z  = --------.
%         /          1 - z
%        -----
%        n = 0
%
% This formula works for complex values of z
% just as it does for real values of z.


% Let us calculate the finite sum for N = 10.

N = 10;

sum = 0;

for n = 0:N
    sum = sum + z^n;
end

fprintf('\nFinite sum up to %d:\n', N)
fprintf('  %.16f %+.16fj\n', real(sum), imag(sum))


% We can calculate the same finite sum directly
% from the geometric-sum formula.

finite_formula = (1 - z^(N + 1))/(1 - z);

fprintf('\nUsing the finite-sum formula:\n')
fprintf('  %.16f %+.16fj\n', ...
        real(finite_formula), imag(finite_formula))


% These two calculations should give essentially
% the same result. Any very small difference is
% the result of floating-point rounding.


% We are especially interested in what happens
% as N becomes larger and larger.
%
% The finite sum is
%
%             N+1
%        1 - z
%     S = --------.
%         1 - z
%
% If
%
%     |z| < 1,
%
% then
%
%      N+1
%     z     -> 0
%
% as N becomes arbitrarily large.
%
% Why?
%
% Recall that
%
%      N+1       N+1
%     |z    | = |z|   .
%
% If |z| < 1, repeatedly multiplying |z| by
% itself produces values that become smaller
% and smaller, approaching zero.
%
% Consequently,
%
%             N+1
%        1 - z            1
%     -----------  ->  -------.
%          1 - z        1 - z
%
% Therefore,
%
%        infinity
%        --------
%           \       n       1
%            )     z  = --------
%           /           1 - z
%        --------
%          n = 0
%
% provided
%
%     |z| < 1.


% Check whether the infinite geometric series
% converges for our chosen value of z.

if abs(z) >= 1

    fprintf(['\nFor the infinite geometric series ' ...
             'to converge, |z| < 1.\n'])

else

    infinite_sum = 1/(1 - z);

    fprintf('\nThe infinite series converges to:\n')
    fprintf('  %.16f %+.16fj\n', ...
            real(infinite_sum), imag(infinite_sum))

end


% Let us now see numerically how the finite sums
% approach the infinite sum.
%
% Start by calculating
%
%          10
%        -----
%         \     n
%          )   z .
%         /
%        -----
%        n = 0

sum = 0;

for n = 0:10
    sum = sum + z^n;
end

fprintf('\nFinite sum up to 10:\n')
fprintf('  %.16f %+.16fj\n', real(sum), imag(sum))


% Add the next ten terms so that we now have
%
%          20
%        -----
%         \     n
%          )   z .
%         /
%        -----
%        n = 0

for n = 11:20
    sum = sum + z^n;
end

fprintf('\nFinite sum up to 20:\n')
fprintf('  %.16f %+.16fj\n', real(sum), imag(sum))


% Add the next eighty terms so that we now have
%
%         100
%        -----
%         \     n
%          )   z .
%         /
%        -----
%        n = 0

for n = 21:100
    sum = sum + z^n;
end

fprintf('\nFinite sum up to 100:\n')
fprintf('  %.16f %+.16fj\n', real(sum), imag(sum))


% Compare this with the exact value of the
% infinite geometric series.

if abs(z) < 1

    infinite_sum = 1/(1 - z);

    fprintf('\nInfinite sum:\n')
    fprintf('  %.16f %+.16fj\n', ...
            real(infinite_sum), imag(infinite_sum))

end


% We can do more than simply observe that the
% finite sums appear to be approaching the
% infinite sum.
%
% We know exactly what the difference must be.
%
% The infinite sum is
%
%             1
%     S = --------
%          1 - z
%
% while the finite sum through z^N is
%
%                  N+1
%             1 - z
%     S_N = -----------
%               1 - z.
%
% Subtracting gives
%
%                       N+1
%                      z
%     S - S_N = -------------.
%                     1 - z
%
% Thus, the error remaining after adding the
% terms through z^N is exactly
%
%          N+1
%         z
%     ---------.
%       1 - z
%
% Taking magnitudes gives
%
%                         N+1
%                      |z|
%     |S - S_N| = -------------.
%                       |1 - z|
%
% Since |z| < 1, the numerator becomes smaller
% and smaller as N increases.


% Let us compare the calculated error with this
% formula.

if abs(z) < 1

    infinite_sum = 1/(1 - z);

    fprintf('\nErrors in the finite sums:\n')

    for N = [10 20 100]

        finite_sum = (1 - z^(N + 1))/(1 - z);

        error = abs(infinite_sum - finite_sum);

        predicted_error = ...
            abs(z)^(N + 1)/abs(1 - z);

        fprintf('\nN = %d\n', N)
        fprintf('Calculated error: %.16e\n', error)
        fprintf('Predicted error:  %.16e\n', ...
                predicted_error)

    end

end


% There is also a useful geometric
% interpretation of what is happening.
%
% Write the common ratio in polar form:
%
%     z = r exp(j*theta).
%
% Then
%
%      n    n
%     z  = r exp(j*n*theta).
%
% Each successive term therefore does two
% things:
%
%   1. Its magnitude is multiplied by r = |z|.
%
%   2. Its argument is increased by theta.
%
% If |z| < 1, the magnitudes of the terms
% become smaller and smaller.
%
% If the argument of z is nonzero, the
% direction of each successive term also
% rotates.
%
% Thus, the terms can be thought of as vectors
% that rotate while becoming progressively
% shorter.


% For our example,

fprintf('\nMagnitude of the common ratio:\n')
abs(z)

fprintf('\nArgument of the common ratio in radians:\n')
angle(z)

fprintf('\nArgument of the common ratio in degrees:\n')
rad2deg(angle(z))


% Each successive term has its magnitude
% multiplied by |z| and its argument increased
% by angle(z).
%
% Consequently, the partial sums generally do
% not approach the infinite sum along a straight
% line. Instead, they may spiral toward it.


% Finally, consider what happens if |z| >= 1.
%
% If |z| > 1, then the magnitudes
%
%          2     3
%     1, |z|, |z|, |z|, ...
%
% grow rather than approach zero.
%
% If |z| = 1, then
%
%      n
%     |z| = 1
%
% for every n, so the terms again do not
% approach zero.
%
% A necessary condition for any infinite series
% to converge is that its individual terms
% approach zero.
%
% Therefore, the infinite geometric series
%
%        infinity
%        --------
%           \       n
%            )     z
%           /
%        --------
%          n = 0
%
% converges exactly when
%
%     |z| < 1.
%
% When it converges, its value is
%
%          1
%     ---------.
%       1 - z