运行 matlab 代码表明没有提供足够的输入，并且我们的 RHS 函数有问题

Question

此代码的目标是用数值模拟三个人的运动三个维度的身体。所有的物体都有恒定的质量。我们将确定物体位置随时间的变化。然而，当我们运行代码时，RHS 函数由于“没有提供足够的输入”而出现故障，我们无法弄清楚原因。

    % initial conditions and parameters
    clc; 
    clear all; 

    M1=1;   %mass of body 1
    M2=2;   %mass of body 2 
    M3=0.5; %mass of body 3 


    % the position of the three bodies in the initial conditions 

    P1_0=[0,0,0];
    P2_0=[1,0,0];
    P3_0=[0,0,1];

    P1=P1_0;
    P2=P2_0;
    P3=P3_0;

    % the initial speed of the three bodies

    a=0.05;
    V1_0=[a,a,a];
    V2_0=[-a,a,a];
    V3_0=[-a,-a,-a];

    V1=V1_0;
    V2=V2_0;
    V3=V3_0;

    X0=[P1_0,V1_0,P2_0,V2_0,P3_0,V3_0]'; % the initial positions and speeds of the three bodies in 
    one column vector.

    [t, X]= ode45(@(t,X) RHS(t,X), [0 3], X0);


    h=0.1;
    K1=(h*RHS(X0));
    K2=h*(RHS(X0+0.5*K1));
    K3=h*(RHS(X0-K1+(2*K2)));
    X1=X0+(1/6)*(K1+4*K2+K3);

    function Y=RHS(t,X)
    V1=X(4:6);
    V2=X(10:12);
    V3=X(16:18);
    P1=X(1:3);
    P2=X(7:9);
    P3=X(13:15);
    M1=1;
    M2=2;
    M3=0.5;
    Y=[V1;Force(P1,M1,P2,M2,P3,M3)/M1;V2;Force(P2,M2,P1,M1,P3,M3)/M2;V3;Force(P3,M3,P1,M1,P2,M2)/M3];
    end 


    function F=Force(P,M,Pa,Ma,Pb,Mb)

    F=(-M*Ma*(P-Pa))/(norm(P-Pa)^3)-(M*Mb*(P-Pb))/(norm(P-Pb)^3);
    end

Answer 1

托马斯：

我只是运行你的代码，它可以完美地工作，包括对 ode45 的标注。但是，我确实遇到了您在问题中声称的相同错误。该错误由行 K1=(h*RHS(X0)) 触发，我认为这是您尝试设置 Runge-Kutta 方法的块的开头。你得到错误的原因是在该行中你只给你的 RHS 函数一个单一的输入，然而，你已经定义 RHS 接受两个输入（t，和 X):

    function Y=RHS(t,X)
    V1=X(4:6);
    V2=X(10:12);
    V3=X(16:18);
    P1=X(1:3);
    P2=X(7:9);
    P3=X(13:15);
    M1=1;
    M2=2;
    M3=0.5;
    Y=[V1;Force(P1,M1,P2,M2,P3,M3)/M1;V2;Force(P2,M2,P1,M1,P3,M3)/M2;V3;Force(P3,M3,P1,M1,P2,M2)/M3];
    end

由于RHS函数根本没有使用t（ODE的自治系统），下面将规避错误：

    h=0.1;
    K1=(h*RHS([],X0));
    K2=h*(RHS([],X0+0.5*K1));
    K3=h*(RHS([],X0-K1+(2*K2)));
    X1=X0+(1/6)*(K1+4*K2+K3);

MATLAB 只是在寻找额外的输入，仅此而已。 注意： 输入一个空向量 [] 没有什么特别的，实际上这个“修复”适用于任何输入（包括 string 或复数).诀窍是占据占位符输入。

Answer 2

另一种变体不是用伪参数扩展，而是减少到

function Y=RHS(X)

那么RK3步骤可以保持原样，而ode45调用需要适配

[t, X]= ode45(@(t,X) RHS(X), [0 3], X0);

运行 matlab 代码表明没有提供足够的输入，并且我们的 RHS 函数有问题

Runing matlab code gives that not enough inputs have been given and that there is something wrong with our RHS function

matlab

numerical-methods

ode

orbital-mechanics