Equazioni differenziali ordinarie

Esercizio 1 - Modello SIR di diffusione di una epidemia

Specifiche

Modello SIR di diffusione di un'epidemia

Calcolare con la function matlab ode45 la soluzione del seguente problema:

S'(t)=-aS(t)I(t) tє [0,20]

I'(t)=aS(t)I(t)-bI(t)

R'(t)=bI(t)

S(0)=199,I(0)=1;R(0)=0

Dove:

S= suscettibili di infezione

I=infetti

R= immuni(guariti)

L'epidemia si diffonde tramite l'incontro tra S ed I, il numero di persone che passa da S ad

I è proporzionale al numero di incontri secondo una costante a= costante di contagio , il

numero di guariti aumenta quando sono curati secondo una costante b=costante di

guarigione.

Porre b=0.1 e fare il grafico della soluzione, della percentuale di individui infetti, e

determinare dopo quanto tempo si verifica il picco dell'epidemia, per diversi valori di a

(0.005,0.01,0.05,0.1).

Svolgimento

Il modello SIR è stato implementato nei seguenti quattro fuction file diversi, uno per ogni valore di a:

% epidemia_SIR_1   MODELLO SIR EPIDEMIA

Implementa il modello SIR di una epidemia con la costante di contagio

pari a 0.005

Parametri di INPUT

t = vettore tale che 0<=t<=20

x = vettore dei valori iniziali

% Esempio di funzionamento:

[t,x]=ode45(@epidemia_SIR_1,[0 20],[199 1 0]);

function f = epidemia_SIR_1(t,x)

f = zeros(3,1);

a = 0.005;                    % costante di contagio

f(1) = -a*x(1)*x(2);          % possibili infetti S

f(2) = a*x(1)*x(2)-0.1*x(2);  % infetti I

f(3) = 0.1*x(2);              % immuni R

% epidemia_SIR_2 MODELLO SIR EPIDEMIA

Implementa il modello SIR di una epidemia con la costante di contagio

pari a 0.01

Parametri di INPUT

t = vettore tale che 0<=t<=20

x = vettore dei valori iniziali

% Esempio di funzionamento:

[t,x]=ode45(@epidemia_SIR_2,[0 20],[199 1 0]);

function f = epidemia_SIR_2(t,x)

f = zeros(3,1);

a = 0.01;

f(1) = -a*x(1)*x(2);

f(2) = a*x(1)*x(2)-0.1*x(2);

f(3) = 0.1*x(2);

% epidemia_SIR_3   MODELLO SIR EPIDEMIA

Implementa il modello SIR di una epidemia con la costante di contagio

pari a 0.05

Parametri di INPUT

t = vettore tale che 0<=t<=20

x = vettore dei valori iniziali

% Esempio di funzionamento:

[t,x]=ode45(@epidemia_SIR_3,[0 20],[199 1 0]);

function f = epidemia_SIR_3(t,x)

f = zeros(3,1);

a = 0.05;

f(1) = -a*x(1)*x(2);

f(2) = a*x(1)*x(2)-0.1*x(2);

f(3) = 0.1*x(2);

% epidemia_SIR_4   MODELLO SIR EPIDEMIA

Implementa il modello SIR di una epidemia con la costante di contagio

pari a 0.1

Parametri di INPUT

t = vettore tale che 0<=t<=20

x = vettore dei valori iniziali

% Esempio di funzionamento:

[t,x]=ode45(@epidemia_SIR_4,[0 20],[199 1 0]);

function f = epidemia_SIR_4(t,x)

f = zeros(3,1);

a = 0.1;

f(1) = -a*x(1)*x(2);

f(2) = a*x(1)*x(2)-0.1*x(2);

f(3) = 0.1*x(2);

E' possibile eseguire le suddette function attraverso i rispettivi script di test organizzati nel seguente menù:

