ode45

Risolvere equazioni differenziali non rigide - metodo di ordine medio

Sintassi

[t,y] = ode45(odefun,tspan,y0)

[t,y] = ode45(odefun,tspan,y0,options)

[t,y,te,ye,ie] = ode45(odefun,tspan,y0,options)

sol = ode45(___)

Descrizione

[t,y] = ode45(odefun,tspan,y0), dove tspan = [t0 tf] integra il sistema di equazioni differenziali $$y\text{'}=f\left(t,y\right)$$ da t0 a tf con condizioni iniziali y0. Ciascuna riga dell'array della soluzione y corrisponde a un valore restituito nel vettore colonna t.

Tutti i risolutori ODE di MATLAB^® possono risolvere sistemi di equazioni della forma $$y\text{'}=f\left(t,y\right)$$ o problemi che coinvolgono una matrice di massa $$M\left(t,y\right)y\text{'}=f\left(t,y\right)$$. Tutti i risolutori utilizzano sintassi simili. Il risolutore ode23s può risolvere solo problemi con una matrice di massa se la suddetta è costante. ode15s e ode23t possono risolvere problemi con una matrice di massa singolare, note come equazioni differenziali-algebriche (DAE). Specificare la matrice di massa utilizzando l'opzione Mass di odeset.

ode45 è un risolutore ODE versatile ed è il primo risolutore che si dovrebbe provare per la maggior parte dei problemi. Tuttavia, se il problema è rigido o richiede un'elevata precisione, esistono altri risolutori ODE che potrebbero essere più adatti al problema. Per maggiori informazioni vedere Scelta di un risolutore ODE.

[t,y] = ode45(odefun,tspan,y0,options) utilizza inoltre le impostazioni di integrazione definite da options, che è un argomento creato utilizzando la funzione odeset. Ad esempio, utilizzare le opzioni AbsTol e RelTol per specificare le tolleranze di errore assolute e relative oppure l'opzione Mass per fornire una matrice di massa.

[t,y,te,ye,ie] = ode45(odefun,tspan,y0,options) inoltre individua i punti in cui le funzioni di (t,y), chiamate funzioni evento, sono pari a zero. Nell'output, te è il tempo dell'evento, ye è la soluzione al tempo dell'evento ie è l'indice dell'evento attivato.

Per ciascuna funzione evento, specificare se l'integrazione deve terminare a zero e se la direzione del passaggio per lo zero è rilevante. A tal fine, impostare la proprietà 'Events' su una funzione, ad esempio myEventFcn o @myEventFcn, e creare una funzione corrispondente: [value,isterminal,direction] = myEventFcn(t,y). Per maggiori informazioni, vedere ODE Event Location.

sol = ode45(___) restituisce una struttura che è possibile utilizzare con deval per valutare la soluzione in qualsiasi punto dell'intervallo [t0 tf]. È possibile utilizzare qualsiasi combinazione di argomenti di input nelle sintassi precedenti.

Esempi

ODE con componente di soluzione singolo

Le ODE semplici che hanno un componente di soluzione singolo possono essere specificate come funzione anonima nella chiamata al risolutore. La funzione anonima deve accettare due input (t,y), anche se uno degli input non è utilizzato nella funzione.

Risolvere l'ODE

$$y^{\prime }=2t.$$

Specificare un intervallo di tempo di [0 5] e la condizione iniziale y0 = 0.

tspan = [0 5];
y0 = 0;
[t,y] = ode45(@(t,y) 2*t, tspan, y0);

Tracciare la soluzione.

plot(t,y,'-o')

Risoluzione dell'equazione non rigida

L'equazione di van der Pol è un'ODE di secondo ordine

$$y_1^{\prime \prime }-\mu \left(1-y_1^2\right)y_1^{\prime }+y_1=0,$$

dove $$\mu >0$$ è un parametro scalare. Riscrivere questa equazione come un sistema di ODE di primo ordine effettuando la sostituzione $$y_1^{\prime }=y_2$$. Il sistema risultante di ODE di primo ordine è

$$\begin{array}{cl}{y}_1^{\prime }& =y_2\\ y_2^{\prime }& =\mu (1-y_1^2)y_2-y_1.\end{array}$$

Il file di funzione vdp1.m rappresenta l'equazione di van der Pol utilizzando $$\mu =1$$. Le variabili $$y_1$$ e $$y_2$$ sono gli elementi y(1) e y(2) di un vettore a due elementi dydt.

function dydt = vdp1(t,y)
%VDP1  Evaluate the van der Pol ODEs for mu = 1
%
%   See also ODE113, ODE23, ODE45.

%   Jacek Kierzenka and Lawrence F. Shampine
%   Copyright 1984-2014 The MathWorks, Inc.

dydt = [y(2); (1-y(1)^2)*y(2)-y(1)];

Risolvere l'ODE utilizzando la funzione ode45 nell'intervallo di tempo [0 20] con valori iniziali [2 0]. L'output risultante è un vettore colonna di punti temporali t e un array di soluzioni y. Ciascuna riga in y corrisponde a un tempo restituito nella riga corrispondente di t. La prima colonna di y corrisponde a $$y_1$$, mentre la seconda colonna corrisponde a $$y_2$$.

[t,y] = ode45(@vdp1,[0 20],[2; 0]);

Tracciare le soluzioni per $$y_1$$ e $$y_2$$ rispetto a t.

plot(t,y(:,1),'-o',t,y(:,2),'-o')
title('Solution of van der Pol Equation (\mu = 1) with ODE45');
xlabel('Time t');
ylabel('Solution y');
legend('y_1','y_2')

Passaggio di ulteriori parametri alla funzione ODE

ode45 è utilizzabile solo con funzioni che utilizzano due argomenti di input t e y. Tuttavia, è possibile passare ulteriori parametri definendoli al di fuori della funzione e inserendoli quando si specifica l'handle della funzione.

Risolvere l'ODE

$$y^{\prime \prime }=\frac{A}{B}ty.$$

Riscrivendo l'equazione come sistema di primo ordine si ottiene

$$\begin{array}{cl}{y}_1^{\prime }& =y_2\\ y_2^{\prime }& =\frac{A}{B}t{y}_1.\end{array}$$

odefcn, una funzione locale inclusa alla fine di questo esempio, rappresenta questo sistema di equazioni come una funzione che accetta quattro argomenti di input: t, y, A e B.

function dydt = odefcn(t,y,A,B)
  dydt = zeros(2,1);
  dydt(1) = y(2);
  dydt(2) = (A/B)*t.*y(1);
end

Risolvere l'ODE utilizzando ode45. Specificare l'handle della funzione in modo che passi i valori predefiniti per A e B a odefcn.

A = 1;
B = 2;
tspan = [0 5];
y0 = [0 0.01];
[t,y] = ode45(@(t,y) odefcn(t,y,A,B), tspan, y0);

Tracciare i risultati.

plot(t,y(:,1),'-o',t,y(:,2),'-.')

function dydt = odefcn(t,y,A,B)
  dydt = zeros(2,1);
  dydt(1) = y(2);
  dydt(2) = (A/B)*t.*y(1);
end

Risoluzione dell'ODE con più condizioni iniziali

Per i sistemi di ODE semplici con una sola equazione, è possibile specificare y0 come vettore contenente più condizioni iniziali. Questa tecnica crea un sistema di equazioni indipendenti tramite espansione scalare, una per ciascun valore iniziale, e ode45 risolve il sistema per produrre risultati per ciascun valore iniziale.

Creare una funzione anonima per rappresentare l'equazione $f\left(t,y\right)=-2y+2\cos \left(t\right)\sin \left(2t\right)$. La funzione deve accettare due input per t e y.

yprime = @(t,y) -2*y + 2*cos(t).*sin(2*t);

Creare un vettore di condizioni iniziali diverse nell'intervallo $\left[-5,5\right]$.

y0 = -5:5; 

Risolvere l'equazione per ciascuna condizione iniziale nell'intervallo di tempo $\left[0,3\right]$ utilizzando ode45.

tspan = [0 3];
[t,y] = ode45(yprime,tspan,y0);

Tracciare i risultati.

plot(t,y)
grid on
xlabel('t')
ylabel('y')
title('Solutions of y'' = -2y + 2 cos(t) sin(2t), y(0) = -5,-4,...,4,5','interpreter','latex')

Questa tecnica risulta utile per risolvere ODE semplici con condizioni iniziali diverse. Tuttavia, questa tecnica presenta anche alcuni compromessi:

• Non è possibile risolvere sistemi di equazioni con più condizioni iniziali. La tecnica funziona solo quando si risolve un'equazione con più condizioni iniziali.

• Il passo temporale scelto dal risolutore a ogni tempo si basa sull'equazione nel sistema che richiede il passo più piccolo. Questo significa che il risolutore può eseguire piccoli passi per soddisfare l'equazione per una condizione iniziale, ma le altre equazioni, se risolte separatamente, utilizzerebbero passi di dimensioni diverse. Ciononostante, risolvere più condizioni iniziali contemporaneamente è generalmente più veloce che risolvere le equazioni separatamente utilizzando un for loop.

Per maggiori informazioni su questa tecnica, vedere Solve System of ODEs with Multiple Initial Conditions.

ODE con termini dipendenti dal tempo

Si consideri la seguente ODE con parametri dipendenti dal tempo

La condizione iniziale è . La funzione f(t) è definita dal vettore n x 1 f valutato ai tempi ft. La funzione g(t) è definita dal vettore m x 1 g valutato ai tempi gt.

Creare i vettori f e g.

ft = linspace(0,5,25);
f = ft.^2 - ft - 3;

gt = linspace(1,6,25);
g = 3*sin(gt-0.25);

Scrivi una funzione denominata myode che interpoli f e g per ottenere il valore dei termini dipendenti dal tempo al tempo specificato. Salvare la funzione nella cartella attuale per eseguire il resto dell'esempio.

La funzione myode accetta ulteriori argomenti di input per valutare l'ODE a ogni passo temporale, ma ode45 utilizza solo i primi due argomenti di input t e y.

function dydt = myode(t,y,ft,f,gt,g)
f = interp1(ft,f,t); % Interpolate the data set (ft,f) at time t
g = interp1(gt,g,t); % Interpolate the data set (gt,g) at time t
dydt = -f.*y + g; % Evaluate ODE at time t

Risolvere l'equazione nell'intervallo di tempo [1 5] utilizzando ode45. Specificare la funzione utilizzando un handle della funzione in modo che ode45 utilizzi solo i primi due argomenti di input di myode. Ridurre inoltre la restrittività delle soglie di errore utilizzando odeset.

tspan = [1 5];
ic = 1;
opts = odeset('RelTol',1e-2,'AbsTol',1e-4);
[t,y] = ode45(@(t,y) myode(t,y,ft,f,gt,g), tspan, ic, opts);

Tracciare la soluzione y come funzione dei punti temporali t.

plot(t,y)

Valutazione ed estensione della struttura della soluzione

L'equazione di van der Pol è un'ODE di secondo ordine

$$y_1^{\prime \prime }-\mu \left(1-y_1^2\right)y_1^{\prime }+y_1=0.$$

Risolvere l'equazione di van der Pol con $$\mu =1$$ utilizzando ode45. La funzione vdp1.m è fornita con MATLAB® e codifica le equazioni. Specificare un unico output per restituire una struttura contenente informazioni sulla soluzione, come il risolutore e i punti di valutazione.

tspan = [0 20];
y0 = [2 0];
sol = ode45(@vdp1,tspan,y0)

sol = struct with fields:
     solver: 'ode45'
    extdata: [1×1 struct]
          x: [0 1.0048e-04 6.0285e-04 0.0031 0.0157 0.0785 0.2844 0.5407 0.8788 1.4032 1.8905 2.3778 2.7795 3.1285 3.4093 3.6657 3.9275 4.2944 4.9013 5.3506 5.7998 6.2075 6.5387 6.7519 6.9652 7.2247 7.5719 8.1226 8.6122 9.1017 9.5054 … ] (1×60 double)
          y: [2×60 double]
      stats: [1×1 struct]
      idata: [1×1 struct]

Utilizzare linspace per generare 250 punti nell'intervallo [0 20]. Valutare la soluzione in questi punti utilizzando deval.

x = linspace(0,20,250);
y = deval(sol,x);

Tracciare il primo componente della soluzione.

plot(x,y(1,:))

Estendere la soluzione a $$t_f=35$$ utilizzando odextend e aggiungere il risultato al grafico originale.

sol_new = odextend(sol,@vdp1,35);
x = linspace(20,35,350);
y = deval(sol_new,x);
hold on
plot(x,y(1,:),'r')

Argomenti di input

odefun — Funzioni da risolvere

handle della funzione

Funzioni da risolvere, specificate come handle della funzione che definisce le funzioni da integrare.

La funzione dydt = odefun(t,y), per uno scalare t e un vettore colonna y, deve restituire un vettore colonna dydt del tipo di dato single o double che corrisponda a $$f\left(t,y\right)$$. odefun deve accettare entrambi gli argomenti di input t e y, anche se uno degli argomenti non è utilizzato nella funzione.

Ad esempio, per risolvere $$y\text{'}=5y-3$$, utilizzare la funzione:

function dydt = odefun(t,y)
dydt = 5*y-3;
end

Per un sistema di equazioni, l'output di odefun è un vettore. Ogni elemento del vettore è il valore calcolato del lato destro di un'equazione. Ad esempio, si consideri il sistema di due equazioni

$$\begin{array}{l}y{\text{'}}_1=y_1+2y_2\\ y{\text{'}}_2=3y_1+2y_2\end{array}$$

Una funzione che calcola il valore del lato destro di ciascuna equazione a ogni passo temporale è

function dydt = odefun(t,y)
dydt = zeros(2,1);
dydt(1) = y(1)+2*y(2);
dydt(2) = 3*y(1)+2*y(2);
end

Per informazioni su come fornire ulteriori parametri alla funzione odefun, vedere Parameterizing Functions.

Esempio @myFcn

Tipi di dati: function_handle

tspan — Intervallo di integrazione

vettore

Intervallo di integrazione, specificato come vettore. Come minimo, tspan deve essere un vettore a due elementi [t0 tf] che specifica i tempi iniziale e finale. Per ottenere le soluzioni a tempi specifici compresi tra t0 e tf, utilizzare un vettore più lungo della forma [t0,t1,t2,...,tf]. Gli elementi in tspan devono essere tutti crescenti o tutti decrescenti.

Il risolutore impone le condizioni iniziali specificate da y0 al tempo iniziale tspan(1), quindi integra da tspan(1) a tspan(end):

• Se tspan presenta due elementi [t0 tf], il risolutore restituisce la soluzione valutata a ogni passo di integrazione interno all'intervallo.

• Se tspan presenta più di due elementi [t0,t1,t2,...,tf], il risolutore restituisce la soluzione valutata nei punti specificati. Tuttavia, il risolutore non passa esattamente su ciascun punto specificato in tspan. Il risolutore utilizza invece i propri passi interni per calcolare la soluzione, quindi valuta la soluzione nei punti richiesti in tspan. Le soluzioni prodotte nei punti specificati hanno lo stesso ordine di precisione delle soluzioni calcolate a ogni passo interno.

  La specifica di diversi punti intermedi ha un effetto minimo sull'efficienza computazionale, ma può influire sulla gestione della memoria nei sistemi di grandi dimensioni.

I valori di tspan sono utilizzati dal risolutore per calcolare i valori appropriati per InitialStep e MaxStep:

• Se tspan contiene diversi punti intermedi [t0,t1,t2,...,tf], i punti specificati forniscono un'indicazione della scala del problema, che può influire sul valore di InitialStep utilizzato dal risolutore. Pertanto, la soluzione ottenuta dal risolutore potrebbe essere diversa a seconda che tspan venga specificato come vettore a due elementi o come vettore con punti intermedi.

• I valori iniziale e finale in tspan sono utilizzati per calcolare la dimensione massima del passo MaxStep. Pertanto, la modifica dei valori iniziali o finali in tspan può indurre il risolutore a utilizzare una sequenza di passi diversa, che potrebbe modificare la soluzione.

Esempio [1 10]

Esempio [1 3 5 7 9 10]

Tipi di dati: single | double

y0 — Condizioni iniziali

vettore

Condizioni iniziali, specificate come vettore. y0 deve avere la stessa lunghezza del vettore di output di odefun, in modo che y0 contenga una condizione iniziale per ciascuna equazione definita in odefun.

Tipi di dati: single | double

options — Struttura opzioni

array di struttura

Struttura opzione, specificata come array di struttura. Utilizzare la funzione odeset per creare o modificare la struttura options. Per un elenco delle opzioni compatibili con ciascun risolutore, vedere Riepilogo delle opzioni ODE.

Esempio options = odeset('RelTol',1e-5,'Stats','on','OutputFcn',@odeplot) specifica una tolleranza di errore relativo di 1e-5, attiva la visualizzazione delle statistiche del risolutore e specifica la funzione di output @odeplot per tracciare la soluzione man mano che viene calcolata.

Tipi di dati: struct

Argomenti di output

t — Punti di valutazione

vettore colonna

Punti di valutazione, restituiti come vettore colonna.

• Se tspan contiene due elementi [t0 tf], t contiene i punti di valutazione interni utilizzati per eseguire l'integrazione.

• Se tspan contiene più di due elementi, t equivale a tspan.

y — Soluzioni

array

Soluzioni, restituite come array. Ciascuna riga in y corrisponde alla soluzione calcolata al valore restituito nella riga corrispondente di t.

te — Tempo degli eventi

vettore colonna

Tempo degli eventi, restituito come vettore colonna. I tempi degli eventi in te corrispondono alle soluzioni restituite in ye, mentre ie specifica quale evento si è verificato.

ye — Soluzione al tempo degli eventi

array

Soluzione al tempo degli eventi, restituita come array. I tempi degli eventi in te corrispondono alle soluzioni restituite in ye, mentre ie specifica quale evento si è verificato.

ie — Indice della funzione dell'evento trigger

vettore colonna

Indice della funzione dell'evento trigger, restituito come vettore colonna. I tempi degli eventi in te corrispondono alle soluzioni restituite in ye, mentre ie specifica quale evento si è verificato.

sol — Struttura per la valutazione

array di struttura

Struttura per la valutazione, restituita come array di struttura. Utilizzare questa struttura con la funzione deval per valutare la soluzione in qualsiasi punto dell'intervallo [t0 tf]. L'array di struttura sol include sempre questi campi:

Campo della struttura	Descrizione
sol.x	Vettore riga dei passi scelti dal risolutore.
sol.y	Soluzioni. Ciascuna colonna sol.y(:,i) contiene la soluzione al tempo sol.x(i).
sol.solver	Nome del risolutore.

Inoltre, se si specifica l'opzione Events di odeset e vengono rilevati eventi, sol include anche questi campi:

Campo della struttura	Descrizione
sol.xe	Punti in cui si sono verificati eventi. sol.xe(end) contiene il punto esatto di un evento terminale, se presente.
sol.ye	Soluzioni corrispondenti agli eventi in sol.xe.
sol.ie	Indici nel vettore restituito dalla funzione specificata nell'opzione Events. I valori indicano quale evento è stato rilevato dal risolutore.

Algoritmi

ode45 si basa su una formula di Runge-Kutta esplicita (4,5), la coppia Dormand-Prince. Si tratta di un risolutore a passo singolo: per calcolare y(tn), necessita solo della soluzione al punto temporale immediatamente precedente y(tn-1) [1], [2].

Riferimenti

[1] Dormand, J. R. and P. J. Prince, “A family of embedded Runge-Kutta formulae,” J. Comp. Appl. Math., Vol. 6, 1980, pp. 19–26.

[2] Shampine, L. F. and M. W. Reichelt, “The MATLAB ODE Suite,” SIAM Journal on Scientific Computing, Vol. 18, 1997, pp. 1–22.

Funzionalità estese

Generazione di codice C/C++

Genera codice C e C++ con MATLAB® Coder™.

Note su utilizzo e limitazioni:

• È necessario fornire punti di valutazione e soluzioni come argomenti di output.

• I tipi di input devono essere omogenei: tutti doppi o tutti singoli.

• Il supporto per le dimensioni variabili deve essere abilitato. La generazione di codice richiede l'allocazione dinamica della memoria quando tspan presenta due elementi o si utilizzano funzioni evento.

Generazione di codice GPU

Genera codice CUDA® per GPU NVIDIA® con GPU Coder™.

Consultare le note su utilizzo e limitazioni nella sezione Generazione di codice C/C++. Le stesse note su utilizzo e limitazioni si applicano alla generazione di codice GPU.

Cronologia versioni

Introduzione prima di R2006a

Vedi anche

ode23 | ode78 | ode89 | ode113 | ode15s | odeset | odeget | deval | odextend

Argomenti

• Scelta di un risolutore ODE
• Riepilogo delle opzioni ODE
• Solve Nonstiff ODEs
• Funzioni anonime
• Troubleshoot Common ODE Problems