Risoluzione Equazioni Differenziali Ordinarie

[Gasso1]

Mi sono imbattuto in un esercizio nel quale mi chiede di utilizzare il metodo di Eulero implicito nel caso in cui $t_0$=0, T=10, N=100, $y_0$=$(1,2)^T$ e, se y =$ (x_1,x_2)^T$, f(t,y) = $(-10^3x_1+sinx_1cosx_2, -2x_2+sinx_1cosx_2)^T$.

Se non mi sbaglio il metodo di Eulero implicito è questo: $y_n = y_(n-1)+hf_n$, per n = 1,2....,N = $T-t_0/h$
Come posso applicarlo al mio esercizio?

[apatriarca]

Che linguaggio devi usare per implementare il metodo di eulero implicito? Matlab?

[vict85]

Cosa non capisci di quello che devi fare?

Parti da $t=0$ e esegui ricorsivamente il metodo... Sei sicuro che la formula per N sia corretta?


$T$ l'hai usato sia per la trasposta che per il tempo è meglio se cambi qualcosa...

Usa le formule su tutto e non solo sulle parti più matematiche.

[Gasso1]

"apatriarca":Che linguaggio devi usare per implementare il metodo di eulero implicito? Matlab?

Si devo usare Matlab!

[apatriarca]

Devi per prima cosa calcolarti h, utilizzando T e t0 come nel seguente codice (la definizione di N è sbagliata nel tuo codice):

h = (T - t0)/N;

A questo punto devi creare una matrice Nx2 in cui hai inizializzato la prima colonna a [1; 2] e fai partire un ciclo in cui calcoli il successivo valore di Y risolvendo l'equazione implicita che hai scritto. Per farlo puoi ad esempio usare fzero. Questa è l'idea generale del codice, prova a implementarlo e se hai difficoltà provo a vedere di fornirti qualche delucidazione in più.

[apatriarca]

Ho appena notato che nella definizione della funzione, \(f(t,y)\) non dipende da \(t\), ma solo da \(y = (x_1, x_2)\). Si tratta inoltre di una funzione a 3 variabili e quindi fzero non dovrebbe essere utilizzabile. È corretto? Il testo del tuo problema non è chiarissimo.

[Gasso1]

"apatriarca":Ho appena notato che nella definizione della funzione, \(f(t,y)\) non dipende da \(t\), ma solo da \(y = (x_1, x_2)\). Si tratta inoltre di una funzione a 3 variabili e quindi fzero non dovrebbe essere utilizzabile. È corretto? Il testo del tuo problema non è chiarissimo.

Il testo dell'esercizio preso da un libro di calcolo numerico recita:
"Uno dei metodi per la risoluzione di equazioni differenziali ordinarie, y'(t)=f(t,y(t)), t$in$[$t_0$,T],y($t_0$)=$t_0$ è il metodo di eulero implicito. In cui $y_n$$~~$y($t_n$),$f_n$ =f($t_n$,y(t)), $t_n$=$t_0$+nh.
Utilizzare tale metodo nel caso in cui $t_0$=0, T=10, N=100, $y_0$=$(1,2)^T$ e, se y =$ (x_1,x_2)^T$, f(t,y) = $(-10^3x_1+sinx_1cosx_2, -2x_2+sinx_1cosx_2)^T$."

Ora provo ad implementare il codice, ti faccio risapere. Grazie!

[apatriarca]

L'idea di base non cambia, volevo solo essere sicuro che non avessi incasinato le variabili (la scelta dei nomi \(x_1\) e \(x_2\) per le componenti di \(y\) mi sembrava molto strana). Ovviamente non è più possibile usare fzero, ma puoi sempre risolvere il sistema non lineare corrispondente usando fsolve.

Dovrai quindi:
1. Calcolarti h con il metodo che ti ho già mostrato
2. Per ogni i che va da 1 a N calcoli un nuovo valore di y a partire da quello precedente utilizzando fsolve e l'equazione implicita che hai del metodo e lo memorizzi in un vettore (nota che non hai bisogno di aggiornare t perché non interviene nella derivata).
Appena ho tempo provo a scriverti una funzione per farlo. Nel frattempo prova a farlo tu e se ci sono cose non chiare chiedi pure.

[Gasso1]

function [y] = euleroim(t0,y0,tn,f)
%t0=0;y0=[1;2];tn=10;
%f=inline('[-10^3x+sinx*cosy;-2*y+sinx*cosy]')
%
T=10; N=100; h = (T - t0)/N;
n = fix((tn-t0)/h) + 1;
y(1,:)= y0;
for i=2:n
y(i,:)= y(i-1,:)+ h*fsolve(f,y(i-1,:));
end
end
Questo è il codice che ho implementato, nei commenti ci sono i valori rispettivamente di t0, y0, tn e f. Mi da un errore su fsolve. Forse il problema sta nel fatto che non uso "t", ma la funzione f(t,y) non dipende da t ma solo da y.Giusto?

[apatriarca]

