为什么我的 Runge kutta 代码不能正常工作？答案

【问题标题】：Why my Runge kutta code is not working correctly?为什么我的 Runge kutta 代码不能正常工作？
【发布时间】：2022-01-02 22:15:34
【问题描述】：

下面是我的 4 阶 Runge-Kutta 代码，但我一直得到错误的答案，我看不出代码有什么问题。我希望有人能帮助我指出正确的方向来修复我的代码。

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

【问题讨论】：

我的 RK4 在这里：github.com/sjhalayka/mercury_orbit_glut/blob/…
你能抛出你的数学问题陈述吗？每个函数的初始值？当您递归查找值时，我看到了问题 - 即您的函数是 f1 f2 f3 f4 但您找到的是 A1 A2 R1 R2。尝试 CTRL+C -> 来自 wiki en.wikipedia.org/wiki/Runge–Kutta_methods 的 CTRL+V 公式
是的，我称 A1=theta 所以 A2=dtheta, R1=r, R2= dr。我的初始条件：dr/dt=0，dtheta/dt=0，theta=pi/6，r=10。输入：(71,1000,8,10,pi/6,0,0,1,10)

标签： matlab

【解决方案1】：

您的 k42 和 k44 计算中似乎缺少代码。你有这个：

k42=h*f2((A1(i)),(A2(i)+k32),(R1(i)+k33),(R2(i)+k34));
:
k44=h*f4((A1(i)),(A2(i)+k32),(R1(i)+k33));

但看起来你需要这个来代替 A1 部分：

k42=h*f2((A1(i)+k31),(A2(i)+k32),(R1(i)+k33),(R2(i)+k34));
:
k44=h*f4((A1(i)+k31),(A2(i)+k32),(R1(i)+k33));

如果您发布您正在求解的微分方程，这样我们可以检查 f1、f2、f3 和 f4 函数句柄，这将有所帮助。

看起来您正在为步长传递一个 h=1 值。对于看起来像某种钟摆问题的情况，这方式太大了。你应该把它降到 0.01 左右。

最后，如果您对微分方程代码进行矢量化处理，将来会有所帮助。这意味着制作一个函数句柄，它接受一个 4 元素状态向量并输出一个 4 元素导数向量。这样，您的 RK4 代码只有一组矢量化方程，而不是像您目前拥有的四组单独的标量方程。由于您可以使用相同的派生函数，因此还可以轻松地根据内置 ode45( ) 函数检查您的代码。

【讨论】：