---------- ---------- ---------- ---------- ---------- ---------- ---------- ----------
#
# il simbolo Inf rappresenta l'infinito. Viene assunto il "valore" Inf quando un numero
# eccede un certo valore.
factorial(170); factorial(171)
# 7.257416e+306  Inf
# Analogamente viene assunto il valore 0 quando il valore assoluto di un numero è molto
# piccolo. Per dettagli aziona:
.Machine
#
h <- function(x) exp(-x)
Piano(0,10, -0,1); grafic(h, 0,10, "brown")
integrale(h,0, Inf)
# 1
  
# il limite per x -> Inf di h è 0:
h(1); h(10); h(100); h(Inf)
# 0.3678794  4.539993e-05  3.720076e-44  0
#
# Attenzione: il limite di f(x) per x -> Inf è stimabile con f(Inf) non in
# tutti i casi. È beve valutare f(x) per x che cresce per controllare. Es.:
f <- function(x)  (1+1/x)^x; f(Inf)
# 1
# mentre il limite è e:
n <- 1:8; f(10^n)
# 2.593742 2.704814 2.716924 2.718146 2.718268 2.718280 2.718282 2.718282
#
# I limiti di q(x) per q che tende ad 1 da destra e da sinistra possono essere diversi
q <- function(x) (x^2-1)/(x^2-2*x+1); Piano(-2,2, -15,15); grafico(q,-2,2, "red")
 
q(1-1E-2);q(1-1E-5);q(1-1E-10)
# -199  -199999  -Inf
q(1+1E-2);q(1+1E-5);q(1+1E-10)
# 201  200001  Inf
#
# Vi sono successioni [cioè funzioni F che hanno come input valori interi crescenti a
# partire da un certo valore, ad es. 0,1,2,3,… o -2,-1,0,1,… o 6,7,8,9,…] tali che F(n)
# non può essere espresso mediante un semplice termine. Consideriamo ad esempio
#   n((n+1)·(n+2)·...·(n+n))/n   per n = 1, 2, 3, …
# Vediamo come studiare il limite di F(n) per n che tende all'infinito. Mi conviene
# innanzi tutto scrivere il termine nel modo più semplice per effettuare i calcoli. In
# questo caso particolare ho un prodotto di termini man mano crescenti sia in valore
# che in quantità. Mi conviene riscrivere F(n) come n√(n+1)·n√(n+2)·...·n√(n+n)/n
# Traduco tutto ciò in R, usando un ciclo for per esprimere il prodotto di più termini:
F = function(n) {x=1; for(i in 1:n) x = x*(n+i)^(1/n); x/n}
# Studio quello che accade. Provo con n potenza di 2. Prima provo anche con n+1 per
# controllare che non ci siano differenze per n dispari (con le successioni, specie se
# ci sono delle potenze, questo può accadere, anche se potrei capire che non è questo
# il caso).
n=2; c(F(n), F(n+1) )
# 1.732051 1.644141
n=2^2; c(F(n), F(n+1) )
# [1] 1.600543 1.574513
n=2^3; c(F(n), F(n+1) )
# [1] 1.535668 1.528503   tutto OK, elimino F(n+1)
n=2^18; F(n)
# [1] 1.47152
n=2^19; F(n)
# [1] 1.471519
n=2^20; F(n)
# [1] 1.471518
n=2^21; F(n)
# [1] 1.471518
# Vedo con WolframAlpha se posso esprimere questo termine in modo più sintetico.
# Ottengo 4/e. Verifica:
4/exp(1)
# [1] 1.471518   OK!
#
# Lo studio di una serie:  1-1/2+1/3-1/4+1/5...
# Il grafico di N -> 1-1/2+ ... 1/N
a <- function(n) (-1)^(n+1)/n  # a(n) elemento n-esimo della sommatoria
N <- function(n) seq(1,n,1)    # N = 1 2 ... n
S <- function(n) sum(a(N(n)))  # somma a(1)+...a(n)
NMax <- 20; Piano(0,NMax, 0,1)
for(i in 1:NMax) PUNTI(i,S(i), "brown")
for(i in 1:(NMax-1)) spezza(c(i,i+1),c(S(i),S(i+1)), "brown")
  
