#
# Risolvere graficamente e numericamente EQ. DIFFERENZIALI del 1º ORDINE
#
## ======== Copia da qui a "FINE" =========
# CAMPO DIREZIONALE
# Tracciamento in (x1,x2)*(y1,y2) del c. direz. di y' = f(x,y) con f messo in Dy:
# Dy <- function(x,y) ... ; prova con m=25, n=25 (m+1 e n+1 sono il numero dei
# trattini in orizz. e in vert.; se metti m=0, n=0 puoi tracciare solo il
# reticolo e limitarti ad usare soledifc per tracciare la curva soluz.)
col1 <<- "orange"; col2 <<- "grey52" # colori cambiabili di griglia e campo
soledif <- function(x1,x2,y1,y2,m,n) {
plot(c(x1,x2), c(y1,y2), type="n",xlab="",ylab="")
abline(v=axTicks(1),h=axTicks(2),lty=3,col=col1)
abline(h=0,v=0,col="blue"); if(n*m > 0) {
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=col2)}}}
#
# CURVA SOLUZIONE
# Soluzione di y' = f(x,y) tale che y(x0) = y0 fino a xf col metodo di
# Eulero in N passi: sia h = (xf-x0)/N, approssimo y(x) col segmento che parte
# da (x0,y0), ha pendenza f(x0,y0), arriva a x1=x0+h; sia y1 la corrispondente
# ordinata. Poi col segmento da (x1,y1) con pendenza f(x1,y1) fino a x2=x1+h; …
### Il tracciamento del poligono approssimante il graf. della sol. y(x)
### 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)}
### Il valore di y(x) in p approssimato in n passi
soledifv <- function(x0,y0, p, n) {xxx <<- NULL; xx <- 0; D2 <- 0
x<- NULL; y<- NULL; h<- (p-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; D1 <- D2; D2 <- y[i]-y[i-1];
if(xx==0) {if(sign(D1)+sign(D2)==0){xx<- 1; xxx<<- c(x[i-1],x[i],y[i-1],y[i]) } } }
# in xxx sono messi i punti in cui cambia segno Δy; vedi Esempio 3
y[n+1] } ### In breve i comandi:
# Dy<- function(x,y) … soledif(,, ,, ,) soledifc(,, , ,) soledifp(,, , ,) soledifv(,, ,)
## =============== "FINE" ================
#
# ESEMPIO 1
# y'(x) = x - y(x) Trovare il campo direzionale in (-1,5)*(-2,5)
# e tracciare le curve soluzioni tali che y(1/2)=3 e che y(2)=1/2.
#
# Definisco l'eq. differenziale [nota: metto solo x-y, non x-y(x)]:
#
Dy <- function(x,y) x-y
#
# Scelgo le dimensioni della finestra e il margine (cosa non necessaria)
dev.new(width=4,height=3); par( mai = c(0.5,0.5,0.2,0.2) )
#
col1 <<- "orange"; col2 <<- "grey52"
soledif(-1,5, -2,4, 30,25)
# Traccio la soluz. approssimata fino a x=5 in 3 passi
soledifc(1/2,3, 5, 3, "red")
# Traccio i pallini estremi dei segmenti
soledifp(1/2,3, 5, 3, "blue")
# Aumento il numero dei passi (7)
soledifc(1/2,3, 5, 7, "green2")
soledifp(1/2,3, 5, 7, "blue")
# Con un numero alto di passi ho il grafico OK
soledifc(1/2,3, 5, 1e4, "green4")
# Traccio la soluzione anche fino a x=-1
soledifc(1/2,3, -1, 1e4, "brown")
#
# Per avere solo i grafici delle soluzioni definisco la scala e uso
# soledifc. Ma posso sceglere la scala anche usando soledif con m=n=0
soledif(-1,5, -2,4, 0,0)
soledifc(1/2,3, 5, 1e4, "green4")
soledifc(1/2,3, -1, 1e4, "brown")
soledifc(2,1/2, 5, 1e4, "green4")
soledifc(2,1/2, -1, 1e4, "brown")
#
# Per avere il valore della soluzione in un qualunque punto uso solediv
# Vediamo ad esempio come trovare il valore che la soluzione passante
# per (0,1) assume nel punto 2
soledif(-1,5, -2,4, 0,0)
soledifc(0,1, 5, 1e4, "green4")
soledifc(0,1, -1, 1e4, "brown")
# Faccio un po' di prove e assumo 1.2707 come valore
soledifv(0,1, 2,1e4)
# 1.270616
soledifv(0,1, 2,2e4)
# 1.270643
soledifv(0,1, 2,4e4)
# 1.270657
points(2,1.2707, col="red",pch=20)
#
## Se so risolvere l'equazione - la soluzione in questo caso è:
f <- function(x) x-1+(C+1)*exp(-x); C <- 1
## - posso verificare che il valore è corretto:
f(2)
# 1.270671
# E sovrapporre il grafico:
curve(f,col="yellow",add=TRUE,lty=3,lwd=2)
#
# ESEMPIO 2
# Vediamo come trovare con più precisione il valore che la soluzione della
# eq. diff.le precedente passante per (0,1) assume nel punto 2
v0 <- soledifv(0,1, 2, 5e3); v0
# 1.270562
v1 <- soledifv(0,1, 2, 1e4); c(v1, v1-v0)
# 1.270616e+00 5.413953e-05
v2 <- soledifv(0,1, 2, 2e4); c(v2, v2-v1)
# 1.270643e+00 2.706841e-05
v3 <- soledifv(0,1, 2, 4e4); c(v3, v3-v2)
# 1.270657e+00 1.353387e-05
## In questo caso le variazioni man mano si dimezzano: 1/2+1/4+1/16+... = 1
## Quindi l'ultima variazione fornisce una stima di quanto manca al valore esatto
v3 + (v3-v2)
# 1.270671
## Ritrovo il valore trovato risolvendo l'equazione differenziale
#
#
# ESEMPIO 3
# y'(x) = x/y Trovare il campo direzionale in (-3,3)*(-3,3)
# e tracciare le curve soluzioni tali che y(-2)=3 e che y(-3)=2.
#
Dy <- function(x,y) x/y
col1 <<- "red"; col2 <<- "grey50"
soledif(-3,3, -3,3, 25,25)
soledifc(-2,3, 3, 1e4, "brown")
v0 <- soledifv(-2,3, 3, 5e3); v0
# 3.741161
v0 <- soledifv(-2,3, 2, 5e3); v0
# 2.999565
v1 <- soledifv(-2,3, 2, 1e4); c(v1, v1-v0)
# 2.9997824339 0.0002175728
v2 <- soledifv(-2,3, 2, 2e4); c(v2, v2-v1)
# 2.9998912178 0.0001087839
v3 <- soledifv(-2,3, 2, 4e4); c(v3, v3-v2)
# 2.999946e+00 5.439132e-05
# Anche in questo caso le variazioni si dimezzano; quindi assumo
v3 + (v3-v2)
# 3 Come era ovvio, data la simmetria della curva soluzione
#
# Vediamo l'altro problema:
soledifc(-3,2, 3, 1e4, "green4")
# Era evidente dal campo direzionale e dalla eq. stessa che la curva non è
# definita per x≥0. Studiamo meglio il dominio (potremmo farlo diversamente,
# date le simmetrie del problema, ma facciamolo prescindendo dalle particolarità
# di esso).
# Sfrutto il fatto che, chiamato soledifv (vedi il commento presente in soledifv),
# xxx mi dà i punti in cui la pendenza della curva soluzione cambia segno;
# Capisco dal grafico verde ottenuto che sarà circa 2.25.
soledifv(-3,2, 3, 1e4); xxx
# -2.234400000 -2.233800000 -0.034917451 0.003477104
# Tra (-2.2344,-0.034917451) e (-2.2338,0.003477104) cambia segno la pendenza; preciso
# l'intervallo
soledifv(-3,2, -2.2, 2e4); xxx
# -2.23592000 -2.23588000 -0.02072017 -0.01640376. x è circa -2.2359
# Evidenzio il punto
xp <- mean(xxx[1],xxx[2]); yp <- mean(xxx[3],xxx[4])
points(xp,yp, pch=19,col="red")
# Traccio il grafico della soluzione e il grafico di quella passante per (-3,-2):
soledifc(-3,2, xp, 1e4, "black")
soledifc(-3,-2, xp, 1e4, "black")
points(xp,yp, pch=19,col="red")
#
# Riferendoci alla stessa equazione illustriamo un altro problema:
# se vi sono dei punti in cui l'equazione non è definita può capitare che
# essi non vengano individuati dall'algoritmo (anche in questo caso occorre
# usare le nostre conoscenze matematiche per controllare le uscite):
soledif(-3,3, -3,3, 25,25)
soledifc(-3,3, 3, 1e4, "brown")
soledifc(-3,-3, 3, 1001, "green4")
