Runge kutta 4 probelm

1 次查看（过去 30 天）
Netanel Malihi
Netanel Malihi 2022-1-2
评论： James Tursa 2022-1-31
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
  3 个评论
显示 1更早的评论
Rena Berman
Rena Berman 2022-1-25

(Answers Dev) Restored edit

James Tursa
James Tursa 2022-1-31
k42 and k44 do not have the correct A1 terms.

请先登录，再进行评论。

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

回答（2 个）

Walter Roberson
Walter Roberson 2022-1-2
Ran in:
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
  1 个评论
Netanel Malihi
Netanel Malihi 2022-1-2
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.

请先登录，再进行评论。

Torsten
Torsten 2022-1-2
编辑：Torsten 2022-1-2
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

类别

在 Help Center 和 File Exchange 中查找有关 ARM Cortex-A Processors 的更多信息

标签

产品


版本

R2019b

Community Treasure Hunt

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

Start Hunting!

Translated by