MATLAB Answers

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

조회 수: 18(최근 30일)
ela mti
ela mti 24 Dec 2020
댓글: James Tursa 25 Dec 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;
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;
tout(j) = t;
Aout(j) = A;
Index exceeds the number of array elements (1).
Error in ivp_RK4 (line 7)
k1 = A(i+1)-exp(-2*t);

  댓글 수: 0

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

채택된 답변

James Tursa
James Tursa 24 Dec 2020
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;
After the loop finishes, the answers are in the Aout and tout vectors. So you would
I am confused what the i_out is supposed to be used for. Can you clarify?

  댓글 수: 3

ela mti
ela mti 25 Dec 2020
first of all thank you so much !
i_out is frequency to output x and y to [xout, yout]
ela mti
ela mti 25 Dec 2020
actually A should be reduced over the time I tried to fix the code with your suggestion code but it cant be right

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

추가 답변(0개)


Community Treasure Hunt

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

Start Hunting!

Translated by