scelta=menu('scegli una costante di contagio' 'a=0.005' 'a=0.01' 'a=0.05' 'a=0.1' 'fine'

if scelta == 1

run_epidemia_SIR_1;

elseif scelta==2

run_epidemia_SIR_2;

elseif scelta == 3

run_epidemia_SIR_3;

elseif scelta==4

run_epidemia_SIR_4;

elseif scelta==5

close all

return

end

clear scelta

menu_epidemia % nome dello script

Riportiamo di seguito soltanto un0 script di test ("run_epidemia_SIR_1") poichè gli altri differiscono solo nella chimata alla specifica function (riga 2 dello script):

x0=[199 1 0];

[t,x]=ode45(@epidemia_SIR_1,[0 20],x0);

% max numero di infetti

maxInfetti=max(x(:,2));

% istante di tempo corrispondente al max num di infetti

%tInf=t(find((x(:,2))==maxInfetti));

tInf=t((x(:,2))==maxInfetti);

% grafico

figure(1)

plot(t,x,tInf,maxInfetti,'O'

xlabel('tempo'

ylabel('misure'

title('Modello SIR con a=0.005'

legend('Suscettibili' 'Infetti' 'Guariti'

text(7,70,'leftarrow S' 'HorizontalAlignment' 'left' 'Color' 'b'

text(5,55,'leftarrow I' 'HorizontalAlignment' 'left' 'Color' 'g'

text(10,60,'leftarrow R' 'HorizontalAlignment' 'left' 'Color' 'r'

% tempo del picco di infetti

A=['leftarrow Picco di infetti=' num2str(tInf)];

text(tInf,maxInfetti,A,'HorizontalAlignment' 'left' 'Color' 'magenta'

Di seguito i risulati del modello SIR di una epidemia mostrati tramite grafico. Si può notare come aumentando la costante di contagio a mantenendo fissata la costante di guarigione b la percentuale di popolazione infetta aumenta con un picco sempre più vicino nel tempo):

Figure Modello SIR di una epidemia con sostante di contagio pari a 0.005

Figure Modello SIR di una epidemia con sostante di contagio pari a 0.01

Figure Modello SIR di una epidemia con sostante di contagio pari a 0.05

Figure Modello SIR di una epidemia con sostante di contagio pari a 0.1

Esercizio 2 - Problema del paracadutista

Specifiche

Calcolare con la function matlab ode45 la soluzione del seguente problema di Cauchy che

descrive il lancio di un paracadutista di massa m=80Kg da un aereo dall'altezza di 600m,

con un paracadute che si apre dopo 5 secondi.

y''(t) = -g + a(t)/m tє[0,60]

y(0)= 600

y'(0)= 0

dove (accelerazione di gravità) g= 9.81; a(t) è la resistenza dell'aria:

a(t)=K1y'(t)2 t<5 (prima dell'apertura del paracadute)

a(t)=K2y'(t)2 t>5 (dopo l'apertura del paracadute)

Testare con:

Caso 1 K1=K2=0 (caduta libera)

Caso 2 K1=1/15, K2=4/15.

A quale altezza si apre il paracadute?(caso 2)

Quando si verifica l'impatto al suolo e con che velocità ?(caso 1 e 2)

Visualizzare in un grafico le traiettorie ed i punti precedentemente calcolati

Svolgimento

Iniziamo ad esaminare il caso 1, cioè quello della caduta libera.

Il seguente function file contenenete l'implementazione del problema nel caso di caduta libera:

% paracadutista   MODELLO PARACADUTISTA

Implementa il modello modello matematico che rappresenta il lancio di un

paracadutista di massa m=80 Kg da una aereo di altezza 600 m e senza

paracadute.

Parametri di INPUT

t = vettore tale che 0<=t<=60

x = vettore dei valori iniziali

% Esempio di funzionamento:

options=odeset('events',@g);

[t,y,tfinal,yfinal]=ode45(@caduta_libera,[0 60],[600;0],options);

function pr = caduta_libera(t,y)

% caduta libera

pr = [y(2);-9.81];

Lo script che permette di testare tale funzione è "run_cadutalibera.m"

options=odeset('events',@g);

y0=[600;0];

[t,y,tfinal,yfinal]=ode45(@caduta_libera,[0 60],y0,options);

%grafico

plot(t,y(:,1),[0 tfinal],[600 0],'o'

xlabel('Tempo'

ylabel('Spazio percorso'

title('Caso 1: caduta libera'

A=[num2str(tfinal) ' = Istante di impatto al suolo rightarrow '

text(tfinal,0,A,'HorizontalAlignment' 'right' 'Color' 'black'

disp(['Istante di impatto in caduta libera:  ' num2str(tfinal)]);

disp(['Velocità di impatto in caduta libera: ' num2str(abs(yfinal(2)))]);

Eseguendo dalla console del Matlab il comando 'run run_cadutalibera' si ottiene:

Figure Grafico nel caso di caduta libera

Si nota che (output dello script):

Istante di impatto in caduta libera: 11.06

Velocità di impatto in caduta libera: 108.4988

Passiamo ora ad esaminare il caso 2, cioè quello nel quale si apre il paracadute.

Il function file contenenete l'implementazione del problema del paracadutista è "paracadutista.m":

% paracadutista   MODELLO PARACADUTISTA

Implementa il modello modello matematico che rappresenta il lancio di un

paracadutista di massa m=80 Kg da una aereo di altezza 600 m e con un

paracadute che si apre dopo 5 sec.

Parametri di INPUT

t = vettore tale che 0<=t<=60

x = vettore dei valori iniziali

% Esempio di funzionamento:

options=odeset('events',@g);

[t,y,tfinal,yfinal]=ode45(@paracadutista,[0 60],[600 0],options);

function pr = paracadutista(t,y)

% caduta con paracadute

if t<5

K = 1/15; % quando t<5

else

K = 4/15; % quando t>=5

end

pr = [y(2);-9.81+1/80*K*y(2)^2];

Mentre lo script che permette di testare tale funzione è "run_paracadutista.m":

% Caso2: calcolo altezza apertura paracadute

[t y]=ode45(@paracadutista,[0 5],[600 0]);

y1=y(:,1); % vettore dello spazio percorso

altAperParac=y1(t==5); %punto di apertura paracadute

disp(['Altezza apertura paracadute: ' num2str(altAperParac)]);

% Caso2: calcolo tempo di impatto al suolo e velocità con K1=1/15 e K2=4/15

options=odeset('events',@g);

[t,y,tfinal,yfinal]=ode45(@paracadutista,[0 60],[600 0],options);

%grafico

plot(t,y(:,1),[0 tfinal],[600 0],'o',5,altAperParac,

xlabel('Velocità'

ylabel('Spazio percorso'

title('Caso 2: caduta con paracadute'

A=['Velocità di atterraggio = ' num2str(tfinal) ' rightarrow '

text(tfinal,0,A,'HorizontalAlignment' 'right' 'Color' 'black'

B=[' leftarrow Apertura paracadute = ' num2str(altAperParac)];

text(5,altAperParac,B,'HorizontalAlignment' 'left' 'Color' 'black'

disp(['Istante di impatto con paracadute:  ' num2str(tfinal)]);

disp(['Velocità di impatto con paracadute: ' num2str(abs(yfinal(2)))]);

Eseguendo dalla console del Matlab il comando 'run run_paracadutista' si ottiene:

Figure Grafico nel caso di apertura del paracadute

Si nota inoltre che (output dello script)::

Altezza apertura paracadute: 481.3375

Istante di impatto con paracadute: 14.2739

Velocità di impatto con paracadute: 53.9514

Come ci si poteva aspettare, nel caso di caduta libera l'impatto al suolo avviene prima e con maggiore velocità rispetto alla caduta con paracadute.

Sia la funzione caduta_libera che paracadutista utilizzano la seguente function events g:

% funzione events del paracadutista

% Parametrei di INPUT

ystop = funzione di cui calcolare lo zero

isterminal = 1 se si vuole che il solver termini quando ystop p zero

direction = [] default (determina tutti dli zeri)

function [ystop,isterminal,direction] = g(t,y)

ystop = y(1);   % funzione di cui determinare lo zero

isterminal = 1;

direction = [];

Esercizio 3 - Modello di combustione

Specifiche

E' il modello di combustione di un combustibile, ad esempio, quando si accende un

fiammifero, la palla di fuoco cresce rapidamente fino a raggiungere un punto critico e

rimane di queste dimensioni fino a che non si consuma l'ossigeno. Il modello è:

y'(t)=y2(t)-y3(t) t=[0,2/d]

y(0)=d

y(t) =raggio, d= raggio iniziale(piccolo). Più è piccolo il raggio più il problema è stiff.

Risolvere con d=0.01;d=0.0001, d=0.00001, con ode45 e ode15s con RelTol=10-4 . Fare il

grafico della soluzione e determinare il numero di punti utilizzato. Cosa si osserva?

Svolgimento

Il file contenente la funzione che rappresenta il modello di combustione di un combustibile è "combustione.m":

% combustione   MODELLO di COMBUSTIONE

Implementa il modello matematico di combustione di un combustibile.

Parametri di INPUT

t = vettore tale che 0<=t<=60

x = vettore dei valori iniziali

% Esempio di funzionamento:

t=0.01;

tspan=[0 2/d];

[t,y]=ode45(@conbustione,tspan,d);

function c = conbustione(t,y)

c = y^2-y^3;

Test di funzionamento

Il test di funzionamento della funzione conbustione può essere agevolmente eseguito mediante il seguente script "test_combustione.m" il quale mette in evidenza le differenze tra il solver ode45 e ode15s:

% ODE45

d=0.01;

tspan=[0 2/d];

options=odeset('RelTol'

disp('------ ODE45 -------'

temp=cputime;

[t,y]=ode45(@conbustione,tspan,d,options);

temp=cputime-temp;

disp( d = 0.01'

disp(['Numero di punti     = ' num2str(length(y))]);

disp(['Tempo di esecuzione = ' num2str(temp)]);

d=0.0001;

tspan=[0 2/d];

temp=cputime;

[t2,y2]=ode45(@conbustione,tspan,d,options);

temp=cputime-temp;

disp( d = 0.0001'

disp(['Numero di punti     = ' num2str(length(y2))]);

disp(['Tempo di esecuzione = ' num2str(temp)]);

d=0.00001;

tspan=[0 2/d];

temp=cputime;

[t3,y3]=ode45(@conbustione,tspan,d,options);

temp=cputime-temp;

disp( d = 0.00001'

disp(['Numero di punti     = ' num2str(length(y3))]);

disp(['Tempo di esecuzione = ' num2str(temp)]);

% grafico

figure(1)

semilogx(t,y,'-g',t2,y2,'-r',t3,y3,'-b'

xlabel('Tempo'

ylabel('Dimensione palla di fuoco'

title('Combustione con ode45'

legend('d=0.01' 'd=0.0001' 'd=0.00001'

% ODE15S

d=0.01;

tspan=[0 2/d];

options=odeset('RelTol'

disp( '

disp('------ ODE15S -------'

temp=cputime;

[t,y]=ode15s(@conbustione,tspan,d,options);

temp=cputime-temp;

disp( d = 0.01'

disp(['Numero di punti     = ' num2str(length(y))]);

disp(['Tempo di esecuzione = ' num2str(temp)]);

d=0.0001;

tspan=[0 2/d];

temp=cputime;

[t2,y2]=ode15s(@conbustione,tspan,d,options);

temp=cputime-temp;

disp( d = 0.0001'

disp(['Numero di punti     = ' num2str(length(y2))]);

disp(['Tempo di esecuzione = ' num2str(temp)]);

d=0.00001;

tspan=[0 2/d];

temp=cputime;

[t3,y3]=ode15s(@conbustione,tspan,d,options);

temp=cputime-temp;

disp( d = 0.00001'

disp(['Numero di punti     = ' num2str(length(y3))]);

disp(['Tempo di esecuzione = ' num2str(temp)]);

% grafico

figure(2)

semilogx(t,y,'-g',t2,y2,'-r',t3,y3,'-b'

xlabel('Tempo'

ylabel('Dimensione palla di fuoco'

title('Combustione con ode15s'

legend('d=0.01' 'd=0.0001' 'd=0.00001'

Eseguendo tramite console Matlab il suddetto script, si ottengono i grafici in figura1 e figura2 ed il seguente output:

>> run run_conbustione

------ ODE45 -------

d = 0.01

Numero di punti   = 185

Tempo di esecuzione = 0.0468

d = 0.0001

Numero di punti     = 12161

Tempo di esecuzione = 0.546

d = 0.00001

Numero di punti     = 120701

Tempo di esecuzione = 5.9592

------ ODE15S -------

d = 0.01

Numero di punti     = 86

Tempo di esecuzione = 0.0312

d = 0.0001

Numero di punti     = 141

Tempo di esecuzione = 0.078001

d = 0.00001

Numero di punti     = 140

Tempo di esecuzione = 0.078001:

Figure Grafico della combustione con ode45 Figure Grafico della combustione con ode15s

Si può notare come per valori di d decrescenti il solver ode45 impieghi sia più tempo di esecuzione che punti per calcolare la funzione di combustione. Ciò evidenzia come tael problema sia di tipo stiff.

Esercizio 4 - Algoritmo RKF45

Specifiche

Function file che implementa :

Algoritmo RKF45 a passo variabile con il parametro TOL opzionale.

Svolgimento

Il file contenente l'implementazione dell'algoritmo RKF45 è "rkf45.m":

% RKF45   algoritmo RKF45

La funzione RFK45 risolve un sistema di equazioni differenziali del tipo

y'=f(t,y(t)) utilizzando le formule di Runge-Kutta Fehlberg di ordine 4 e 5

nell'intervallo [t,T] con passo di discretizzazione adattativo, valore

iniziale y0 e tolleranza assegnata TOL.

% Paremtri di INPUT

f = function_handle che descrive il sistema come yi'=fi(t,y(t)).

Deve avere 2 parametri in ingresso t ed y

Tspan = intervallo di ricerca della soluzione

y0 = vettore delle condizioni iniziali

TOL (opzionale) = tolleranza definita dall'utente; di default vale 1e-6

a,b,c (opzionali) = costanti definite dall'utente; di default valgono

rispettivamente 5, 20000 e 100.

verbose (opzionale) = stampa a video messaggi

Parametri di OUTPUT

t : vettore dei punti generati dalla discretizzazione

y : matrice della soluzione approssimata in cui la colonna i-esima rappresenta la soluzione dell'i-esima equazione differenziale valutata in t

function [t,y] = rkf45(f,Tspan,y0,TOL,a,b,c,verbose)

% controllo numero parametri di input

if nargin <3

error('Parametri obbligatori mancanti'

elseif nargin == 3

% valori di default

TOL = 1.e-6;

a = 5;

b = 20000;

c = 100;

verbose = 0;

elseif nargin == 4

a = 5;

b = 20000;

c = 100;

verbose = 0;

elseif nargin == 5

b = 20000;

c = 100;

verbose = 0;

elseif nargin == 6

c = 100;

elseif nargin == 7

verbose = 0;

end

% controllo f

if (~(isobject(f) || isa(f, 'function_handle') || ischar(f)))

error('f deve essere un function_handle, una funzione inline o una stringa'

end

% controllo il vettore y0

outcheck = checkVector(y0);

if ~isempty(outcheck)

error(['Parametro y0 invalido: ', outcheck]);

end

% controllo Tspan

if isnumeric(Tspan)~=1 || isvector(Tspan)~=1 || length(Tspan)~=2

error('Tspan deve essere un vettore numerico di 2 componenti'

elseif Tspan(1)>=Tspan(2)

error('Intervallo [t,T] non valido'

end

% controllo validità t

outcheck = isPositiveNumber(Tspan(1));

if ~isempty(outcheck)

error(['Parametro t invalido: ', outcheck]);

end

% controllo validità T

outcheck = isPositiveNumber(Tspan(2));

if ~isempty(outcheck)

error(['Parametro T invalido: ', outcheck]);

end

% controllo validità TOL

outcheck = isPositiveNumber(TOL);

if ~isempty(outcheck)

error(['Parametro TOL invalido: ', outcheck]);

end

% correggo la tolleranza TOL

TOL = max(TOL,eps);

if ischar(a) || ischar(b) || ischar(c) || isnan(a) || isnan(b) || isnan(c)

|| isinf(a) || isinf(b) || isinf(c)

error('Le costanti a b c devono essere dei numeri'

end

t = Tspan(1);

T = Tspan(2);

hmax = (T-t)/a;

hmin = (T-t)/b;

h = (T-t)/c;

% visualizzazione messaggi utente

if verbose == 1

disp(['hmin = ' num2str(hmin)]);

disp(['hmax = ' num2str(hmax)]);

disp(['h = ' num2str(h)]);

disp(['a = ' num2str(a)]);

disp(['b = ' num2str(b)]);

disp(['c = ' num2str(c)]);

disp(['TOL = ' num2str(TOL)]);

end

i = 1;

t(i) = t;

y(1,:) = y0;

while t(i)<T && h>=hmin

if t(i)+h > T

h = T-t(i);

end

% calcolo Ki, i=16

K1 = h*(feval(f,t(i),y(i,:)));

K2 = h*(feval(f,t(i)+h/4,y(i,:)+K1'/4));

K3 = h*(feval(f,t(i)+3*h/8,y(i,:)+(3*K1'+9*K2')/32));

K4 = h*(feval(f,t(i)+12*h/13,y(i,:)+(1932*K1'-7200*K2'+7296*K3')/2197));

K5 = h*(feval(f,t(i)+h,y(i,:)+439*K1'/216-8*K2'+3680*K3'/513-845*K4'/4104));

K6 = h*(feval(f,t(i)+h/2,y(i,:)-8*K1'/27+2*K2'-3544*K3'/2565+1859*K4'/4104-11*K5'/40));

% calcolo errore R

R = max(abs((1/360)*K1-(128/4275)*K3-(2197/75240)*K4+(1/50)*K5+(2/55)*K6));

if verbose == 1

disp(['K1 = ' num2str(K1)]);

disp(['K2 = ' num2str(K2)]);

disp(['K3 = ' num2str(K3)]);

disp(['K4 = ' num2str(K4)]);

disp(['K5 = ' num2str(K5)]);

disp(['K6 = ' num2str(K6)]);

disp(['R  = ' num2str(R)]);

end

% controllo dell'accuratezza

if R<TOL*h     

t(i+1) = t(i)+h;

% calcolo nuova approssimazione y(k+1)

y(i+1,:) = y(i,:)+((25/216)*K1'+(1408/2565)*K3'+(2197/4104)*K4'-(K5'/5));

i = i+1;

end

% calcolo nuovo passo

if R~=0

q = 0.84*(TOL*h/R)^0.25;

h = min(hmax,q*h);

end

end

% condizione che si verifica quando h<hmin

if t<T

warning(1'probabile singolarità'

end

Test di funzionamento

Il test sull'algoritmo rkf45 verrà effettuato tramite la funzione epidemia_SIR_1 presentata nell'esercizio 1 così da poter confrontare i risulati di rkf45 con quelli calcolati mediante ode45.

Lo script seguente ("test_rkf45") permette di eseguire facilemnte la funzione rkf45:

% test con il modello SIR

b=0.1;

a=0.005;

% vettore valori iniziali

x0=[199 1 0];

[t,x]=rkf45(@epidemia_SIR_1,[0 20],x0);

%grafico soluzione

figure;

plot(t,x);

text=['Soluzione SIR con a=', num2str(a)];

title(text);

legend('Suscettibili' 'Infetti' 'IGuariti'

xlabel('Giorni'

ylabel('Popolazione'

L'output generato è il seguente:

Figure Grafico del modello SIR calcolato con rkf45

Come si può notare il grafico di figura 9 è identico a Figure 1 dell'esercizio 1.

Passiamo infine a controllare la robustezza della funzione rkf45 nei confronti dei parametri di ingresso.

A tale scopo consideriamo la seguente funzione di ingresso

>> f=inline('-5*y','t','y')

f =

Inline function:

f(t,y) = -5*y

• Omissione parametri di input

>> [t,y]=rkf45(f)

??? Error using ==> rkf45 at 25

Parametri obbligatori mancanti

>> [t]=rkf45(f)

??? Error using ==> rkf45 at 25

Parametri obbligatori mancanti

• Vettore y0 contiene caratteri

>> [t]=rkf45(f,[1 2],['ciao'])

??? Error using ==> rkf45 at 56

Parametro y0 invalido: deve contenere solo numeri

• Il vettore y0 contiene Nan

>> [t]=rkf45(f,[1 2],[nan])

??? Error using ==> rkf45 at 56

Parametro y0 invalido: non deve contenere valori NaN

• Il vettore y0 contiene Inf

>> [t]=rkf45(f,[1 2],[inf])

??? Error using ==> rkf45 at 56

Parametro y0 invalido: non deve contenere valori Inf

• Valori dt Tspan errati

>> [t]=rkf45(f,[3 1],10^-6)

??? Error using ==> rkf45 at 63

Intervallo [t,T] non valido

• Tspan contiene caratteri

>> [t]=rkf45(f,['1' 2],10^-6)

??? Error using ==> rkf45 at 61

Tspan deve essere un vettore numerico di 2 componenti

• Tspan è un vettore con più di 3 elementi

>> [t]=rkf45(f,[1 2 3],10^-6)

??? Error using ==> rkf45 at 61

Tspan deve essere un vettore numerico di 2 componenti

• Tspan contiene Inf

>> [t]=rkf45(f,[inf 3],[1])

??? Error using ==> rkf45 at 63

Intervallo [t,T] non valido

• La tolleranza TOL è un carattere

>> [t]=rkf45(f,[1 3],[1], 'tol')

??? Error using ==> rkf45 at 79

Parametro TOL invalido: non deve contenere caratteri

• La tolleranza TOL è un numero negativo

>> [t]=rkf45(f,[1 3],[1], -1)

??? Error using ==> rkf45 at 79

Parametro TOL invalido: deve essere positivo

• La tolleranza TOL contiene NaN

>> [t]=rkf45(f,[1 3],[1], nan)

??? Error using ==> rkf45 at 79

Parametro TOL invalido: non deve essere NaN