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

3 views (last 30 days)
TUSHAR SEN on 11 Sep 2022
Edited: Sam Chak on 14 Sep 2022
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')
##### 2 CommentsShowHide 1 older comment
James Tursa on 14 Sep 2022
Edited: James Tursa on 14 Sep 2022
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( ).

Walter Roberson on 11 Sep 2022
Cyr(r*s/V)
is indexing Cyr at the location in the expression. You probably wanted to multiply.

Torsten on 11 Sep 2022
%%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')
##### 2 CommentsShowHide 1 older comment
Sam Chak on 12 Sep 2022
Edited: Sam Chak on 12 Sep 2022
@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?

Sam Chak on 14 Sep 2022
Edited: Sam Chak on 14 Sep 2022
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

### Categories

Find more on 3DOF in Help Center and File Exchange

R2022a

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!