Implementazione Runge-Kutta

[cooper1]

Ciao a tutti. dovrei implementare un codice per risolvere EDO con un RK esplicito di ordine qualunque.

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!

[feddy]

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

 t=linspace(t0,tf,n+1) 

Questo genera un array da t0 a tf di n nodi equi-spaziati. Il "while (t

Cita

Segnala

[feddy]

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}$.

[cooper1]

ciao e grazie per la risposta.

"feddy":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":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":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":vedo che quando calcoli y nelle ultime righe c'è f che è una variabile non dichiarata.

certo perchè sono fesso io. dovrebbe essere una k

"feddy":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.

[feddy]

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

[cooper1]

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

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 [inline]plot(tt,y)[/inline] per disegnare mi vengono gli assi ma non la funzione.

[feddy]

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

[cooper1]

ho provato a modificare in questo modo ma non tornano le dimensioni. il problema è che non riesco a capire come correggere

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 [inline]k=y*zeros(1,s)[/inline]. se però cambio la dimensione di y ad una matrice $d xx 1$ da problemi quando calcola [inline]k(:,1) = feval(fun,t,y);[/inline].
come posso risolvere?
sempre grazie per l'aiuto

[feddy]

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"

[cooper1]

ritorno con una nuova versione ed un nuovo problema: grazie per la pazienza
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!

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

[feddy]

L'ho fatto partire sul problema $y'=-y, y(0)=1$, e produce una soluzione evidentemente errata.
[Regola generale in questo tipo di problemi: testalo su un problema di cui conosci la soluzione analitica, controlla se differiscono e cerca l'errore].

Vedo che il tuo tableau corrisponde al metodo di Heun... vedo di scriverti il codice per questo caso particolare, poi puoi adattarlo sicuramente senza troppa fatica.

function [t,y]=Heun(f,t0,tf,y0,k)


%Metodo di Heun- RK esplicito
%f=funzione del problema di Cauchy
%[t0 tf] vettore tspan
%y0 dato iniziale
%k passo (costante)

c=[0 1];
b=[1/2 1/2]; %tableau corrispondente
A=[0 0;
   1 0];
ts=(tf-t0)/k; %timesteps
y=NaN(length(y0),ts+1);
t=linspace(t0,tf,ts+1);

y(:,1)=y0;
K=NaN(length(y0),2); %Heun metodo a due stadi !
for n=1:ts
	K(:,1)=f(t(n)+k*c(1),y(:,n)+k);
	K(:,2)=f(t(n)+k*c(2),y(:,n)+k*A(2,1)*K(:,1));
	y(:,n+1)=y(:,n)+k*(b(1)*K(:,1)+b(2)*K(:,2));	
endfor

plot(t,y,'-o')
	
end

Ora da terminale (o command window su MatLab) prova a eseguire

 [t,y]=Heun(f,0,1,y0,k)

dopo aver ovviamente definito prima

f=@(t,y) -y

,

y0=1

, $k$ ad esempio 0.01;

[cooper1]

ho provato a far partire il tuo codice ma mi dava l'errore seguente:

Array indices must be positive integers or logical values.

Error in Heun (line 21)
K(:,1)=f(t(n)+k*c(1),y(:,n)+k);

ho allora provato a modificare il calcolo di [inline]ts[/inline] cambiandolo con [inline]ts=ceil((tf-t0)/k)[/inline] per assicurarmi fosse un intero, ma l'errore è rimasto.

[feddy]

Mmm strano a me partiva senza problemi, sicuro di aver passato i parametri correttamente? Io uso Octave, ma la sintassi è sicuramente compatibile in questo caso