# A quale valore tende questa serie?
# Stampo i valori della sommatoria e del termine successivo che verrebbe sommato
n <- 10; piu(S(n)); a(n+1)
# 0.645634920634921  0.09090909
n <- 100; piu(S(n)); a(n+1)
# 0.688172179310195  0.00990099
n <- 1e6; piu(S(n)); a(n+1)
# 0.693146680560195  9.99999e-07
> n <- 1e7; piu(S(n)); a(n+1)
# 0.693147130559948  9.999999e-08
> n <- 1e8; piu(S(n));a (n+1)
# 0.693147175559945  1e-08
> piu(log(2))
# 0.693147180559945
#
# È facilmente realizzabile lo studio di serie di Fourier (vedi qui e qui)
        
# (esitono anche specifici comandi di R per studiarle). Un esempio:
#      4/π*sin(x) + 4/π*sin(3x)/3 + 4/π*sin(5x)/5 + ...
p <- function(x) {f <- 0;
  for(i in 0:n) f <- f+sin((2*i+1)*x)/(2*i+1); f <- f*4/pi}
BF=6.5; HF=2.3
n=0; graficF(p, -5,10, "black")
 
n=1; graficF(p, -5,10, "blue")
n=2; grafic(p, -5,10, "red")
n=3; grafic(p, -5,10, "seagreen")
 
BF=6.5; HF=2.3
n=20; graficF(p, -5,10, "brown")
 
BF=6.5; HF=2.3
n=2000; graficoF(p, -5,10, "brown")
 
#
# Per ottenere da una funzione F la sua serie (o i polinomi) di Fourier puoi
# semplicemente riccorrere a WolframAlpha (vai qui e batti "fourier series")
# Nel caso in questione se metto in  function to expand
#       Piecewise[ { {-1, -pi < x < 0}, {1, 0 < x < pi} } ]
# (e lo scalino rappresentato sopra, che poi si ripete periodicamente)
# e in  oder  1, 3 e 21, ottengo i grafici seguenti e i relativi polinomi

# Ad esempio per order=3 ottengo il polinomio considerato sopra:
        
# Se rappresentiamo l'ampiezza delle sinusoidi, invece che in funzione di "x"
# come fatto inizialmente, in funzione della frequenza otteniamo dei segmenti
# verticali. Ecco che cosa otteniamo il polinomio di Fourier costituito da
# quattro addendi:
 
Piano(0,0.2, -1.3,1.3)
segm(1/(2*pi),0, 1/(2*pi),4/pi, "black")
segm(1/(2*pi*3),0, 1/(2*pi*3),4/(pi*3), "blue")
segm(1/(2*pi*5),0, 1/(2*pi*5),4/(pi*5), "red")
segm(1/(2*pi*7),0, 1/(2*pi*7),4/(pi*7), "seagreen")
# Questo è il cosiddetto spettro del polinomio.
# Il concetto di "spettro" viene usato in tutte le situazioni in cui una funzione
# può essere rappresentata come una somma di funzioni sinusoidali:   la funzione
# "rossa" somma delle tre sinusoidi rappresentate, e lo spettro corrispondente:
 

BF=3.5; HF=2.3
g=function(t) sin(1200*pi*t)
h=function(t) 2*sin(800*pi*t)
k=function(t) 3*sin(400*pi*t)
a=-0.005; b=0.005; Piano(a,b, -5,5)
grafi(g, a,b, "black"); grafi(h, a,b, "blue"); grafi(k, a,b, "green3")
f=function(t) g(t)+h(t)+k(t);  grafi(f, a,b, "red")
#
Piano(0,800, -5,5)
Segm(600,0, 600,1, "black"); Segm(400,0, 400,2, "blue")
Segm(200,0, 200,3, "green3")
# le freq. sono 1/periodo, ossia 1200π/(2π) = 600, …