Elastic pendulum with runge kutta 4

3 ビュー (過去 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 件のコメント
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
Sulaymon Eshkabilov 2022 年 1 月 1 日
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 ExchangeProgramming についてさらに検索

タグ

製品


リリース

R2019b

Community Treasure Hunt

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

Start Hunting!

Translated by