# Tracciamento e calcolo dei polinomi di Taylor di G attorno ad un punto A.
# Definisci una funzione G, scegli A, un intervallo [B,C] e batti
#  taylor(G, A, B,C, H,K)  per ottenere i grafici di G e, via via che CLICCHI
# lo schermo, dei pol. di T. fino a quello di grado 10  in  [B,C] x [H,K]; prova
# con una scelta di H e K e se è il caso con una scelta diversa.
# Vengono poi stampati i coeff. del pol., dal grado 0 al grado 10.
# Se vuoi vedere G e le sue derivate 1^,2^,... (non semplificate) batti
#   Fx   poi   F1x   poi   F2x   poi   F3x   ...
# Se vuoi vedere l'eventuale forma frazionaria del coeff. ad es. di grado 3 e 4 se
# A=1 batti    A <- 1; fractions( c( F3(A)/factorial(3),F4(A)/factorial(4) )
# [fractions(Fn(A)/factorial(n)) è approssimato; può differire dal valore corretto]
# Se vuoi fissare l'immagine fino al polinomio ad es. di grado 4 vai sulla
# finestra di testo e premi ESC.
#
# Copia e incolla le righe da ####### a #######
# Quindi vedi gli esempi successivi.
#######
# Sottoprogramma che traccia i polinomi di Taylor fino a quello di ordine 10
taylor <- function(f,A,x1,x2,y1,y2) {
library(MASS)      # Il grafico di f:
plot(f, x1,x2, ylim=c(y1,y2),lwd=3,n=1000)
abline(v=axTicks(1),h=axTicks(2),lty=3, col="grey50")
abline(v=0, h=0, lty=2)
# le derivate di f e i polinomi di Taylor
c <- c("blue","grey50","brown","red","orange","green3","red","violet","cyan","yellow4")
j <- 0
Fx <<- body(f); F0 <<- function(x) f(x)
p <- locator(1)
F1x <<- D(Fx,"x"); F1 <<- function(x) eval(F1x)
 p1 <- function(x) f(A)+F1(A)*(x-A)
 j <- j+1; curve(p1,add=TRUE,col=c[j],lty=2,lwd=2)
p <- locator(1)
F2x <<- D(F1x,"x"); F2 <<- function(x) eval(F2x)
 p2 <- function(x) p1(x)+F2(A)*(x-A)^2/factorial(2)
 j <- j+1; curve(p2,add=TRUE,col=c[j],lty=2,lwd=2)
p <- locator(1)
F3x <<- D(F2x,"x"); F3 <<- function(x) eval(F3x)
 p3 <- function(x) p2(x)+F3(A)*(x-A)^3/factorial(3)
 j <- j+1; curve(p3,add=TRUE,col=c[j],lty=2,lwd=2)
p <- locator(1)
F4x <<- D(F3x,"x"); F4 <<- function(x) eval(F4x)
 p4 <- function(x) p3(x)+F4(A)*(x-A)^4/factorial(4)
 j <- j+1; curve(p4,add=TRUE,col=c[j],lty=2,lwd=2)
p <- locator(1)
F5x <<- D(F4x,"x"); F5 <<- function(x) eval(F5x)
 p5 <- function(x) p4(x)+F5(A)*(x-A)^5/factorial(5)
 j <- j+1; curve(p5,add=TRUE,col=c[j],lty=2,lwd=2)
p <- locator(1)
F6x <<- D(F5x,"x"); F6 <<- function(x) eval(F6x)
 p6 <- function(x) p5(x)+F6(A)*(x-A)^6/factorial(6)
 j <- j+1; curve(p6,add=TRUE,col=c[j],lty=2,lwd=2)
p <- locator(1)
F7x <<- D(F6x,"x"); F7 <<- function(x) eval(F7x)
 p7 <- function(x) p6(x)+F7(A)*(x-A)^7/factorial(7)
 j <- j+1; curve(p7,add=TRUE,col=c[j],lty=2,lwd=2)
p <- locator(1)
F8x <<- D(F7x,"x"); F8 <<- function(x) eval(F8x)
 p8 <- function(x) p7(x)+F8(A)*(x-A)^8/factorial(8)
 j <- j+1; curve(p8,add=TRUE,col=c[j],lty=2,lwd=2)
p <- locator(1)
F9x <<- D(F8x,"x"); F9 <<- function(x) eval(F9x)
 p9 <- function(x) p8(x)+F9(A)*(x-A)^9/factorial(9)
 j <- j+1; curve(p9,add=TRUE,col=c[j],lty=2,lwd=2)
p <- locator(1)
F10x <<- D(F9x,"x"); F10 <<- function(x) eval(F10x)
 p10 <- function(x) p9(x)+F10(A)*(x-A)^10/factorial(10)
 j <- j+1; curve(p10,add=TRUE,col=c[j],lty=2,lwd=2)
# il punto ( A, f(A) )
points(A,f(A),pch=19,col="red")
# i coefficienti del polinomio f(A) + F1(A)*(x-1) + F2(A)*(x-1)^2 + F3(A)*(x-1)^3 + ...
cat(f(A),' ',F1(A),' ',F2(A),'/',factorial(2),' ',F3(A),'/',factorial(3),' ',F4(A),'/',
    factorial(4),' ',F5(A),'/',factorial(5),' ',F6(A),'/',factorial(6),' ',F7(A),'/',
    factorial(7),' ',F8(A),'/',factorial(8),' ',F9(A),'/',factorial(9),' ',F10(A),'/',
    factorial(10),'\n')
}
##  taylor(G, A, B,C, H,K)
##  fractions( F3(A)/factorial(3) )
#######
#
## sin in 0
f <- function(x) sin(x); dev.new(); taylor(f, 0, -5,5, -5,5)
# deduco che attorno a 0 sin(x) = x - 1/6*x^3 + 1/120*x^5 - 1/5040*x^7 + 1/362880*x^9 + …
Fx; F1x; F2x; F3x; F4x
## log in 1
g <- function(x) log(x); dev.new(); taylor(g, 1, 0,4, -4,3)
# deduco che attorno a 1 log(x) = x-1 - (x-1)^2*1/2 + (x-1)^3*2/6 - (x-1)^4*6/24 + (x-1)^5*24/120 + …
A <- 1; fractions( c(F3(A)/factorial(3),F4(A)/factorial(4),F5(A)/factorial(5),F6(A)/factorial(6) ) )
#  1/3 -1/4  1/5 -1/6
# 
## cos in 1
h <- function(x) cos(x); dev.new(); taylor(h, 1, -5,5, -2,2)
## exp in 0
k <- function(x) exp(x); dev.new(); taylor(k, 0, -5,5, -2,15)
Fx; F1x