如何在此 matlab 中用 runge-kutta 替换 ode45 方法？

Question

我尝试了所有方法并到处寻找，但找不到解决我问题的方法。

clc
clear all 


%% Solving the Ordinary Differential Equation 
G = 6.67408e-11; %Gravitational constant 
M = 10; %Mass of the fixed object 
r = 1; %Distance between the objects 

tspan = [0 100000]; %Time Progression from 0 to 100000s 
conditions = [1;0]; %y0= 1m apart, v0=0 m/s

F=@(t,y)var_r(y,G,M,r);

[t,y]=ode45(F,tspan,conditions); %ODE solver algorithm

%%part1: Plotting the Graph 
% plot(t,y(:,1)); %Plotting the Graph 
% xlabel('time (s)') 
% ylabel('distance (m)')

%% part2: Animation of Results 
plot(0,0,'b.','MarkerSize', 40); 
hold on    %to keep the first graph 
for i=1:length(t) 
k = plot(y(i,1),0,'r.','MarkerSize', 12); 
pause(0.05); 
axis([-1 2 -2 2]) %Defining the Axis 
xlabel('X-axis') %X-Axis Label 
ylabel('Y-axis') %Y-Axis Label 
delete(k)
end 

function yd=var_r(y,G,M,r) %function of variable r 
g = (G*M)/(r + y(1))^2; 
yd = [y(2); -g]; 
end

这是我尝试用 runge kutta 方法替换 ode45 的代码，但它给了我错误。我的 runge kutta 函数：

function y = Runge_Kutta(f,x0,xf,y0,h)

n= (xf-x0)/h;
y=zeros(n+1,1);
x=(x0:h:xf);
y(1) = y0;


for i=1:n
  k1 = f(x(i),y(i));
  k2= f(x(i)+ h/2 , y(i) +h*(k1)/2);
  y(i+1) = y(i)+(h*k2);
end

plot(x,y,'-.M')
legend('RKM')
title ('solution of y(x)');
xlabel('x');
ylabel('y(x)')
hold on
end

Answer 1

y(1) 是 y 数据结构中的第一个标量元素（按 column-first 顺序计数）。您希望在 y 中生成一个列向量列表，因为您的 ODE 是一个状态维数为 2 的系统。因此您需要生成具有该格式的 y，y=zeros(length(x0),n+1); 然后解决将条目列为矩阵列 y(:,1)=x0 并在您提取或分配列表条目的每个地方进行相同的修改。

Matlab 介绍了各种 short-cuts，如果因此使用，会导致矛盾（我认为 script-hater rant（德语）在很大程度上仍然有效）。本质上，与其他系统不同，Matlab 可以直接访问矩阵的底层数据结构。 y(k) 是底层平面数组的元素（在 Matlab 中被解释为 column-first，就像在 Fortran 中一样，不像 Numpy，它是 row-first）。

只有 two-index 才能访问矩阵及其维度。所以 y(:,k) 是 k-th 矩阵列，y(k,:) 是 k-th 矩阵行。 single-index 访问适用于行或列向量，但在列表中收集此类向量时会立即导致问题，因为这些列表自动为矩阵。

Answer 2

在将您的 ode45( ) 解决方案转换为手动编写的 RK 方案之前，您的 ode45( ) 解决方案看起来甚至都不正确。看起来你设置了一个初始速度为 0 的引力问题，所以一个小物体只会落入一条直线上的大质量 M（直线运动），这就是你有标量位置和速度的原因。

按照这个假设，r 是您应该即时计算的东西，而不是用作导数函数的固定输入。例如，我会期待这样的事情：

F=@(t,y)var_r(y,G,M); % get rid of r
    :

function yd=var_r(y,G,M) % function of current position y(1) and velocity y(2) 
g = (G*M)/y(1)^2; % gravity accel based on current position
yd = [y(2); -g]; % assumes y(1) is positive, so acceleration is negative
end

小对象必须以正初始位置开始，这样衍生代码才有效，正如您编写的那样。当小物体落入大质量 M 中时，上述情况只会持续到它撞击 M 的表面或大气层。或者如果您将 M 建模为点质量，那么该方案将变得越来越难以正确集成，因为加速度变得大而无界，因为小质量非常接近点质量 M。在这种情况下，您肯定需要可变步长方法。如果解决方案“通过”质量 M，则解决方案将变得无效。事实上，一旦速度变得太大，由于相对论效应，整个设置将变得无效。

也许你可以更详细地解释一下你的系统是否应该这样设置，以及集成的目的是什么。如果真的应该是2D或者3D的问题，那就需要增加更多的状态了。

对于您的手动 Runge-Kutta 代码，您完全忘记了对速度进行积分，因此这将失败得很惨。您需要一步一步地携带一个 2 元素状态，而不是像您目前所做的那样携带一个标量。例如，像这样：

y=zeros(2,n+1); % 2-element state as columns of the y variable
x=(x0:h:xf);
y(:,1) = y0; % initial state is the first 2-element column
% change all the scalar y(i) to column y(:,i)
for i=1:n
  k1 = f(x(i),y(:,i));
  k2= f(x(i)+ h/2 , y(:,i) +h*(k1)/2);
  y(:,i+1) = y(:,i)+(h*k2);
end

plot(x,y(1,:),'-.M') % plot the position part of the solution

这一切都假设传入的 f 与原始代码中的 F 相同。

如何在此 matlab 中用 runge-kutta 替换 ode45 方法？

how to replace the ode45 method with the runge-kutta in this matlab?

matlab

runge-kutta

ode45