Runge kutta 4 probelm

Question

Netanel Malihi 2022-1-2

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1621060-runge-kutta-4-probelm

评论： 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更早的评论隐藏 1更早的评论

Rena Berman 2022-1-25

(Answers Dev) Restored edit

James Tursa 2022-1-31

k42 and k44 do not have the correct A1 terms.

请先登录，再进行评论。

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

Answer 1

Walter Roberson 2022-1-2

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1621060-runge-kutta-4-probelm#answer_866775

在 MATLAB Online 中打开

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 个评论
显示 -1更早的评论隐藏 -1更早的评论

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.

请先登录，再进行评论。

Answer 2

Torsten 2022-1-2

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1621060-runge-kutta-4-probelm#answer_866800

编辑：Torsten 2022-1-2

在 MATLAB Online 中打开

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 