I am writing a program for lateral- direction motion of an airvraft using Runga Kutta 4th order method.

at the time of compilation i am facing following error.

Array indices must be positive integers or logical values.

Error in ld>@(b,p,r,phi)-r+(g/V)*sin(phi(i))+((qdp*S/m*V)*(Cy0+Cyp*(p*s/V)+Cyr(r*s/V)+Cyb*b+Cyda*dela(i)+Cydr*delr(i))) (line 65)

bdot= @(b,p,r,phi) -r+(g/V)*sin(phi(i))+((qdp*S/m*V)*(Cy0+Cyp*(p*s/V)+Cyr(r*s/V)+Cyb*b+Cyda*dela(i)+Cydr*delr(i)));

Error in ld (line 71)

b1=h*bdot(b(i),p(i),r(i),phi(i));

%%this is my program:
%write the computer programme in Matlab to solve equations of motion of
%lateral-directional motion of an aircraft using Runga Kutta 4th order method(ts1519092001)

clc
clear all
close all
%A/D derivatives  
Cyp=0.303;
Cyr=0.727;
Cyb=-1.333;
Cyda=0.029;
Cydr=0.191;
Clp=-0.978;
Clr=0.418;
Clb=-0.126;
Clda=-0.247;
Cldr=0.046;
Cnp=-0.115;
Cnr=-0.495;
Cnb=0.281;
Cnda=0.0;
Cndr=-0.166;
V=106.7; %velocity
c=3.159; %mean aerodynamic system
rho=0.731; %density
s=10.75; %half wing span
S=64;%wing area
g=9.81; %Gravity
%parameters
m=19633.23; %mass
wf=5023; %Fuel weight 
% Moment of inertia
Ix=104431+6.31626.*wf+0.0021277.*wf.*wf;
Iz=331864+5.53291.*wf+0.0021277.*wf.*wf;
Iy=252687;
Ixz=11442.0;
Cy0=0;
Cl0=0;
Cn0=0;
qdp=rho.*V.*V/2;
h=0.01; % Time step
%Initial Condition
b(1)=0;
p(1)=0;
r(1)=0;
phi(1)=0;
t=0:h:7.0;
dela(1:300)=2.0*pi/180;
dela(301:500)=-2.0*pi/180;
dela(501:600)=1.0*pi/180;
dela(601:700)=-1.0*pi/180;
delr(1:300)=2.0*pi/180;
delr(301:500)=-2.0*pi/180;
delr(501:600)=1.0*pi/180;
delr(601:700)=-1.0*pi/180;
for i=1:(length(t)-1) 
    bdot= @(b,p,r,phi) -r+(g/V)*sin(phi(i))+((qdp*S/m*V)*(Cy0+Cyp*(p*s/V)+Cyr(r*s/V)+Cyb*b+Cyda*dela(i)+Cydr*delr(i)));
    pdot= @(p,phi,b,r) qdp*S*s*((Iz*(Cl0+Clp*(p*s/V)+Cyr*(r*s/V)+Clb*b+Clda*dela(i))+Cldr*delr(i))+(Ixz*(Cn0+Cnp*(p*s/V)+Cnr*(r*s/V)+Cnb*b+Cnda*dela(i)+Cndr*delr(i))))/(Ix*Iz-Ixz*Ixz);
    rdot= @(r,phi,b,p) qdp*S*s*(Ix*(Cn0+Cnp*(p*s/V)+Cnr*(r*s/V)+Cnb*b+Cnda*dela(i)+Cndr*delr(i))+Ixz*(Cl0+Clp*(p*s/V)+Cyr*(r*s/V)+Clb*b+Clda*dela(i))+Cldr*delr(i))/(Ix*Iz-Ixz*Ixz);
    phidot= @(phi,b,p,r) p(i);
    
    
    b1=h*bdot(b(i),p(i),r(i),phi(i));
    p1=h*pdot(p(i),r(i),phi(i),b(i));
    r1=h*rdot(r(i),phi(i),b(i),p(i));
    phi1=h*phidot(phi(i),b(i),p(i),r(i));
    
    b2=h*bdot(b(i)+0.5*h,p(i)+0.5*p1,r(i)+0.5*r1,phi(i)+0.5*phi1);
    p2=h*pdot(p(i)+0.5*h,r(i)+0.5*r1,phi(i)+0.5*phi1,b(i)+0.5*b1);
    r2=h*rdot(r(i)+0.5*h,phi(i)+0.5*phi1,b(i)+0.5*b1,p(i)+0.5*p1);
    phi2=h*phidot(phi(i)+0.5*h,b(i)+0.5*b1,p(i)+0.5*p1,r(i)+0.5*r1);
    
    b3=h*bdot(b(i)+0.5*h,p(i)+0.5*p2,r(i)+0.5*r2,phi(i)+0.5*phi2);
    p3=h*pdot(p(i)+0.5*h,r(i)+0.5*r2,phi(i)+0.5*phi2,b(i)+0.5*b2);
    r3=h*rdot(r(i)+0.5*h,phi(i)+0.5*phi2,b(i)+0.5*b2,p(i)+0.5*p2);
    phi3=h*phidot(phi(i)+0.5*h,b(i)+0.5*b2,p(i)+0.5*p2,p(i)+0.5*p2);
    
    b4=h*bdot(b(i)+h,p(i)+p3,r(i)+r3,phi(i)+phi3);
    p4=h*pdot(p(i)+h,r(i)+r3,phi(i)+phi3,b(i)+b3);
    r4=h*rdot(r(i)+h,phi(i)+phi3,b(i)+b3,p(i)+p3);
    phi4=h*phidot(phi(i)+h,b(i)+b3,p(i)+p3,r(i)+r3);
    b(i+1)=b(i)+(1/6)*(b1+b2+b3+b4);
    p(i+1)=p(i)+(1/6)*(p1+p2+p3+p4);
    r(i+1)=r(i)+(1/6)*(r1+r2+r3+r4);
    phi(i+1)=phi(i)+(1/6)*(phi1+phi2+phi3+phi4);
    
    
