Elastic pendulum with runge kutta 4

4 views (last 30 days)
Netanel Malihi
Netanel Malihi on 1 Jan 2022
Edited: Netanel Malihi on 3 Jan 2022
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
  7 Comments
Stephen23
Stephen23 on 3 Jan 2022
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

Sign in to comment.

Answers (1)

Sulaymon Eshkabilov
Sulaymon Eshkabilov on 1 Jan 2022
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
  1 Comment
Netanel Malihi
Netanel Malihi on 2 Jan 2022
thnks, but i believe it's not the correct way because the amplitude of the theta in the graph is growing

Sign in to comment.

Tags

Products


Release

R2019b

Community Treasure Hunt

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

Start Hunting!