implementazione Runge-Kutta

Messaggioda cooper » 03/11/2018, 20:39

Ciao a tutti. dovrei implementare un codice per risolvere EDO con un RK esplicito di ordine qualunque.
Codice:
function [yout,tt] = RK2_esplicito(fun,t0,tf,y0,tab,h)
t=t0;
%tableau
[c,A,b]=tab();
y = y0(:);

%numero passi
s=length(b);
k = y*zeros(1,length(b)+1);
k(:,1) = feval(fun,t,y);
%variabile per contare i tempi per plottare
n=ceil((tf-t0)/h)+1;

%creo un vettore di tempi in cui valuto la funz
tt=zeros(1,n);
%la prima colonna è il tempo iniziale
tt(1)=t0;
%ciclo sul tempo
for i=2:(n-1)
    tt(i)=tt(i-1)+h;
end
%avevo ceil + 1, quindi aggiungo l'ultima colonna
tt(n)=tf;

while(t < tf)
if((t+h) > tf)
    h=tf-t;
end
%calcolo le k
k(:,1) = feval(fun,t,y);
for j=1:s
    k(:,j+1) = feval(fun, t+c(j)*h, y+h*f*A(:,j));
end

y = y + h*f*b(:,1);

t=t+h;

end
yout=y;
end


i tableau (tab nelle dipendenze della funzione) sono in un file e non so come fare per dirgli questa cosa. ho provato a passare 'nomefile.m' nelle dipendenze ma non esce. premetto che non ho mai seguito nessun corso di matlab e quindi anche i comandi più semplici mi sono oscuri.
non potendolo provare dato che non so come passare i tableau, vedete già qualche errore nel codice?
grazie in anticipo a tutti!
cooper
Cannot live without
Cannot live without
 
Messaggio: 2244 di 4642
Iscritto il: 25/07/2014, 09:19

Re: implementazione Runge-Kutta

Messaggioda feddy » 03/11/2018, 22:19

Ciao cooper, una cosa al volo: il tableau lo deve inserire l'utente come "matrice", esatto? Perché se così fosse a partire dalla routine che usi non ci dovrebbero essere grossi problemi.
Se è così e non va, allora leggo il codice.

MatLab è ottimizzato per evitare, quando possibile, cicli (googla BLAS MatLab). Per cui, per creare il vettore dei tempi basta solamente scrivere
Codice:
 t=linspace(t0,tf,n+1)
Questo genera un array da t0 a tf di n nodi equi-spaziati. Il "while (t<tf)" non credo serva in realtà... un ciclo for sugli $n$ forse è più chiaro. Quel while sarebbe più adatto nel caso in cui utilizzi un passo adattivo e non sai a priori la lunghezza di ogni passo.
Avatar utente
feddy
Moderatore
Moderatore
 
Messaggio: 2267 di 5934
Iscritto il: 26/06/2016, 00:25
Località: SISSA

Re: implementazione Runge-Kutta

Messaggioda feddy » 03/11/2018, 22:28

Altra cosa... vedo che quando calcoli $y$ nelle ultime righe c'è $f$ che è una variabile non dichiarata. Inoltre, nei parametri passati dall'utente $h$ è un parametro che non serve, tanto se il passo è costante è definito automaticamente come $\frac{tf-t_0}{n}$.
Avatar utente
feddy
Moderatore
Moderatore
 
Messaggio: 2268 di 5934
Iscritto il: 26/06/2016, 00:25
Località: SISSA

Re: implementazione Runge-Kutta

Messaggioda cooper » 04/11/2018, 00:41

ciao e grazie per la risposta.
feddy ha scritto:il tableau lo deve inserire l'utente come "matrice", esatto?

i tableau sono dati come matrici, inserite a mano nel file a parte di cui parlavo. in pratica la professoressa ci passa i file con i tableau già scritti e noi li dobbiamo dare in pasto al metodo evitando così di scrivere noi le matrici
feddy ha scritto:per creare il vettore dei tempi basta solamente scrivere

questa particolarità di matlab non la sapevo. per quanto riguarda invece quel codice: con $n$ intendi esattamente la $n$ che ho definito io nel codice?
feddy ha scritto:Quel while sarebbe più adatto nel caso in cui utilizzi un passo adattivo e non sai a priori la lunghezza di ogni passo.