end
dela(length(dela)+1)=dela(length(dela));
delr(length(delr)+1)=delr(length(delr));
Cy= (Cy0.*ones(1,701))+(Cyp.*(p.*s/V))+(Cyr.*(r.*s/V))+(Cyb.*b)+(Cyda.*dela)+(Cydr.*delr);
Cl= (Cl0.*ones(1,701))+(Clp.*(p.*s/V))+(Clr.*(r.*s/V))+(Clb.*b)+(Clda.*dela)+(Cldr.*delr);
Cn= (Cn0.*ones(1,701))+(Cnp.*(p.*s/V))+(Cnr.*(r.*s/V))+(Cnb.*b)+(Cnda.*dela)+(Cndr.*delr);
subplot(6,1,1)
plot(t,(b*180)/pi,'.'); 
hold on 
xlabel('Time')
ylabel('beta')
subplot(6,1,2)
plot(t,(p*180)/pi,'.');
hold on
xlabel('Time')
ylabel('p')
subplot(6,1,3)
plot(t,(r*180)/pi,'.');
hold on
xlabel('Time')
ylabel('r')
subplot(6,1,4)
plot(t,(phi*180)/pi,'.');
hold on
xlabel('Time')
ylabel('phi')
subplot(6,1,5)
plot(t,(dela*180/pi),'.');
hold on 
xlabel('Time')
ylabel('delta-a')
subplot(6,1,6)
plot(t,(delr*180/pi),'.');
hold on 
xlabel('Time')
ylabel('delta-r')

[Torsten]

And why don't you use one of MATLAB's ODE integrators (e.g. ODE45) where Runge-Kutta-4-5 is professionally implemented ?

[James Tursa]

I agree with Torsten. It would be much better to set up your code using a state vector approach. This will greatly simplify your code and allow you to use MATLAB supplied integrators such as ode45( ).

Answer 1

[Walter Roberson]

Cyr(r*s/V)

is indexing Cyr at the location in the expression. You probably wanted to multiply.

Answer 2

[Torsten]

%%this is my program:
%write the computer programme in Matlab to solve equations of motion of
%lateral-directional motion of an aircraft using Runga Kutta 4th order method(ts1519092001)

