######### EQ. DIFFERENZIALI del 1º ORDINE ######### # f'(x) = -x/f(x), ovvero y' = -x/y # Traccio tanti trattini con pendenza -x/y (in questo caso le soluzioni della # equazione hanno tutte come grafico dei mezzi cerchi con centro nell'origine). # Ottengo il "campo direzionale" associato alla eq. differenziale. # Dy <- function(x,y) -x/y # riduco i margini della pagina par( mai = c(0.5,0.5,0.2,0.2) ) # # tracciamento in (x1,x2)*(y1,y2) del campo direzionale di y' = f(x,y) # con f messo in Dy; prova con m=25, n=25 # (m e n indicano quanto sono i trattini in orizz. e in vert.) soledif <- function(x1,x2,y1,y2,m,n) { plot(c(x1,x2), c(y1,y2), type="n",xlab="",ylab="") abline(h=0,v=0,col="blue") dx <- (x2-x1)/m; dy <- (y2-y1)/n a <- 0; b <- 0; f <- function(x) Dy(a,b)*(x-a)+b for(i in 0:m) for(j in 0:n) {a <- x1+dx*i; b <- y1+dy*j Dx <- dx; for(k in 1:100) {if(abs(f(a+Dx/2.5)-f(a-Dx/2.5)) > dy) Dx <- Dx*0.9 else k <- 100} plot(f,a-Dx/2.5,a+Dx/2.5,ylim=c(b-dy/2.5,b+dy/2.5),add=TRUE,col="grey55")}} # soledif(-20,20,-20,20, 25,25) # # il grafico di alcune soluzioni h <- function(x) sqrt(a^2-x^2); k <- function(x) -h(x) for(a in c(10,20)) for(t in c(h,k)) plot(t,-a,a,add=TRUE,col="red") # # Altro esempio: y'(x) = x-y(x) Dy <- function(x,y) x-y x1 <- -3; x2 <- 4 soledif(x1,x2,-3,3, 25,25) # # Alcune soluzioni (risolvendo, in questo caso, l'eq. differenziale; # e' facile verificare che le seguenti sono soluzioni) f <- function(x) x-1+(C+1)*exp(-x) C <- -2; points(0,C,pch=20); plot(f,x1,x2,col="red",add=TRUE) C <- -1; points(0,C,pch=20); plot(f,x1,x2,col="red",add=TRUE) C <- 0; points(0,C,pch=20); plot(f,x1,x2,col="red",add=TRUE) C <- 1; points(0,C,pch=20); plot(f,x1,x2,col="red",add=TRUE) # # Se soledif(a,b,c,d, m,n) si arresta per l'incontro di punti in cui c'e' # una divisione per zero (o altri "fuori dominio") si puo' tenere # conto di questo e, poi, riavviare soledif(a,b*1.01,c,d*1.01, m,n) # # (E' bene controllare le soluzioni proposte da libri di esercizi di # analisi con questo programmino: a volte esse sono errate o parziali) # # Altro esempio # # Circuito RL L*D(ic)(t)+R*ic-V(t)=0 (L induttanza in Henry, # R resist. in Hom, ic intens. corrente in Ampere, V tensione in Volt) Dy <- function(x,y) -R*y/L+V(x)/L L <- 1; R <- 2; V <- function(t) 6+t-t # metto 6+t-t per avere 6 par( mai = c(0.5,0.5,0.2,0.2) ) soledif(0,10,-1,4, 25,25) # dev.new(); par( mai = c(0.5,0.5,0.2,0.2) ) V <- function(t) 6*cos(2*t) soledif(0,15,-3,3,2*pi*15,30) # # Controllo della soluzione per (0,0) [trovata con # tecniche di analisi] k <- function(t) 3/2*cos(2*t)+sin(2*t)-exp(-2*t)*3/2 plot(k,0,15,add=TRUE) axis(1,pos=0,col="red",label=FALSE,at=seq(0,15,pi/2)) abline(h=3/2*cos(pi/4)+sin(pi/4),lty=3) abline(h=-(3/2*cos(pi/4)+sin(pi/4)),lty=3) # # Come risolvere graficamente un'equazione differenziale del # 1° ord. col metodo di Eulero: # Il tracciamento del poligono approssimante il graf. della sol. # (quella passante per (x0,y0), da x0 a xf). soledifc <- function(x0,y0,xf,n,col) { x <- NULL; y <- NULL; h <- (xf-x0)/n; y[1] <- y0 for (i in 1:(n+1)) x[i] <- x0+h*(i-1) for (i in 2:(n+1)) y[i] <- y[i-1]+Dy(x[i-1],y[i-1])*h lines(x[1:(n+1)],y[1:(n+1)],lwd=2,col=col); points(x0,y0,pch=21,bg="black") } # Il tracciamento dei suoi vertici soledifp <- function(x0,y0,xf,n,col) { x <- NULL; y <- NULL; h <- (xf-x0)/n; y[1] <- y0 for (i in 1:(n+1)) x[i] <- x0+h*(i-1) for (i in 2:(n+1)) y[i] <- y[i-1]+Dy(x[i-1],y[i-1])*h for(i in 1:(n+1)) points(x[i],y[i],pch=20,col=col)} # Esempio: l'equazione y' = x-y Dy <- function(x,y) x-y soledif(-3,4,-3,3, 25,25) # Usando la stessa scala, il grafico di alcune solzioni: soledifp(-3,2,4,10,"red") soledifc(-3,2,4,10,"blue") soledifc(-1/2,1,4,1e4,"brown") soledifc(-1/2,1,-3,1e4,"brown") soledifc(4,2,-3,1e4,"green") soledifc(4,2.5,-3,1e4,"red") soledifc(0,-2,4,1e4,"orange") # Altro esempio (che evidenzia come al crescere dei passi si arriva # ad una soluzione grafica accettabile - esistono metodi in cui la # convergenza al grafico OK è più veloce) Dy <- function(x,y) sin(x)*(1-y) plot(c(0,50), c(0,2), type="n",xlab="",ylab="") abline(h=0,v=0,col="blue") soledifc(0,2,50,100,"red") soledifc(0,2,50,500,"blue") soledifc(0,2,50,2000,"orange") soledifc(0,2,50,50000,"black") # Posso risolvere, per controllo, l'eq, con WoframAlpha: # y'(x) = sin(x)*(1-y(x)), y(0) = 2 # Ottengo: 1+e^(-1+cos(x)). Poi eseguo plot y = 1+e^(-1+cos(x)), x = 0..50, y =0..2 # e ho: