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