implementazione Runge-Kutta

[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.

Codice:

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

Codice:

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

dopo aver ovviamente definito prima

Codice:

f=@(t,y) -y

,

Codice:

y0=1

, $k$ ad esempio 0.01;

[cooper]

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 ts cambiandolo con ts=ceil((tf-t0)/k) 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