Central difference method to obtain solutions for derived Equations of Motions of UAV aircraft

Having trouble find a way to use the central difference method to obtain solutions for nonlinear derived Equations of Motion (EOM) of a UAV aircraft. I have already derived the Equations of Motions but just not sure if I am doing the coding correctly to obtain the solution for the EOM using the time integration method (Central difference method). Can anyone help please

[John D'Errico]

Are you doing the coding correctly? Is it possible to know, if you don't show what you did? If you want help, then show what you are doing. Otherwise, this just asks someone to write complete code for you.

[Taikiko Taukei]

Hi @John D'Errico so this is my code for the central difference method with example inputs, the code following this is the fourth-order Runge-Kutta method (RK4). Both codes only involves 3-DOF Equations of Motion for Longitudunal dynamics with same parameters for both codes however I am obtaining different results. Why is this?

UAS 3DOF EOM solutions obtained using Central Difference Method code:

% Define parameters

g = 9.81; % gravitational acceleration (m/s^2)

m = 1000; % mass (kg)

S = 10; % wing area (m^2)

c = 1; % mean aerodynamic chord (m)

Iyy = 1000; % moment of inertia (kg*m^2)

CM0 = -0.05; % Pitching moment coefficient at zero angle of attack

CM_alpha = -0.5; % Pitching moment coefficient with respect to angle of attack

CMq = -15; % Pitching moment coefficient with respect to pitch rate

CL0 = 0.1; % Lift coefficient at zero angle of attack

CL_alpha = 5; % Lift coefficient with respect to angle of attack

CLq = 0; % Lift coefficient with respect to pitch rate

CD0 = 0.05; % Drag coefficient at zero angle of attack

CD_alpha = 0.3; % Drag coefficient with respect to angle of attack

CDq = 0;

rho = 1.225;

% Initial conditions

u0 = 15; % initial velocity (m/s)

alpha0 = 0.1; % initial angle of attack (radians)

q0 = 0; % initial pitch rate (rad/s)

% Time parameters

T = 10; % total simulation time (s)

dt = 0.01; % time step (s)

N = ceil(T/dt);

% Initialize arrays

t = zeros(num_steps, 1); % time array

u = zeros(num_steps, 1); % velocity array

alpha = zeros(num_steps, 1);% angle of attack array

q = zeros(num_steps, 1); % pitch rate array

% Set initial conditions

t(1) = 0;

u(1) = u0;

alpha(1) = alpha0;

q(1) = q0;

% Main loop

for i = 2:N

% Calculate derivatives using central difference method

du_dt = (1/(2*dt)) * (u(i-1) - u(i+1));

dalpha_dt = (1/(2*dt)) * (alpha(i-1) - alpha(i+1));

dq_dt = (1/(2*dt)) * (q(i-1) - q(i+1));

% Update state variables using derivatives

u(i) = u(i-1) + dt * du_dt;

alpha(i) = alpha(i-1) + dt * dalpha_dt;

q(i) = q(i-1) + dt * dq_dt;

% Forces and moments

L = 0.5 * rho * u(i)^2 * S * (CL0 + CL_alpha * alpha(i)); % lift force

D = 0.5 * rho * u(i)^2 * S * (CD0 + CD_alpha * alpha(i)); % drag force

M = 0.5 * rho * u(i)^2 * S * c * (CM0 + CM_alpha * alpha(i)); % pitching moment

% Kinematics

du_dt = (1/m) * (L - m*g*cos(alpha(i))); % acceleration

dalpha_dt = q(i); % angular velocity

dq_dt = (1/Iyy) * M; % angular acceleration

% Update state variables

u(i) = u(i-1) + dt * du_dt;

alpha(i) = alpha(i-1) + dt * dalpha_dt;

q(i) = q(i-1) + dt * dq_dt;

% Update time

t(i) = t(i-1) + dt;

end

% Plot results

figure;

subplot(3,1,1);

plot(t, u);

xlabel('Time (s)');

ylabel('Velocity (m/s)');

title('Velocity vs Time');

subplot(3,1,2);

plot(t, alpha * 180/pi);

xlabel('Time (s)');

ylabel('Angle of Attack (degrees)');

title('Angle of Attack vs Time');

subplot(3,1,3);

plot(t, q * 180/pi);

xlabel('Time (s)');

ylabel('Pitch Rate (deg/s)');

title('Pitch Rate vs Time');

[Taikiko Taukei]

@John D'Errico This is the fourth-order Runge-Kutta method (RK4) method code for solving the equations of motion for 3 DOF Longitudinal

function UAS_3DOF_Longitudinal_RK4()

% Parameters

m = 1000; % Mass (kg)

Iyy = 1000; % Moment of inertia about y-axis (kg*m^2)

S = 10; % Wing area (m^2)

c = 1; % Mean aerodynamic chord (m)

rho = 1.225; % Air density (kg/m^3)

g = 9.81; % Gravity (m/s^2)

Cm0 = -0.05; % Pitching moment coefficient at zero angle of attack

Cma = -0.5; % Pitching moment coefficient with respect to angle of attack

Cmq = -15; % Pitching moment coefficient with respect to pitch rate

CL0 = 0.1; % Lift coefficient at zero angle of attack

CLa = 5; % Lift coefficient with respect to angle of attack

CLq = 0; % Lift coefficient with respect to pitch rate

CD0 = 0.05; % Drag coefficient at zero angle of attack

CDa = 0.3; % Drag coefficient with respect to angle of attack

CDq = 0; % Drag coefficient with respect to pitch rate

% Initial conditions

u0 = 30; % Initial velocity (m/s)

alpha0 = 0.1; % Initial angle of attack (radians)

q0 = 0; % Initial pitch rate (rad/s)

% Simulation time

dt = 0.01; % Time step (s)

t_end = 10; % End time (s)

N = ceil(t_end/dt); % Number of time steps

% Initialize arrays to store results

t = zeros(1, N+1);

u = zeros(1, N+1);

alpha = zeros(1, N+1);

q = zeros(1, N+1);

% Initial conditions

t(1) = 0;

u(1) = u0;

alpha(1) = alpha0;

q(1) = q0;

% RK4 integration loop

for i = 1:N

t(i+1) = t(i) + dt;

[du, dalpha, dq] = dynamics(u(i), alpha(i), q(i), m, Iyy, S, c, rho, g, ...

Cm0, Cma, Cmq, CL0, CLa, CLq, CD0, CDa, CDq);

k1u = dt * du;

k1alpha = dt * dalpha;

k1q = dt * dq;

[du, dalpha, dq] = dynamics(u(i) + 0.5*k1u, alpha(i) + 0.5*k1alpha, q(i) + 0.5*k1q, ...

m, Iyy, S, c, rho, g, Cm0, Cma, Cmq, CL0, CLa, CLq, CD0, CDa, CDq);

k2u = dt * du;

k2alpha = dt * dalpha;

k2q = dt * dq;

[du, dalpha, dq] = dynamics(u(i) + 0.5*k2u, alpha(i) + 0.5*k2alpha, q(i) + 0.5*k2q, ...

m, Iyy, S, c, rho, g, Cm0, Cma, Cmq, CL0, CLa, CLq, CD0, CDa, CDq);

k3u = dt * du;

k3alpha = dt * dalpha;

k3q = dt * dq;

[du, dalpha, dq] = dynamics(u(i) + k3u, alpha(i) + k3alpha, q(i) + k3q, ...

m, Iyy, S, c, rho, g, Cm0, Cma, Cmq, CL0, CLa, CLq, CD0, CDa, CDq);

k4u = dt * du;

k4alpha = dt * dalpha;

k4q = dt * dq;

u(i+1) = u(i) + (1/6) * (k1u + 2*k2u + 2*k3u + k4u);

alpha(i+1) = alpha(i) + (1/6) * (k1alpha + 2*k2alpha + 2*k3alpha + k4alpha);

q(i+1) = q(i) + (1/6) * (k1q + 2*k2q + 2*k3q + k4q);

end

% Plot results

figure;

subplot(3,1,1);

plot(t, u);

xlabel('Time (s)');

ylabel('Velocity (m/s)');

title('Velocity vs Time');

subplot(3,1,2);

plot(t, alpha * 180/pi);

xlabel('Time (s)');

ylabel('Angle of Attack (degrees)');

title('Angle of Attack vs Time');

subplot(3,1,3);

plot(t, q * 180/pi);

xlabel('Time (s)');

ylabel('Pitch Rate (deg/s)');

title('Pitch Rate vs Time');

end

function [du, dalpha, dq] = dynamics(u, alpha, q, m, Iyy, S, c, rho, g, ...

Cm0, Cma, Cmq, CL0, CLa, CLq, CD0, CDa, CDq)

% Dynamic equations

CL = CL0 + CLa * alpha + CLq * (c/(2*u)) * q;

CD = CD0 + CDa * alpha + CDq * (c/(2*u)) * q;

Cm = Cm0 + Cma * alpha + Cmq * (c/(2*u)) * q;

du = - (CD * rho * S * u^2)/(2*m) + g * sin(alpha);

dalpha = q;

dq = (1/Iyy) * (Cm * rho * S * c * u^2/2 - m * g * c * sin(alpha));

end

[Torsten]

% Calculate derivatives using central difference method
du_dt = (1/(2*dt)) * (u(i-1) - u(i+1));
dalpha_dt = (1/(2*dt)) * (alpha(i-1) - alpha(i+1));
dq_dt = (1/(2*dt)) * (q(i-1) - q(i+1));

You don't know u(i+1), alpha(i+1) and q(i+1) at this stage of your computation. So you just insert the initialized values for them (namely zero) which is of course incorrect.

Can you give a link to this method for solving ODEs ?

Why don't you use a professional MATLAB ODE integrator, e.g. ode45 ?

[Taikiko Taukei]

@Torsten my investigation is to use the central difference method as a time integration method to solve the 6-DOF of the derived nonlinear equation and compare to a more professional MATLAB ODE integrator such as ode45 but seems like the results aren't showing the same

[Taikiko Taukei]

@Torsten and @John D'Errico the problem I am facing is that every time I try to run the central differece method code this message appears:

Index exceeds the number of array elements. Index must not exceed 1000.

Error in UAS_3DOF_Longitudinal_CDM (line 44)

du_dt = (1/(2*dt)) * (u(i-1) - u(i+1));

[Torsten]

Didn't you read my response ?

You try to use a value in the computation of du_dt that is not yet known, namely u(i+1).

If you hadn't initialized the arrays as

u = zeros(num_steps, 1);    % velocity array
alpha = zeros(num_steps, 1);% angle of attack array
q = zeros(num_steps, 1);    % pitch rate array

the above error message would have popped up already for i = 2. Now it appears for i = N because the above arrays are only initialized up to index N, but you try to access u(N+1).

It would be interesting to know how the central difference method to solve ODEs works correctly. I cannot find it anywhere in the web.