in realtà è stata proprio la prof ha consigliarci il while, forse in vista proprio di quei metodi (che per altro cominciamo luneidì)
feddy ha scritto:vedo che quando calcoli y nelle ultime righe c'è f che è una variabile non dichiarata.

certo perchè sono fesso io. #-o dovrebbe essere una k
feddy ha scritto:nei parametri passati dall'utente h è un parametro che non serve, tanto se il passo è costante è definito automaticamente come tf−t0n.

in effetti l'ho sempre passato anche se inutile.
cooper
Cannot live without
Cannot live without
 
Messaggio: 2245 di 4642
Iscritto il: 25/07/2014, 09:19

Re: implementazione Runge-Kutta

Messaggioda feddy » 04/11/2018, 00:47

Sì, con $n$ intendo esattamente il numero di time-steps ! Sinceramente questa cosa di passare il tableau in un file a parte non l'avevo mai vista. Per questa cosa penserei alla fine, intanto controlla di aver scritto correttamente la routine con esempi semplici di cui conosci la soluzione analitica :)
Avatar utente
feddy
Moderatore
Moderatore
 
Messaggio: 2269 di 5934
Iscritto il: 26/06/2016, 00:25
Località: SISSA

Re: implementazione Runge-Kutta

Messaggioda cooper » 04/11/2018, 02:00

ho provato a farlo partire modificando in questo modo il codice ma non mi disegna il grafiico; ne deduco quindi che ci sia qualcosa di sbagliato

Codice:
function [yout,tt] = RK2_esplicito(fun,t0,tf,y0,h)
t=t0;
%tableau
%[c,A,b]=tab();
c=[0 1/2]';
b=[1 1]';
A=[0 0;
   1/2 0];
y = y0(:);

%numero passi
s=length(b);
k = y*zeros(1,length(b));
n=ceil((tf-t0)/h)+1;
tt=linspace(t0,tf,n+1);
while(t < tf)
if((t+h) > tf)
    h=tf-t;
end
%calcolo le k
k(:,1) = feval(fun,t,y);
for j=1:(s-1)
    k(:,j+1) = feval(fun, t+c(j)*h, y+h*k*A(:,j));
end

y = y + h*k*b(:,1);

t=t+h;

end
yout=y;
end


la soluzione viene anche, ma quando però poi faccio plot(tt,y) per disegnare mi vengono gli assi ma non la funzione. :?:
cooper
Cannot live without
Cannot live without
 
Messaggio: 2246 di 4642
Iscritto il: 25/07/2014, 09:19

Re: implementazione Runge-Kutta

Messaggioda feddy » 04/11/2018, 02:12

Beh sicuramente non potrà plottare nulla: $t$ e $y$ come li hai scritti tu sono solo i valori della soluzione all'ultimo timestep. Devi registrarli in un array per plottare
Avatar utente
feddy
Moderatore
Moderatore
 
Messaggio: 2270 di 5934
Iscritto il: 26/06/2016, 00:25
Località: SISSA

Re: implementazione Runge-Kutta

Messaggioda cooper » 04/11/2018, 17:02

ho provato a modificare in questo modo ma non tornano le dimensioni. il problema è che non riesco a capire come correggere
Codice:
function [yout,tt] = RK2_esplicito(fun,t0,tf,y0,h)
t=t0;
%tableau
%[c,A,b]=tab();
c=[0 1/2]';
b=[0 1]';
A=[0 0;
   1/2 0];
%lunghezza del problema
d=length(y0);
%var x numero tempi
n=ceil((tf-t0)/h)+1;
%creo vettore per la soluzione
y=zeros(d,n);

%la prima colonna della soluzione è il punto di partenza iniziale
y(:, 1)=y0;


%numero passi
s=length(b);
k = y*zeros(1,s);

%creo un vettore di tempi in cui valuto la funz
tt=zeros(1,n);
%la prima colonna è il tempo iniziale
tt(1)=t0;

%ciclo sul tempo
for i=2:(n-1)
    tt(i)=tt(i-1)+h;