Utilizza i tag code per inserire il codice in modo da renderlo maggiormente leggibile e scrivi sempre l'errore che ottieni eseguendo o compilando un determinato codice. Entrambe le cose aiutano tantissimo nel trovare gli errori. Credo che il problema nel tuo codice sia che la funzione che hai passato ad fsolve è una funzione a due variabili, ma ne hai passato solo uno. O almeno questo è quello che credo significasse l'errore che ho ottenuto.

Non ho mai scritto le funzioni inline in quel modo ma ho sempre usato la forma seguente (qui trovi una guida su come usare questa sintassi):

f = @(y) [-10^3 * y(1) + sin(y(1)) * cos(y(2)); -2 * y(2) + sin(y(1)) * cos(y(2))];

che si può forse anche riscrivere nella forma alternativa e più sintetica:

f = @(y) [-10^3; -2] .* y + sin(y(1))*cos(y(2));

Ovviamente puoi anche scrivere la stessa formula con il tuo metodo, ricordandoti però di scrivere x1 e y2 come elementi di un array perché dovranno essere passati in questo modo. Volendo però definire il metodo di eulero nel modo più generico possibile, dovrai allora definire la funzione in modo che f abbia due parametri, il primo che rappresenta il tempo ed è uno scalare mentre il secondo è un vettore colonna che rappresenta gli altri parametri. Dovrai quindi definire f come segue:

f = @(t, y)  [-10^3; -2] .* y + sin(y(1))*cos(y(2));

Nota come in questo caso sia molto più comoda questa sintassi.

fsolve serve per risolvere un sistema non-lineare. Il seguente codice

fsolve(f,y(i-1,:))

dovrebbe infatti trovare uno zero della funzione f, ma non è questo il nostro scopo. Il nostro scopo è quello di trovare uno zero della funzione \( y_{n} - y_{n+1} + h \, f(t, y_{n+1} \) in funzione di \( y_{n+1} \). Quello che devi quindi scrivere è qualcosa come il seguente:

y(i,:) = fsolve(@(x) y(i-1,:) - x + h*f(t0 + ih, x), y(i-1,:));

Un altro errore che ho notato nel tuo codice è l'uso di y(i,:) invece che y(:,i). I punti sono infatti vettori colonna e non riga. In ogni caso il seguente è la mia versione del codice (che però non ho avuto modo di provare):

function [y] = euleroim(f, t0, tn, N, y0)

h = (tn - t0)/N;
y = [y0 zeros(2, N)];

for i=2:(N+1)
    y(:,i) = fsolve(@(x) y(:,i-1) - x + h*f(t0 + i*h, x), y(:,i-1));
end

end

[Gasso1]

In base alle tue indicazioni ho "ricostruito" il codice:

function [y] = euleroim(t0,y0,tn,f)
N=100;h = (tn - t0)/N;
y = [y0 zeros(2, N)];
for i=2:(N+1)
y(:,i) = fsolve(@(x) y(:,i-1) - x + h*f(t0 + i*h, x), y(:,i-1));
end
end

Ho notato che scrivendo la funzione f nei modi seguenti:

f=inline('[-10^3y(1)+siny(1)*cosy(2);-2*y(2)+siny(1)*cosy(2)]')
f = @(y) [-10^3; -2] .* y + sin(y(1))*cos(y(2));
f = @(y) [-10^3 * y(1) + sin(y(1)) * cos(y(2)); -2 * y(2) + sin(y(1)) * cos(y(2))];

generava il seguente codice di errore:

??? Error using ==> inline.subsref at 17
Too many inputs to inline function.

Oppure
??? Error using ==> @(y)[-10^3;-2].*y+sin(y(1))*cos(y(2))
Too many input arguments.

Oppure
??? Error using ==> @(y)[-10^3*y(1)+sin(y(1))*cos(y(2));-2*y(2)+sin(y(1))*cos(y(2))]
Too many input arguments.

Seguito Da:
Error in ==> euleroim>@(x)y(:,i-1)-x+h*f(t0+i*h,x) at 11
y(:,i) = fsolve(@(x) y(:,i-1) - x + h*f(t0 + i*h, x), y(:,i-1));

Error in ==> fsolve at 254
            fuser = feval(funfcn{3},x,varargin{:});

Error in ==> euleroim at 11
y(:,i) = fsolve(@(x) y(:,i-1) - x + h*f(t0 + i*h, x), y(:,i-1));

Caused by:
    Failure in initial user-supplied objective function evaluation. FSOLVE cannot continue.

Mentre definendo f come:

f = @(t, y)  [-10^3; -2] .* y + sin(y(1))*cos(y(2));

il metodo viene eseguito, producendo in output y.

Mi compare questo messaggio:

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the default value of the function tolerance, and
the problem appears regular as measured by the gradient.

<stopping criteria details>

ma penso che sia normale.
Attendo conferma, grazie ancora per l'aiuto!

[apatriarca]

Usando Octave non mi ha mostrato quel messaggio. In ogni caso è un messaggio di successo per cui suppongo abbia restituito il risultato corretto.

[Gasso1]

Ok, allora grazie ancora