- Looks like you are creating i as a vector and then using that for indexing in A(i)
- You are not using the k's properly
- You are not calculating k4 correctly (it should be a full step not a half step)
- You are using A as both a vector and a scalar when updating
- You are not saving the intermediate results in Aout and tout inside the loop
hi I am trying to use this code to solve RK4 equation and I keep receiving this error . what should I do?
1 visualización (últimos 30 días)
Mostrar comentarios más antiguos
ela mti
el 24 de Dic. de 2020
Comentada: James Tursa
el 25 de Dic. de 2020
function [tout,Aout] = ivp_RK4(t0,A0,h,n_out,i_out)
tout = zeros(n_out+1,1); tout(1) = 0;
Aout = zeros(n_out+1,1); Aout(1) = 1;
t = t0; A = A0;
for j=2:n_out+1;
i=1:i_out;
k1 = A(i)-exp(-2*t);
k2 =(A(i)+h/2)*k1-exp(-2*(t+h/2)) ;
k3 =(A(i)+h/2)*k2-exp(-2*(t+h/2));
k4 = (A(i)+h/2)*k3-exp(-2*(t+h/2));
A = A + (h/6)*(k1 + 2*k2 + 2*k3 + k4);
t = t + h;
end
tout(j) = t;
Aout(j) = A;
plot(t,A)
end
Index exceeds the number of array elements (1).
Error in ivp_RK4 (line 7)
k1 = A(i+1)-exp(-2*t);
0 comentarios
Respuesta aceptada
James Tursa
el 24 de Dic. de 2020
Your implementation has several errors:
To make things simpler and easier to read and debug, I would suggest making a function handle for the derivative. Then construct your code using this function handle, indexing both Aout and tout inside the loop. E.g.,
f = @(t,A) A-exp(-2*t);
tout = zeros(n_out+1,1); tout(1) = t0;
Aout = zeros(n_out+1,1); Aout(1) = A0;
for i=1:n_out
k1 = f(tout(i) , Aout(i) );
k2 = f(tout(i)+h/2, Aout(i)+k1*h/2);
k3 = f(tout(i)+h/2, Aout(i)+k2*h/2);
k4 = f(tout(i)+h , Aout(i)+k3*h );
Aout(i+1) = Aout(i) + (h/6)*(k1 + 2*k2 + 2*k3 + k4);
tout(i+1) = tout(i) + h;
end
After the loop finishes, the answers are in the Aout and tout vectors. So you would
plot(tout,Aout)
I am confused what the i_out is supposed to be used for. Can you clarify?
3 comentarios
Más respuestas (0)
Ver también
Categorías
Más información sobre Loops and Conditional Statements en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!