Search code examples
matlabnumerical-methodsodeorbital-mechanics

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


the goal with this code is to simulate numerically the movement of three bodies in three dimensions. All of the bodies have a constant masses. We are going to determinate the evolution of the position of the bodies in tn relation to time. However when we run the code the RHS function malfuctions becasue "not enough inputs have been given" and we cant figure out why.

    % 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 

Solution

  • Thomas:

    I just ran your code on my end and it works flawlessly up to and including the callout to ode45. I did, however, get the same error that you claimed in the question. The error is triggered by the line K1=(h*RHS(X0)), which I believe is the beginning of a block where you're trying to set up a Runge-Kutta method. The reason you get the error is that in that line you only give your RHS function a single input, yet, you've defined RHS to take two inputs (t, and 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 
    

    Since the RHS function does not use t at all (autonomous system of ODE's), the following will circumvent the error:

        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 was just looking for that additional input, that's all. NOTE: There was nothing special about inputting an empty vector [], really this "fix" would have worked with any inputs (including a string or a complex number). The trick was to occupy that placeholder input.