end
%avevo ceil + 1, quindi aggiungo l'ultima colonna
tt(n)=tf;

i=2;

while(t < tf)
if((t+h) > tf)
    h=tf-t;
end
%calcolo le k
k(:,1) = feval(fun,t,y);
for j=1:(s-1)
    k(:,j+1) = feval(fun, t+c(j)*h, y+h*k*A(:,j));
end

t=t+h;

y(:,i) = y(:,i-1) + h*k*b(:,1);

%aggiorno il contatore
i=i+1;

end
yout=y;
end


già all'inizio ovviamente non riesce a calcolare k=y*zeros(1,s). se però cambio la dimensione di y ad una matrice $d xx 1$ da problemi quando calcola k(:,1) = feval(fun,t,y);.
come posso risolvere?
sempre grazie per l'aiuto :)
cooper
Cannot live without
Cannot live without
 
Messaggio: 2247 di 4642
Iscritto il: 25/07/2014, 09:19

Re: implementazione Runge-Kutta

Messaggioda feddy » 04/11/2018, 17:13

Mi spiace ma in questo momento non riesco a controllare per bene il codice... prova a farti un esempio su un foglio a mano, dove magari la matrice è quella di un semplcie RK a due stadi, e poi cerca di capire dove stanno i problemi. Dovresti riuscirci.

Riguardo ai due pezzi che hai citato: il primo non ne capisco il senso, $K$ cambia ad ogni time-step, ti basta definirlo come una matrice di length(y0) righe e #stadi colonne, e sovrascriverla ogni volta. Riguardo alla seconda cosa che hai citato, magari da problemi perché "feval(...)" ritorna una matrice/vettore le cui dimensioni non coinicidono con l'assegnamento.

Della serie... "non sono bello ma debuggo"
Avatar utente
feddy
Moderatore
Moderatore
 
Messaggio: 2271 di 5934
Iscritto il: 26/06/2016, 00:25
Località: SISSA

Re: implementazione Runge-Kutta

Messaggioda cooper » 04/11/2018, 18:25

ritorno con una nuova versione ed un nuovo problema: grazie per la pazienza :oops:
adesso un risultato lo produce, non so se corretto ma lo produce. il problema è che non riesco nuovamente a fare questo benedetto disegno perchè il vettore dei tempi e quello della soluzione non hanno la stessa dimensione: quello dei tempi ha una colonna in meno! il codice è il seguente, guardalo ovviamente quando hai tempo! :-D
Codice:
function [yout,tt] = RK2_esplicito(fun,t0,tf,y0,h)
t=t0;
%tableau
%[c,A,b]=tab();
c=[0 1]';
b=[1/2 1/2]';
A=[0 0;
   1 0];
%lunghezza del problema
d=length(y0);
%variabile per contare i tempi per plottare
n=ceil((tf-t0)/h)+1;
%creo vettore per la soluzione
y=zeros(d,n);
%la prima colonna della solzione è il punto di partenza iniziale
y(:, 1)=y0;


%numero passi
s=length(b);
k = zeros(d,s);
%k(:,1) = feval(fun,t,y);

%creo un vettore di tempi in cui valuto la funz
tt=zeros(1,n);
%la prima colonna è il tempo iniziale
tt(1)=t0;


%ciclo sul tempo
for i=2:(n-1)
    tt(i)=tt(i-1)+h;
end
%avevo ceil + 1, quindi aggiungo l'ultima colonna
tt(n)=tf;
%

i=2;

while(t < tf)
if((t+h) > tf)
    h=tf-t;
end
%calcolo le k
k(:,1) = feval(fun,t,y(:,1));
for j=1:(s-1)
    k(:,j+1) = feval(fun, t+c(j)*h, y(:,j)+h*k*A(:,j));
end

t=t+h;

y(:,i) = y(:,i-1) + h*k*b(:,1);

%aggiorno il contatore
i=i+1;

end
yout=y;
end
cooper
Cannot live without
Cannot live without
 
Messaggio: 2250 di 4642
Iscritto il: 25/07/2014, 09:19

Prossimo

Torna a Analisi Numerica e Ricerca Operativa

Chi c’è in linea

Visitano il forum: Nessuno e 1 ospite