hi I am trying to use this code to solve RK4 equation and I keep receiving this error . what should I do?

조회 수: 1 (최근 30일)
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);

채택된 답변

James Tursa
James Tursa 2020년 12월 24일
Your implementation has several errors:
  • 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
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
ela mti
ela mti 2020년 12월 25일
actually A should be reduced over the time I tried to fix the code with your suggestion code but it cant be right

댓글을 달려면 로그인하십시오.

추가 답변 (0개)

카테고리

Help CenterFile Exchange에서 Matrices and Arrays에 대해 자세히 알아보기

태그

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by