Runge kutta 4 probelm

What's the problem with this runge kutta code ?

for this parametrs (71,1000,8,10,pi/6,0,0,0.1,10)

function [Theta,r]=test3rk4(m,k,l_0,R1,A1,R2,A2,h,t_end)
t=0:h:t_end;
% A1 = zeros(1,numel(t)) % A2 = A1; R1 = A1; R2 = A1;
% R1 = zeros(1,numel(t))
g=9.81;
f1=@(A2) A2;
f2=@(A1,A2,R1,R2) -2*(R2./R1).*A2-g*sin(A1);
f3=@(R2) R2;
f4=@(A1,A2,R1) R1.*(A2.*A2)+g*cos(A1)-k*(R1-l_0)/m;
for i=1:(length(t)-1)
    k11=h*f1(A2(i));
    k12=h*f2(A1(i),A2(i),R1(i),R2(i));
    k13=h*f3(R2(i));
    k14=h*f4(A1(i),A2(i),R1(i));
     
    k21=h*f1((A2(i)+k12/2));
    k22=h*f2((A1(i)+k11/2),(A2(i)+k12/2),(R1(i)+k13/2),(R2(i)+k14/2));
    k23=h*f3((R2(i)+k14/2));
    k24=h*f4((A1(i)+k11/2),(A2(i)+k12/2),(R1(i)+k13/2));
    
    k31=h*f1((A2(i)+k22/2));
    k32=h*f2((A1(i)+k21/2),(A2(i)+k22/2),(R1(i)+k23/2),(R2(i)+k24/2));
    k33=h*f3((R2(i)+k24/2));
    k34=h*f4((A1(i)+k21/2),(A2(i)+k22/2),(R1(i)+k23/2));
      
    k41=h*f1((A2(i)+k32));
    k42=h*f2((A1(i+1)),(A2(i)+k32),(R1(i)+k33),(R2(i)+k34));
    k43=h*f3((R2(i)+k34));
    k44=h*f4((A1(i)),(A2(i)+k32),(R1(i)+k33));
     
      A1(i+1) = A1(i) + (k11+2*k21+2*k31+k41)/6;
      A2(i+1) = A2(i) + (k12+2*k22+2*k32+k42)/6;
      R1(i+1) = R1(i) + (k13+2*k23+2*k33+k43)/6;
      R2(i+1) = R2(i) + (k14+2*k24+2*k34+k44)/6;
end
Theta=A1;
r=R1;
subplot(2,2,1)
plot(t,A1,'LineWidth',1)
title('theta vs t')
xlabel('t','fontsize',14,'fontweight','bold')
ylabel('Theta','fontsize',14,'fontweight','bold')
subplot(2,2,2)
plot(t,R1,'LineWidth',1)
xlabel('t','fontsize',14,'fontweight','bold')
ylabel('r','fontsize',14,'fontweight','bold')
subplot(2,2,[3,4])
 x=R1.*sin(A1)
 y=-R1.*cos(A1)
 plot(x,y)
end

