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

3 views (last 30 days)
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 Comments
James Tursa
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( ).

Sign in to comment.

Answers (3)

Walter Roberson
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
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 Comments
Sam Chak
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?

Sign in to comment.


Sam Chak
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

Tags

Products


Release

R2022a

Community Treasure Hunt

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

Start Hunting!