MATLAB/Divergent series investigations

Shown below are the 15 terms of the series for ε equal to 0, .001, and .01, respectively. The series converges to -1/12 in the limit that ε>0 vanishes.

S(ε=.000) = 1.0 + 2.0 + 3.0 + 4.0 + 5.0 + 6.0 + 7.0 + 8.0 + 9.0 + 10.0 + 11.0 + 12.0 + 13.0 + 14.0 + 15.0 + ...
S(ε=.001) = 1.0 + 2.0 + 3.0 + 4.0 + 5.0 + 6.0 + 7.0 + 7.9 + 8.9 + 09.9 + 10.9 + 11.9 + 12.8 + 13.8 + 14.8 + ...
S(ε=.010) = 1.0 + 2.0 + 2.9 + 3.8 + 4.8 + 5.6 + 6.5 + 7.4 + 8.2 + 09.0 + 09.8 + 10.6 + 11.3 + 12.1 + 12.8 + ...	

This code runs on both MATLAB and Octave. It sums,

${\displaystyle S(N,\epsilon )=\sum _{n=1}^{N}nz^{n}}$ where, ${\displaystyle z=e^{-\epsilon +i\epsilon }\approx 1\ \mathrm {if} \ \epsilon <<1}$.

The numerical simulation suggests that for real ${\displaystyle \epsilon }$:

${\displaystyle \lim _{\epsilon \to 0}\sum _{n=1}^{\infty }ne^{-\epsilon n}\cos(\epsilon n)=-{\frac {1}{12}}}$

Warning: This result does not hold if ${\displaystyle \epsilon }$ is a complex number that goes to zero, unless the real and imaginary parts are constrained to be equal in magnitude. This would be a significant simulation if the series converges to -1/12 for all ratios of real to imaginary parts of ε. Nevertheless, it is curious that a real to imaginary ratio of one yields the well known value, 1/12≈.083333

1. Verify that ${\displaystyle \int _{0}^{\infty }xe^{-\epsilon x}\cos(\epsilon x)dx=0}$ if ε > 0
2. Why is this interesting?
3. Can Octave do symbolic integration?

The code that follows gave this output:

If z= 
0.990000332+0.009900333i the sum from 1 to 4000 of n*z^n is 
 -8.333333332e-02 + i 4.999999999e+03

% This code is designed to investigate the "astonishing" claim that  
%   
% 1 + 2 + 3 + 4 + 5 + ... = -1/12  
%  
% It numerically calculates  
%  
%   1*z + 2*z^2 + 3*z^3 + 4*z^4 + 5*z^5 + ... N*z^N 
%   
% where N is large and z is nearly equal to, 1.0 + 0.0i  
clear;clc; close all;% clear memory and command window  
N=4000; % sets the number of terms in the expansion
shortN=20; % how far we plot to show initial behaviour of the series
alpha=.01; %the real part of z
beta=.01;
epsilon = -alpha+1*i*beta;
% The result is astonishing if alpha=beta
z=exp(epsilon);
theta=atan(abs(imag(z)/real(z)));%is the argument of z
period=round(2*pi/theta);
%We estimate the limit as N->infty by averaging over one period
n=1:N; % the partial sums shall be plotted against these integers  
seq=linspace(0,0,N);  % preallocate elements of the sequence {1,2,3...,N}  
parsum=linspace(0,0,N);% preallocate the partial sums: {1, 1+2, 1+2+3, ...}  
seq(1)=1*z;  % first element of sequence  
parsum(1)=seq(1); % first partial sum  
%All the partial sums are stored in parsum, all the terms in seq
for count=2:N  
    seq(count)=count*z^count;  
    parsum(count)=parsum(count-1)+seq(count);  
end  
parsumR=real(parsum); % we plot only the real part  
subplot(3,1,1)
plot(n,parsumR,'.k',[0, N],[-1/12,-1/12],'MarkerSize',6.0,'Linewidth',1.5); 
axis([0,shortN,0,shortN*(shortN+1)/2]);
%title('(a) The initial partial sums proceed almost as 1+2+3+4+...');
% The titles do not print properly on Octave and had to be commented out
%plots the first terms (should go as n*(n+1)/2) Verify with parabola:
hold on;
plot(n,n.*(n+1)/2,'-r',[0, N],[-1/12,-1/12]);
hold off;
subplot(3,1,2)
plot(n,parsumR,'.k',[0, N],[-1/12,-1/12],'MarkerSize',8.0); 
%title('(b) The small deviation from z=1 introduces decay and oscillation'); %plots all N terms
subplot(3,1,3)
plot(n,parsumR,'.k',[0, N],[-1/12,-1/12],'MarkerSize',8.0); 
axis([0,N,-.25,.25]);
%title('(c) The last terms show convergence to approximately -1/12'); %plots the partial sums for the largest values of n.
final=parsum(N);
fprintf('If z= \n%11.9f+%11.9fi ',real(z),imag(z) );
fprintf('the sum from 1 to %d of n*z^n is \n',N ); 
fprintf(' %11.9e + i %11.9e\n',real(final),imag(final));%this works
%%%%Save the real part of the sequence
for count=1:N  
    Rseq(count)=real(seq(count));    
end