[John D'Errico]

When you edit away your question, you are being rude to the person who spent their time answering your question. You make it impossible for others to get help from that same problem. Someone else might be able to learn something. Now your question is useless to others.

If your plan was to do exactly this, then you should never ask the question in the first place.

[Rena Berman]

(Answers Dev) Restored edit

[James Tursa]

k42 and k44 do not have the correct A1 terms.

Answer 1

[Walter Roberson]

k42=h*f2((A1(i+1)),(A2(i)+k32),(R1(i)+k33),(R2(i)+k34));

That line requires that A1(i+1) exists before A1(i+1) has been assigned to.

[t,r] = test3rk4(71,1000,8,10,pi/6,0,0,0.1,10)

Index exceeds the number of array elements. Index must not exceed 1.

Error in solution>test3rk4 (line 29)
k42=h*f2((A1(i+1)),(A2(i)+k32),(R1(i)+k33),(R2(i)+k34));

function [Theta,r]=test3rk4(m,k,l_0,R1,A1,R2,A2,h,t_end)
t=0:h:t_end;
% A1 = zeros(1,numel(t)) % A2 = A1; R1 = A1; R2 = A1;
% R1 = zeros(1,numel(t))
g=9.81;
f1=@(A2) A2;
f2=@(A1,A2,R1,R2) -2*(R2./R1).*A2-g*sin(A1);
f3=@(R2) R2;
f4=@(A1,A2,R1) R1.*(A2.*A2)+g*cos(A1)-k*(R1-l_0)/m;
for i=1:(length(t)-1)
    k11=h*f1(A2(i));
    k12=h*f2(A1(i),A2(i),R1(i),R2(i));
    k13=h*f3(R2(i));
    k14=h*f4(A1(i),A2(i),R1(i));
     
    k21=h*f1((A2(i)+k12/2));
    k22=h*f2((A1(i)+k11/2),(A2(i)+k12/2),(R1(i)+k13/2),(R2(i)+k14/2));
    k23=h*f3((R2(i)+k14/2));
    k24=h*f4((A1(i)+k11/2),(A2(i)+k12/2),(R1(i)+k13/2));
    
    k31=h*f1((A2(i)+k22/2));
    k32=h*f2((A1(i)+k21/2),(A2(i)+k22/2),(R1(i)+k23/2),(R2(i)+k24/2));
    k33=h*f3((R2(i)+k24/2));
    k34=h*f4((A1(i)+k21/2),(A2(i)+k22/2),(R1(i)+k23/2));
      
    k41=h*f1((A2(i)+k32));
    k42=h*f2((A1(i+1)),(A2(i)+k32),(R1(i)+k33),(R2(i)+k34));
    k43=h*f3((R2(i)+k34));
    k44=h*f4((A1(i)),(A2(i)+k32),(R1(i)+k33));
     
      A1(i+1) = A1(i) + (k11+2*k21+2*k31+k41)/6;
      A2(i+1) = A2(i) + (k12+2*k22+2*k32+k42)/6;
      R1(i+1) = R1(i) + (k13+2*k23+2*k33+k43)/6;
      R2(i+1) = R2(i) + (k14+2*k24+2*k34+k44)/6;
end
Theta=A1;
r=R1;
subplot(2,2,1)
plot(t,A1,'LineWidth',1)
title('theta vs t')
xlabel('t','fontsize',14,'fontweight','bold')
ylabel('Theta','fontsize',14,'fontweight','bold')
subplot(2,2,2)
plot(t,R1,'LineWidth',1)
xlabel('t','fontsize',14,'fontweight','bold')
ylabel('r','fontsize',14,'fontweight','bold')
subplot(2,2,[3,4])
 x=R1.*sin(A1)
 y=-R1.*cos(A1)
 plot(x,y)
end

[Netanel Malihi]

yea it was my mistake i left the code with that mistake but still it's not working.well it's working but the answere in not correct.

Answer 2

[Torsten]

function main
  g = 9.81;
  m = 71;
  k = 1000;
  l_0 = 8;
  tstart = 0.0;
  tend = 10.0;
  h = 0.01;
  T = (tstart:h:tend).';
  Y0 = [pi/6, 0, 10, 0]; 
  f = @(t,y) [y(2), -2*y(4)/y(3)*y(2)-g*sin(y(1)), ...
              y(4), y(3)*y(2)^2+g*cos(y(1))-k*(y(3)-l_0)/m];
  Y = runge_kutta_RK4(f,T,Y0);
  plot(T,Y)                           
end

function Y = runge_kutta_RK4(f,T,Y0)
  N = numel(T);
  n = numel(Y0);
  Y = zeros(N,n);
  Y(1,:) = Y0;
  for i = 2:N
    t = T(i-1);
    y = Y(i-1,:);
    h  = T(i) - T(i-1);
    k0 = f(t,y);
    k1 = f(t+0.5*h,y+k0*0.5*h);
    k2 = f(t+0.5*h,y+k1*0.5*h);
    k3 = f(t+h,y+k2*h);
    Y(i,:) = y + h/6*(k0+2*k1+2*k2+k3);
  end
end