clc
clear all
close all
%A/D derivatives  
Cyp=0.303;
Cyr=0.727;
Cyb=-1.333;
Cyda=0.029;
Cydr=0.191;
Clp=-0.978;
Clr=0.418;
Clb=-0.126;
Clda=-0.247;
Cldr=0.046;
Cnp=-0.115;
Cnr=-0.495;
Cnb=0.281;
Cnda=0.0;
Cndr=-0.166;
V=106.7; %velocity
c=3.159; %mean aerodynamic system
rho=0.731; %density
s=10.75; %half wing span
S=64;%wing area
g=9.81; %Gravity
%parameters
m=19633.23; %mass
wf=5023; %Fuel weight 
% Moment of inertia
Ix=104431+6.31626.*wf+0.0021277.*wf.*wf;
Iz=331864+5.53291.*wf+0.0021277.*wf.*wf;
Iy=252687;
Ixz=11442.0;
Cy0=0;
Cl0=0;
Cn0=0;
qdp=rho.*V.*V/2;
h=0.01; % Time step
%Initial Condition
b(1)=0;
p(1)=0;
r(1)=0;
phi(1)=0;
t=0:h:7.0;
dela(1:300)=2.0*pi/180;
dela(301:500)=-2.0*pi/180;
dela(501:600)=1.0*pi/180;
dela(601:700)=-1.0*pi/180;
delr(1:300)=2.0*pi/180;
delr(301:500)=-2.0*pi/180;
delr(501:600)=1.0*pi/180;
delr(601:700)=-1.0*pi/180;
for i=1:(length(t)-1) 
    bdot= @(b,p,r,phi) -r+(g/V)*sin(phi)+((qdp*S/m*V)*(Cy0+Cyp*(p*s/V)+Cyr*(r*s/V)+Cyb*b+Cyda*dela(i)+Cydr*delr(i)));
    pdot= @(p,phi,b,r) qdp*S*s*((Iz*(Cl0+Clp*(p*s/V)+Cyr*(r*s/V)+Clb*b+Clda*dela(i))+Cldr*delr(i))+(Ixz*(Cn0+Cnp*(p*s/V)+Cnr*(r*s/V)+Cnb*b+Cnda*dela(i)+Cndr*delr(i))))/(Ix*Iz-Ixz*Ixz);
    rdot= @(r,phi,b,p) qdp*S*s*(Ix*(Cn0+Cnp*(p*s/V)+Cnr*(r*s/V)+Cnb*b+Cnda*dela(i)+Cndr*delr(i))+Ixz*(Cl0+Clp*(p*s/V)+Cyr*(r*s/V)+Clb*b+Clda*dela(i))+Cldr*delr(i))/(Ix*Iz-Ixz*Ixz);
    phidot= @(phi,b,p,r) p;
    
    
    b1=h*bdot(b(i),p(i),r(i),phi(i));
    p1=h*pdot(p(i),r(i),phi(i),b(i));
    r1=h*rdot(r(i),phi(i),b(i),p(i));
    phi1=h*phidot(phi(i),b(i),p(i),r(i));
    
    b2=h*bdot(b(i)+0.5*h,p(i)+0.5*p1,r(i)+0.5*r1,phi(i)+0.5*phi1);
    p2=h*pdot(p(i)+0.5*h,r(i)+0.5*r1,phi(i)+0.5*phi1,b(i)+0.5*b1);
    r2=h*rdot(r(i)+0.5*h,phi(i)+0.5*phi1,b(i)+0.5*b1,p(i)+0.5*p1);
    phi2=h*phidot(phi(i)+0.5*h,b(i)+0.5*b1,p(i)+0.5*p1,r(i)+0.5*r1);
    
    b3=h*bdot(b(i)+0.5*h,p(i)+0.5*p2,r(i)+0.5*r2,phi(i)+0.5*phi2);
    p3=h*pdot(p(i)+0.5*h,r(i)+0.5*r2,phi(i)+0.5*phi2,b(i)+0.5*b2);
    r3=h*rdot(r(i)+0.5*h,phi(i)+0.5*phi2,b(i)+0.5*b2,p(i)+0.5*p2);
    phi3=h*phidot(phi(i)+0.5*h,b(i)+0.5*b2,p(i)+0.5*p2,p(i)+0.5*p2);
    
    b4=h*bdot(b(i)+h,p(i)+p3,r(i)+r3,phi(i)+phi3);
    p4=h*pdot(p(i)+h,r(i)+r3,phi(i)+phi3,b(i)+b3);
    r4=h*rdot(r(i)+h,phi(i)+phi3,b(i)+b3,p(i)+p3);
    phi4=h*phidot(phi(i)+h,b(i)+b3,p(i)+p3,r(i)+r3);
    b(i+1)=b(i)+(1/6)*(b1+b2+b3+b4);
    p(i+1)=p(i)+(1/6)*(p1+p2+p3+p4);
    r(i+1)=r(i)+(1/6)*(r1+r2+r3+r4);
    phi(i+1)=phi(i)+(1/6)*(phi1+phi2+phi3+phi4);
    
    
