Elastic pendulum with runge kutta 4

1 次查看（过去 30 天）
Netanel Malihi
Netanel Malihi 2022-1-1
编辑： Netanel Malihi 2022-1-3
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 个评论
显示 5更早的评论
Netanel Malihi
Netanel Malihi 2022-1-3
I will edit it back asap, it was a mistake . My Sencier apology.
Stephen23
Stephen23 2022-1-3
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

请先登录，再进行评论。

请先登录，再回答此问题。

回答（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
  1 个评论
Netanel Malihi
Netanel Malihi 2022-1-2
thnks, but i believe it's not the correct way because the amplitude of the theta in the graph is growing

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

产品


版本

R2019b

Community Treasure Hunt

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

Start Hunting!

Translated by