Considera l'equazione differenziale  y'(x) = y(x)^2.  Sia y(0) = 1.  Vogliamo calcolare sperimentalmente quanto valgono y(0.4), y(0.8), y(1).  Usando R (vedi) o WolframAlpha, traccia il campo direzionale della equazione in  [-1/2, 2]×[-2, 6]  e traccia il grafico della soluzione tale che y(0)=1.  Prova a trovare i valori delle soluzioni arrotondati a 7 cifre.  Che cosa ottieni?  Verifica, poi, che  y(x) = 1/(1−x)  è la soluzione esatta dell'equazione.

Con R. Si vede dalle uscite, sotto riportate che y(0.4) = 5/3, y(0.8) = 5 mentre y(1) non è definito:  al raddoppiare di "n" nei primi due casi la variazione di dimezza  (per cui posso approssimare il valore aggiungendo all'ultima uscita la variazione rispetto alla precedente, in quanto 1/2+1/4+1/8+… = 1),  nell'ultimo caso la variazione cresce. Possiamo controllare la cosa utilizzando la soluzione proposta (è facile verificare che y' = y²):  1/(1−x) non è definito per x=1  (mentre per x=0.4 e x=0.8 ritrovo i valori trovati sperimentalmente).

Dy = function(x,y) y^2
BF=4; HF=3
Plane(-0.5,2, -2,6)
diredif(-0.5,2, -2,6, 19,19)
soledif(0,1, 2,1e5,"brown"); soledif(0,1, -1,1e5,"seagreen")
A=soledifv(0,1, 0.4, 1e5); A
# 1.666661
B=A; A=soledifv(0,1, 0.4, 2e5); c(A,A-B)
# 1.666664e+00 2.837866e-06
B=A; A=soledifv(0,1, 0.4, 4e5); c(A,A-B)
# 1.666665e+00 1.418947e-06   Ad ogni passo la variazione si dimezza. Quindi aggiungo
A+A-B                         # ad A la variazione precedente in quanto 1/2+1/4+... = 1
# 1.666667
fraction(A+A-B)
# 5/3
A=soledifv(0,1, 0.8, 1e5); A
# 4.999678
B=A; A=soledifv(0,1, 0.8, 2e5); c(A,A-B)
# 4.9998390654 0.0001609162
B=A; A=soledifv(0,1, 0.8, 4e5); c(A,A-B)
# 4.999919530 0.000080465    Ad ogni passo la variazione si dimezza. Quindi aggiungo
A+A-B                        # ad A la variazione precedente
# 5
A=soledifv(0,1, 1, 1e5); A
# 10823.81
B=A; A=soledifv(0,1, 1, 2e5); c(A,A-B)
# 20265.285  9441.477
B=A; A=soledifv(0,1, 1, 4e5); c(A,A-B)
# 38085.36 17820.07                    non converge! in 1 c'e' asintoto verticale
f = function(x) 1/(1-x)
deriv(f,"x")
# 1/(1 - x)^2  OK: il quadrato di f(x)

Con WolframAlpha

slope field of dy/dx = y^2 for -1/2 < x < 2, -2 < y < 6

dy/dx = y^2, y(0)=1

plot y = 1/(1-x), y = 1, -1/2 < x < 2, -2 < y < 6

Per altri commenti: Modelli differenziali.