end
dela(length(dela)+1)=dela(length(dela));
delr(length(delr)+1)=delr(length(delr));
Cy= (Cy0.*ones(1,701))+(Cyp.*(p.*s/V))+(Cyr.*(r.*s/V))+(Cyb.*b)+(Cyda.*dela)+(Cydr.*delr);
Cl= (Cl0.*ones(1,701))+(Clp.*(p.*s/V))+(Clr.*(r.*s/V))+(Clb.*b)+(Clda.*dela)+(Cldr.*delr);
Cn= (Cn0.*ones(1,701))+(Cnp.*(p.*s/V))+(Cnr.*(r.*s/V))+(Cnb.*b)+(Cnda.*dela)+(Cndr.*delr);
subplot(6,1,1)
plot(t,(b*180)/pi,'.'); 
hold on 
xlabel('Time')
ylabel('beta')
subplot(6,1,2)
plot(t,(p*180)/pi,'.');
hold on
xlabel('Time')
ylabel('p')
subplot(6,1,3)
plot(t,(r*180)/pi,'.');
hold on
xlabel('Time')
ylabel('r')
subplot(6,1,4)
plot(t,(phi*180)/pi,'.');
hold on
xlabel('Time')
ylabel('phi')
subplot(6,1,5)
plot(t,(dela*180/pi),'.');
hold on 
xlabel('Time')
ylabel('delta-a')
subplot(6,1,6)
plot(t,(delr*180/pi),'.');
hold on 
xlabel('Time')
ylabel('delta-r')

[TUSHAR SEN]

thank you man, i understood the error now and now its working.

[Sam Chak]

@TUSHAR SEN, just want to add that the results should not blow up to .

Which are the aircraft control surfaces?

What mathematical equations did you use to describe the actions of the control surfaces?

Answer 3

[Sam Chak]

Hi @TUSHAR SEN

This is just a test. The angles of deflection for the Ailerons and Rudder are set to 0. And it turns out that the system is stable. If you want to make {} converge to some desired values, then you need to carefully design the tracking laws for and . However, you only have two actuators to manipulate the 3 main states {}.

So, the system is underactuated, and this indicates that it is a little challenging design problem for humans, unless you use some optimizers, or iterative learning, or reinforcement learning to search in the solution spaces and .

% Settings
tspan  = linspace(0, 30, 3001);
b0     = 0.1;
p0     = 0.2;
r0     = 0.3;
phi0   = 0.4;
x0     = [b0 p0 r0 phi0];

% ODE solver
[t, x] = ode45(@odefcn, tspan, x0);

% Plot
plot(t, x), grid on, xlabel('\it{t}')
legend('\it{b}(\it{t})', '\it{p}(\it{t})', '\it{r}(\it{t})', '\phi(\it{t})', 'location', 'east')

function dxdt = odefcn(t, x)

    dxdt = zeros(4, 1);
    % parameters
    Cyp  =  0.303;
    Cyr  =  0.727;
    Cyb  = -1.333;
    Cyda =  0.029;
    Cydr =  0.191;
    Clp  = -0.978;
    Clr  =  0.418;
    Clb  = -0.126;
    Clda = -0.247;
    Cldr =  0.046;
    Cnp  = -0.115;
    Cnr  = -0.495;
    Cnb  =  0.281;
    Cnda =  0.0;
    Cndr = -0.166;
    
    V    =  106.7;      % velocity
    c    =  3.159;      % mean aerodynamic system
    rho  =  0.731;      % density
    s    =  10.75;      % half wing span
    S    =  64;         % wing area
    g    =  9.81;       % Gravity
    m    =  19633.23;   % mass
    wf   =  5023;       % Fuel weight 
    
    % Moment of inertia
    Ix   =  104431 + 6.31626*wf + 0.0021277*wf^2;
    Iz   =  331864 + 5.53291*wf + 0.0021277*wf^2;
    Iy   =  252687;
    Ixz  =  11442.0;
    
    Cy0  = 0;
    Cl0  = 0;
    Cn0  = 0;
    qdp  = rho*(V^2)/2;
    
    % flight control surfaces
    % https://en.wikipedia.org/wiki/Flight_control_surfaces
    dela    = 0;    % Aileron angle of deflection 
    delr    = 0;    % Rudder angle of deflection 
    
    b       = x(1);
    p       = x(2);
    r       = x(3);
    phi     = x(4);
    dxdt(1) = - r + (g/V)*sin(phi) + (qdp*S/m*V)*(Cy0 + Cyp*(p*s/V) + Cyr*(r*s/V) + Cyb*b + Cyda*dela + Cydr*delr);
    dxdt(2) = qdp*S*s*((Iz*(Cl0 + Clp*(p*s/V) + Cyr*(r*s/V) + Clb*b + Clda*dela) + Cldr*delr) + (Ixz*(Cn0 + Cnp*(p*s/V) + Cnr*(r*s/V) + Cnb*b + Cnda*dela  + Cndr*delr)))/(Ix*Iz - Ixz*Ixz);
    dxdt(3) = qdp*S*s*( Ix*(Cn0 + Cnp*(p*s/V) + Cnr*(r*s/V) + Cnb*b + Cnda*dela  + Cndr*delr) +  Ixz*(Cl0 + Clp*(p*s/V) + Cyr*(r*s/V) + Clb*b + Clda*dela) + Cldr*delr)  /(Ix*Iz - Ixz*Ixz);
    dxdt(4) = p;

end