我试了所有的方法，到处找，但都找不到解决我问题的方法。
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
在将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时，加速度会变得没有边界。在这种情况下，你肯定需要一个变步长方法。如果解决方案"通过"，那么它就无效了。事实上，一旦速度太大，由于相对论效应，整个设置就失效了。
也许您可以更详细地解释您的系统是否应该以这种方式设置，以及集成的目的是什么。如果它真的应该是一个2D或3D问题，那么就需要添加更多的状态。
对于你的手工龙格-库塔代码，你完全忘记了积分速度，所以这将会惨败。您需要一步一步地携带一个双元素状态，而不是像您目前所做的那样携带一个标量。例如，像这样:
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相同。
y(1)是y数据结构中的第一个标量元素(按列优先顺序计数)。您希望在y中生成列向量列表，因为您的ODE是具有状态维数2的系统。因此，您需要用这种格式生成y,y=zeros(length(x0),n+1);，然后将列表项定位为矩阵列y(:,1)=x0，并在提取或分配列表项的每个地方进行相同的修改。
Matlab引入了各种快捷方式，如果使用这些快捷方式，则会导致矛盾(我认为脚本憎恨者的咆哮(德语)在很大程度上仍然有效)。本质上，与其他系统不同，Matlab提供了对矩阵底层数据结构的直接访问。y(k)是底层平面数组的元素(它在Matlab中像在Fortran中一样被解释为列优先，而不像Numpy，它是行优先)。
只有双索引访问矩阵的维度。所以y(:,k)是矩阵的第k列y(k,:)是矩阵的第k行。单索引访问对于行或列向量很好，但是当在列表中收集这些向量时，会立即导致问题，因为这些列表自动是矩阵。