How do I go about plotting points from ODE45 with this?

Question

0 个投票

I have been trying to use ODE45 to solve the linearised model, my aim is to plot values of y and x displacement but the code does not seem to churn out any right answers. The attached images are the equations I am given to work around with. Are my ODE45 parameters wrong? Any advice would help a lot!

my function is as follows

function dydx =  odefcn(x, y, a, v)
dydx = zeros(2, 1);
dydx(1) = y(2);
dydx(2) = a/v^2;
%========================= MAIN CODE ===========================%
% Get launch angle (theta) from user, then get initial x and y coordinates
% Constants are VELOCITY (v), N, DESIGNATED TIME (Td) 
% Variables are x and y 
% LOS, a, L are always changing, so need to INCLUDE THEM IN A LOOP
% Output will be x and y, ALSO ALWAYS CHANGING
%======================= Get input from user =============================%
theta0 = input('Please input launch angle in degrees: ');
x0 = input('Please enter initial x coordinates: '); 
y0 = input('Please enter initial y coordinates: ');
%=========================== Constants ===================================%
v = 300;
N = 3;                          % Proportionality constant
Td = 40;                        % designated impact time
targetx= 50; targety = 50;      % set targetcoordinates 
Rgo_x = targetx - x0;
Rgo_y = targety - y0;
    theta = theta0; 
    xspan = [x0 targetx];
    yspan = [y0 targety];
    IC = [0 theta0]; 
    
    x=x0; y=y0;
    while Rgo_y ~= 0 & Rgo_x ~= 0
    
        Rgo = sqrt(Rgo_x^2 + Rgo_y^2);                                          
            
        LOS = -(Rgo_y/(Rgo^2));                                                                                                    
        Tgo = (1+ (theta-acosd(Rgo_x/Rgo))/10)*(Rgo/v); 
        a = N*v*LOS - (60*v^5*(Td-Tgo))/(N*v*LOS*(Rgo^3)); 
        theta = a*(Rgo_x)/v^2 + theta; 
        
        [x, y] = ode45(@(x, y)odefcn(x, y, a, v), xspan, yspan, IC);            
        Rgo_x = targetx - x;
        Rgo_y = targety - y;
        
    end
%================ Plot x and Y values wrt to angle changes ===============%
figure
grid on
plot(x, y, 'o')
xlabel('x displacement (m)');
ylabel('y displacement (m)');