##########           EQUAZIONI DIFFERENZIALI del 2º ORDINE           ##########
# Per equazioni differenziali di ordine superiore conviene usare WolframAlpha
# per trovare le soluzioni, e poi R per tracciarne il grafico.
#
# Considero prima un esempio che sono in grado di affrontare facilmente. Posso
# pensare che corrisponda alla situazione di un'auto che si sta muovendo (all'
# istante che assumo come iniziale, in cui dista 2 m dal punto che assumo come
# riferimento ed ha velocità 1 m/s) con un'accelerazione costante di 5 m/s^2.
# Se f è la posizione all'istante x: f"(x)=5, f'(0)=1, f(0)=2.  Posso ricavare
# che  f'(x)=5*x+h  e che 1=5*0+h, da cui h=1.  Da  f'(x)=5*x+1  ricavo  f(x)=
# 2.5*x^2+x+k;  da  f(0)=2  ho  k=2.  Quindi  f(x)=2.5*x^2+x+2.
# Che cosa otterrei con WolframAlpha?  Introduco  f"(x)=5, f'(0)=1, f(0)=2  e
# ottengo   f(x) = 5*x^2/2+x+2.    Eccone il grafico con R:
f <- function(x) 5*x^2/2+x+2; plot(f,0,20,n=3000,lwd=2)
abline(v=axTicks(1), h=axTicks(2), col="blue",lty=3)

#
# Facciamo alcuni altri esempi riferiti sempre alla fisica (vedi).
#
# Secondo esempio. x = tempo, y = posizione di oggetto di massa m sottoposto a
# una forza di -k*y (ad es. esercitata da una molla) con k = 4, m = 1; da F =
# -k*y(x), F = m*y"(t) ho  y"(x) = -4/1*y(x).  Scrivo quanto sopra in
# WolframAlpha. Ottengo:  y(x) = c2*sin(2*x)+c1*cos(2*x).  Se inserisco una
# posizione e una velocità iniziali y"(x) = -4*y(x), y(0)=2, y'(0)=0
# Ottengo  y(x) = 2*cos(2*x)  e i grafici seguenti, che rivedremo con R:

# Eccone il grafico della funzione soluzione con R:
plot(f,0,20,n=3000,lwd=2)
abline(v=axTicks(1), h=axTicks(2), col="blue",lty=3)

# Vediamo come tracciare il diagramma di fase, ossia la curva (y(x),y'(x))
# che rappresenta la relazione tra "y" e la sua "velocità".
g <- function(x) eval(D(body(f),"x"))
# individuo i range in cui variano y(x) e y'(x)
x1 <- 0; x2 <- 20; punti <- 3000; x <- seq(x1,x2,(x2-x1)/punti)
c(min(f(x)),max(f(x)),min(g(x)),max(g(x)))
[1] -1.999998  2.000000 -4.000000  3.999997
dev.new()   # apro una nuova finestra per mantenere il graf. precedente
plot(c(-2,2), c(-4,4), type="n", xlab="", ylab="") 
abline(v=axTicks(1), h=axTicks(2), col="blue",lty=3)
lines(f(x),g(x),lwd=2)

# In questo caso come diagramma di fase ottengo un'ellisse. Il fatto che sia
# una curva chiusa è legato alla periodicità della soluzione
#
# Terzo esempio. Oscillazioni forzate, con forza oscillante esterna, del tipo
# F·cos(Ω·t):  y" = -(k/m)·y + F/m·cos(Ω·t). Consideriamo il caso: F=1, Ω=2
# e gli altri dati come sopra:  y"(x)=-4*y(x)+cos(2*x), y(0)=2, y'(0)=0
# Scrivo quanto sopra in WolframAlpha. Ottengo:    2*cos(2*x)+1/4*x*sin(2*x)
# e i grafici seguenti, che rivedremo con R:

# Eccone il grafico della funzione soluzione con R:
f <- function(x) 2*cos(2*x)+1/4*x*sin(2*x)
plot(f,0,30,n=3000,lwd=2)
abline(v=axTicks(1), h=axTicks(2), col="blue",lty=3)

# Il diagramma di fase, ossia la relazione tra "y" e la sua "velocità".
g <- function(x) eval(D(body(f),"x"))
# individuo i range in cui variano y(x) e y'(x)
x1 <- 0; x2 <- 30; punti <- 3000; x <- seq(x1,x2,(x2-x1)/punti)
c(min(f(x)),max(f(x)),min(g(x)),max(g(x)))
[1]  -7.123881  7.503536    -15.328534  14.564993
dev.new()   # apro una nuova finestra per mantenere il graf. precedente
plot(c(-8,8), c(-16,16), type="n", xlab="", ylab="") 
abline(v=axTicks(1), h=axTicks(2), col="blue",lty=3)
lines(f(x),g(x),lwd=2)

# La parte rossa corrisponde alle prime "x". E' stata tracciata,
# per evidenziare come è stata generata la curva, con:
x1 <- 0; x2 <- 15; punti <- 3000; x <- seq(x1,x2,(x2-x1)/punti)
lines(f(x),g(x), lwd=2,col="red")
# Per quasti valori di F ed Ω si ha il fenomeno della risonanza (le
# oscillazioni aumentano di ampiezza).
#
# Quarto esempio. Oscillazioni smorzate con  y" = -(k/m)·y - k1·y' (con forza
# di attrito proporzionale alla velocità). Esempio con k1 = 1
# y"(x) = -4*y(x)-y'(x), y(0)=1, y'(0)=3
# Con WolframAlpha ottengo:
# y(x) = 1/15*exp(-x/2)*(7*sqrt(15)*sin((sqrt(15)*x)/2)+15*cos((sqrt(15)*x)/2))

# In R posso fare:
S <- sqrt(15)
f <- function(x) 1/15*exp(-x/2)*(7*S*sin((S*x)/2)+15*cos((S*x)/2))
plot(f,0,20,n=3000)
abline(v=axTicks(1), h=axTicks(2), col="blue",lty=3)
g <- function(x) eval(D(body(f),"x"))
dev.new()
c(min(f(x)),max(f(x)),min(g(x)),max(g(x)))
[1] -0.7204572  1.6214072 -2.3073632  3.0000000
plot(c(-0.8,1.7), c(-2.4,3), type="n", xlab="", ylab="")
# aggiungo asp=1 se voglio il sistema monometrico
abline(v=axTicks(1), h=axTicks(2), col="blue",lty=3)
x1 <- 0; x2 <- 20; punti <- 3000; x <- seq(x1,x2,(x2-x1)/punti)
lines(f(x),g(x),lwd=2)