%%  SymsumEsim.m
% Symbolista summausta, Solmu 2/2019, Heikki Apiola
%
clear;close all
format rat  % Pelkkiä rationaalisia laskuja
ind=(1:10)'; % Pystyvektori
jono=(ind.^5);
S=cumsum(jono); % Kumulatiiviset summat
indeksit_jono_ja_osasummat=[ind jono S]
%
%% Osasummien erotukset
format rat  % Pelkkiä rationaalisia laskutoimituksia
ydiff1=[NaN;diff(S)]; % Sama kuin jono
ydiff2=[NaN;1/2*diff(ydiff1)]; % 2. jaetut erotukset
ydiff3=[NaN;1/3*diff(ydiff2)]; % 3. 
ydiff4=[NaN;1/4*diff(ydiff3)]; % ...
ydiff5=[NaN;1/5*diff(ydiff4)]; % ...
ydiff6=[NaN;1/6*diff(ydiff5)]; % ...
erotustaulukko=[S ydiff1 ydiff2 ydiff3 ydiff4 ydiff5 ydiff6];
c=diag(erotustaulukko);c' % Newtonin polynomin kertoimet
%%
syms n
p_n=c(1)+c(2)*(n-1)+c(3)*(n-1)*(n-2)+c(4)*(n-1)*(n-2)*(n-3)+...
    c(5)*(n-1)*(n-2)*(n-3)*(n-4)+c(6)*(n-1)*(n-2)*(n-3)*(n-4)*(n-5)+...
    c(7)*(n-1)*(n-2)*(n-3)*(n-4)*(n-5)*(n-6);
%% Induktioaskel:
p_nplus1=subs(p,n,n+1) % Sijoitetaan p_n:n lausekeessa
            % n:n paikalle n+1
% Päteekö: p(n+1) = p(n)+(n+1)^5 % Luettavammmin merkinnöin
erotus=p_nplus1-p_n
simplify(erotus)
display('MOT')
%% Sama Lagrangella
% Menee lyhyemmin, koska tehtiin valmis funktio: 
% LagrangeInterp
syms n
p=LagrangeInterp(ind,S,n)
p=simplify(p)
latex(p)  % kokeillaan, hieno väline 
% (myös Maplessa, Mathematicassa
%% Induktioaskel:
p_n=p
p_nplus1=subs(p,n,n+1)
% Päteekö: p(n+1) = p(n)+(n+1)^5
% Ts.: Onko p(n+1)-p(n) = (n+1)^k ?
erotus=p_nplus1-p_n
simplify(erotus)