########## 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)