Elastic pendulum with runge kutta 4

Something is wrong with the code because the answer I get is not the correct answer, I checked the code a hundred times but I do not find what the problem is.

Edit: removed the h multiplication in this equations

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;

In this functions do i need to add t as the first varible like this?

f1=@(t,A2) A2;
f2=@(t,A1,A2,R1,R2) -2*(R2./R1).*A2-g*sin(A1);
f3=@(t.R2) R2;
f4=@(t,A1,A2,R1) R1.*(A2.*A2)+g*cos(A1)-k*(R1-l_0)/m;

What about the the half a step h/2 do I need to consider it in this problem?

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)),(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

[Jan]

We cannot reproduce the problem without knowing the input arguments you use. A comparison with a reference solution is possible only, if this solution is shown also.

[Netanel Malihi]

yes of course! sorry about that

(71,1000,8,10,pi/3,0,0,0.01,10)

[Torsten]

Six errors in your code:

k24=h*f4((A1(i)+h*k11/2),(A2(i)+h*k12/2),(R1(i)+h*k13/2));

k34=h*f4((A1(i)+h*k21/2),(A2(i)+h*k22/2),(R1(i)+h*k23/2));

[Netanel Malihi]

thnks ,but the answer is still not correct. im getting diffrent vectors using ode45 and runge kutta .

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

[Netanel Malihi]

I will edit it back asap, it was a mistake . My Sencier apology.

[Stephen23]

Original question by Netanel Malihi retrieved from Google Cache:

Elastic pendulum with runge kutta 4

Something is wrong with the code because the answer I get is not the correct answer, I checked the code a hundred times but I do not find what the problem is.

function [Theta,r]=RK4_1(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)+h*k12/2));
    k22=h*f2((A1(i)+h*k11/2),(A2(i)+h*k12/2),(R1(i)+h*k13/2),(R2(i)+h*k14/2));
    k23=h*f3((R2(i)+h*k14/2));
    k24=h*f4((A1(i)+h*k11),(A2(i)+h*k12),(R1(i)+h*k13));
    k31=h*f1((A2(i)+h*k22/2));
    k32=h*f2((A1(i)+h*k21/2),(A2(i)+h*k22/2),(R1(i)+h*k23/2),(R2(i)+h*k24/2));
    k33=h*f3((R2(i)+h*k24/2));
    k34=h*f4((A1(i)+h*k21),(A2(i)+h*k22),(R1(i)+h*k23));
    
    k41=h*f1((A2(i)+k32));
    k42=h*f2((A1(i)+h*k31),(A2(i)+h*k32),(R1(i)+h*k33),(R2(i)+h*k34));
    k43=h*f3((R2(i)+h*k34));
    k44=h*f4((A1(i)+h*k31),(A2(i)+h*k32),(R1(i)+h*k33));
    
    A1(i+1) = A1(i) + h*(k11+2*k21+2*k31+k41)/6;
    A2(i+1) = A2(i) + h*(k12+2*k22+2*k32+k42)/6;
    R1(i+1) = R1(i) + h*(k13+2*k23+2*k33+k43)/6;
    R2(i+1) = R2(i) + h*(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

Answer 1

[Sulaymon Eshkabilov]

Here are a few corrections made in your code:

function [Theta,r]=RK4_L(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)+h*k11/2));
    k22=h*f2((A1(i)+h*k12/2),(A2(i)+h*k12/2),(R1(i)+h*k12/2),(R2(i)+h*k12/2));
    k23=h*f3((R2(i)+h*k13/2));
    k24=h*f4((A1(i)+h*k14),(A2(i)+h*k14),(R1(i)+h*k14));

    k31=h*f1((A2(i)+h*k21/2));
    k32=h*f2((A1(i)+h*k22/2),(A2(i)+h*k22/2),(R1(i)+h*k22/2),(R2(i)+h*k22/2));
    k33=h*f3((R2(i)+h*k23/2));
    k34=h*f4((A1(i)+h*k24),(A2(i)+h*k24),(R1(i)+h*k24));
      
    k41=h*f1((A2(i)+k31));
    k42=h*f2((A1(i)+h*k32),(A2(i)+h*k32),(R1(i)+h*k32),(R2(i)+h*k32));
    k43=h*f3((R2(i)+h*k33));
    k44=h*f4((A1(i)+h*k34),(A2(i)+h*k34),(R1(i)+h*k34));
     
      A1(i+1) = A1(i) + h*(k11+2*k21+2*k31+k41)/6;
      A2(i+1) = A2(i) + h*(k21+2*k22+2*k32+k42)/6;
      R1(i+1) = R1(i) + h*(k13+2*k23+2*k33+k43)/6;
      R2(i+1) = R2(i) + h*(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]

thnks, but i believe it's not the correct way because the amplitude of the theta in the graph is growing