16.901 Project #1
Sample Solution
1
Background
Modeling the dynamics of an airplane can be very complex. Twelve state variables are required to describe
the motion: three states for position, three for velocity, three for angular orientation, and three angular rates.
However, this system of twelve coupled diﬀerential equations can be simpliﬁed greatly by linearization and
some assumptions. These simpliﬁcations yield two decoupled four-state systems, one describing longitudinal
motion and the other describing lateral motion.
The longitudinal modes can be simpliﬁed even more to give fairly accurate decoupled models of the
familiar pheugoid and short-period modes. The lateral system cannot be simpliﬁed as easily to solve for
the three lateral modes: Dutch roll, spiral mode, and roll mode. The Dutch roll mode is a combination of
yawing and rolling oscillations. In the roll mode, the roll rate reaches steady state quickly. The spiral mode
can be slightly stable or unstable. An unstable spiral mode can lead to divergence from the ﬂight path or
result in a spiral dive. These lateral modes will be the topic of the ﬁrst programming project and will be
the basis for most of the questions in this homework assignment.
1.1
Deﬁnition of coordinate systems and variables
The motion of the aircraft is measured relative to a ﬁxed frame; however, the properties of the aircraft are
often known in a coordinate system relative to the aircraft body. Thus, the states to be determined are the
the relative position, motion, angular orientation, and angular rate of the body frame. Figure 1 shows the
forces (X, Y , Z), moments (L, M , N ), angular rates of rotation (p, q, r) with respect to the body axes (x b ,
yb , zb ).
Figure 1: Deﬁnition of body coordinate system attached to the aircraft
The Table 1 deﬁnes the relevant variables. The linearized equations of motion are written in terms of
stability derivatives, the ﬁrst derivatives of the forces and moments with respect to the states. These are
generally written in subscript notation. For example, Y β ≡ ∂Y
∂β . Often stability derivatives are given in a
1
States
Forces and Moments
Other quantities
Symbol
Δβ
Δp
Δr
Δφ
Y
L
N
Q
S
b
W
m
Ix , Iz , Ixz
g
M0
h0
u0
Description
Side-slip angle perturbation
Roll rate perturbation
Yaw rate perturbation
Roll angle perturbation
Side force (Force in y direction)
Rolling moment (Moment about x axis)
Yawing moment (Moment about z axis)
Dynamic pressure
Planform area
Wingspan
Weight of aircraft
Mass of aircraft
Momentw of inertia
Gravitational acceleration
Mach number
altitude
Initial velocity in x direction
Table 1: Explanation of Variables
Yβ
=
QSCyβ ,
Lβ
=
QSbClβ ,
2
Nβ
=
QSbCnβ
Yp
=
QSb
2u0 Cyp ,
Lp
=
QSb
2u0
Clp ,
Np
=
QSb2
2u0 Cnp
Yr
=
QSb
2u0 Cyr ,
Lr
=
QSb2
2u0 Clr ,
Nr
=
QSb2
2u0 Cnr
Table 2: Stability derivative non-dimensional deﬁnitions
non-dimensional form, such as C yβ as the non-dimensionalized form of Y β . Speciﬁcally, the deﬁnitions in
Table 2 are conventional and will be used in this homework and project. Values of the stability derivatives
and dimensions for the airplanes we will study in this homework and project are given in Table 3.
1.2
Equations of motion
The linearized equations of lateral motion are given in Equation (1)-(4). Equation (1) is the conservation of
y-momentum. Equation (2) is the conservation of x-angular momentum. Equation (3) is the conservation
of z-angular momentum. Finally, Equation (4) is the relation between the roll angle and the roll rate.
Speciﬁcally, the governing equations take the following form:
�
�
d
mu0 − Yβ Δβ − Yp Δp + (mu0 − Yr ) Δr − mgΔφ = 0
(1)
dt
�
�
�
�
d
d
−Lβ Δβ + Ix − Lp Δp − Ixz + Lr Δr = 0
(2)
dt
dt
�
�
�
�
d
d
−Nβ Δβ − Ixz + Np Δp + Iz − Nr Δr = 0
(3)
dt
dt
d
Δφ = Δp
dt
2
(4)
Quantity
Cyβ
Cyp
Cyr
Clβ
Clp
Clr
Cnβ
Cnp
Cnr
S
b
h
M0
W
Ix
Iz
Ixz
Values for 747
-0.96
0
0
-0.221
-0.45
0.101
0.150
-0.121
-0.30
5500 f t 2
195.68 f t
Sea level
0.25
636600 lbs
18.2 × 10 6 slug f t2
49.7 × 10 6 slug f t2
0.97 × 10 6 slug f t2
Values for F-104
-1.17
0
0
-0.175
-0.285
0.265
0.50
-0.14
-0.75
196.1 f t2
21.94 f t
Sea level
0.257
16300 lbs
3549 slug f t 2
59669 slug f t 2
0
Table 3: Stability derivatives and dimensions for 747 and F-104 at sea level
2
Assignment
• Implement the following six integration methods:
– 1st order Adams-Bashforth (AB1). More commonly known as Forward Euler.
– 2nd order Adams-Bashforth (AB2).
– 1st order Adams-Moulton (AM1). More commonly known as Backward Euler.
– 2nd order Adams-Moulton (AM2). More commonly known as Trapezoidal Method.
– 2-stage Runge-Kutta (RK2).
– 4-stage Runge-Kutta (RK4).
Turn in hard copies of your algorithms with your project write-up.
®
Solution: See the Matlab scripts for the diﬀerent methods which were distributed with this solution
on the class website.
We will study the accuracy and eﬃciency of these 6 methods applied to the following problem: Find
u(t) from t = 0s to t = 30s for the following initial condition:
Δβ (0) = 0.1 rad, Δp (0) = 0, Δr (0) = 0, Δφ (0) = 0
NOTE: the ﬁnal time has been shortened compared to Homework #2.
• For this assignment, the measure of accuracy will be the maximum error in the sideslip angle pertur
bation (at any iteration), i.e.,
�
�
�
�
E ≡ max �Δβˆn − Δβ(nΔt)�
n
where Δβˆn is the value of the sideslip angle perturbation calculated by the integration method at
the n-th iteration. For both airplanes and using all integration methods, determine the accuracy for
3
Δt
1.0
0.1
0.01
0.001
AB1
1.23e+1
3.89e-2
2.86e-3
2.78e-4
AB2
1.09e-0
1.60e-3
1.59e-5
1.60e-7
AM1
6.20e-2
2.04e-2
2.68e-3
2.76e-4
AM2
2.93e-2
3.20e-4
3.20e-6
3.20e-8
RK2
7.95e-2
6.40e-4
6.40e-6
6.40e-8
RK4
1.44e-3
1.49e-7
1.49e-11
3.77e-13
Table 4: Accuracy convergence study for Boeing 747
Δt
1.0
0.1
0.01
0.001
AB1
4.04e+9
3.11e+2
5.72e-1
4.26e-2
AB2
1.55e+13
8.57e-1
7.76e-3
7.75e-5
AM1
6.23e-1
6.59e-1
3.07e-1
4.00e-2
AM2
7.22e-1
1.52e-1
1.55e-3
1.55e-5
RK2
2.66e+10
3.21e-1
3.10e-3
3.10e-5
RK4
6.23e-1
6.71e-4
6.70e-8
2.57e-11
Table 5: Accuracy convergence study for F-104
timesteps of Δt = 0.001, 0.01, 0.1, and 1.0 seconds. For both aircraft, please discuss which methods
exhibit their predicted global order of accuracy as the timestep decreases? Explain your answer.
Solution: The results of the accuracy study are given in Tables 4 and 5. For a method with order
of accuracy, p, the error E will be bounded by a quantity which is O(Δt p ). Thus, for two diﬀerent
timesteps Δt1 and Δt2 , we would expect the ratio of the errors to be,
�
�p
E(Δt2 )
Δt2
=
.
(5)
Δt1
E(Δt1 )
For the accuracy study performed, the ratio of timesteps is 1/10, thus, we expect the ratio of the errors
to be (1/10)p . But, it is important to remember that this error ratio will only be observed in the limit
of small enough Δt since the accuracy estimate is based on a Taylor series analysis which is only valid
for small Δt.
We begin looking at the Forward Euler results (AB1) for the Boeing 747 in Table 4. AB1 is a ﬁrst-order
accurate method (p = 1), therefore, we expect the error to reduce by a factor of 10. For the largest
timestep, Δt = 1.0 seconds, we note that the error is quite large. In fact, based on the analysis from
Homework #2, we expect Forward Euler to be unstable for a Δt > 0.1 seconds. Thus, at the largest
timestep, Forward Euler is unstable and we expect the answer to be highly inaccurate as the error
will increase at every iteration. In decreasing the timestep from Δt = 1 to Δt = 0.1, the error is
decreased by a factor greater than 10, which is most likely a direct reﬂection that the largest timestep
was unstable and therefore the error was large. The error reduction from Δt = 0.1 to 0.01 is quite
close to 10, and even more so for the ﬁnest timestep. Thus, AB1 is achieving the asymptotic rate by
the ﬁnest timesteps.
As another example, we consider the AM2 results for the Boeing 747. In this case, the method is p = 2
order-accurate and we expect decrease of 100 in the error. This is clearly observed in the AM2 results
contained in Table 4.
The eﬀect of the ﬁnite precision of the computer can also be observed in the RK4 algorithm. In this
case, the accuracy is fourth-order accurate and we expect the error to reduce by a factor of 10 4 . This
reduction is clearly observed for all but the ﬁnest timestep. In decreasing the timestep from 0.01 to
0.001, the error reduces by only a factor of about 40. Because the computer has only a ﬁnite precision,
eventually the accuracy cannot be improved by lowering the timestep. Clearly, though, the solution is
quite accurate at this point.
The results for the F-104 are quite similar in general to the Boeing 747 results. However, the overall
accuracy is decreased. Also, the AB1 and AM1 results never quite achieve asymptotic behavior. This
4
WU
AB1
1.00
AB2
1.14
AM1
1.05
AM2
1.98
RK2
1.88
RK4
4.45
Table 6: Work Units for a single iteration of each integration method.
slight decrease in accuracy is a reﬂection of the unstable physical system which causes the solution to
have larger amplitudes and hence larger errors than in the stable Boeing 747 simulations.
• In order to compare the work involved in these methods, we need a single ’currency’ which is relatively
independent of the particular computer, network, version of software compiler, etc. We will use the
computational time required to run a single iteration of the Forward Euler method (AB1) as our Work
Unit (WU). To ﬁnd the WU’s for an iteration of the other methods, run each of them for a few hundred
iterations and calculate the average time required for a single iteration. MAKE SURE TO DO THIS
ON THE SAME COMPUTER, IMMEDIATELY AFTER EACH OTHER to reduce the possibility of
®
varying conditions aﬀecting the timings. In Matlab, the command to use to ﬁnd the current time is
®
called cputime (look up this command in the Matlab help to see its usage). Then, normalize these
timings by the Forward Euler timings to calculate the WU for an iteration of each method. Report
the WU for an iteration of each method in a table such as Table 6. Note: refer to the discussion
in the sample solution to Homework #2 to see what the expected behavior should be. If you get
something very diﬀerent from the discussion in that sample solution, you probably have a bug or have
not implemented your method eﬃciently.
Solution: All of the implicit methods were implemented by inverting any matrices once-and-for-all
prior to beginning the iteration loop. This was possible because the system being solved is linear with
a constant matrix A. As a result, the costs of the diﬀerent algorithms would be largely a function of the
number of matrix-vector multiplications required per iteration. For AB1, AB2, and AM1, only a single
matrix-vector is required at any given iteration. Note, AB2 does not require 2 matrix-vector multiplies
in my implementation because the forcing function at n − 1 can be stored and not re-calculated. AM2
requires a matrix-vector multiply for both the implicit and explicit parts of the operator. Finally, each
Runge-Kutta method requires a matrix-vector multiply at every stage. The Work Units reﬂect this
implementation and are shown in Table 6. Note: to better understand this discussion, it may be useful
®
to look at the implementation in the Matlab scripts.
• Based on the results in the tables for the accuracy study and the known computational work costs
of each algorithm, estimate the amount of work required (i.e. in terms of WU) for each method to
achieve an accuracy of 0.0001 radians for β for each airplane (you might use a table like Table 7. Try to
explain the WU results based on your understanding of the methods and their accuracy and stability.
Which method would you recommend for eﬃciently solving these aircraft lateral dynamics equations?
Solution: To estimate the work required to achieve a ﬁxed accuracy, we need to interpolate the
accuracy results from the previous results. To do this, we will assume that all of the methods have
achieved their asymptotic order and then extrapolate from the nearest error level which was just greater
than the 0.0001 requirement. For example, let’s consider the RK2 method for the Boeing 747 results.
As can be observed from Table 4, the RK2 accuracy for a Δt = 0.1 is E = 6.4 × 10 −4 . We can solve
for the Δt to achieve the desired accuracy by re-arranging Equation (5),
�
E(Δt2 )
Δt2 =
E(Δt1 )
�1/p
Δt1 .
So, setting Δt1 = 0.1 seconds, E(Δt1 ) = 6.4 × 10 −4 , E(Δt2 ) = 1 × 10−4 , and p = 2, we ﬁnd that
Δt2 = 0.0395 seconds. Then, to ﬁnd the work required, we calculate the number of iterations required
5
Aircraft
747
F-104
AB1
8.34 × 10 4
1.28 × 10 7
AB2
1.37 × 10 3
3.01 × 10 4
AM1
8.69 × 10 4
1.26 × 10 7
AM2
1.06 × 10 3
2.34 × 10 4
RK2
1.43 × 10 3
3.14 × 10 4
RK4
2.60 × 10 2
2.15 × 10 3
Table 7: Work Unit requirements to achieve an accuracy of 0.0001 rads for the sideslip angle.
and multiply by the work units per iteration for the method. Thus, for this RK2 example, we ﬁnd that
W URK2:747 =
T
30
W URK2 =
1.88W U = 1427W U
Δt
0.0395
The results for the other methods and both aircraft are given in Table 7.
As can be observed from the results in Table 7, the RK4 method is approximately an order of magnitude
less work than any of the other methods, and therefore is the clear choice for eﬃcient integration of
the aircraft lateral dynamics. The reason for this is clearly the higher order accuracy of the method,
with only a slight increase in the work over the other methods. The results would be even more in
favor of RK4 is higher accuracy were required.
MATLAB® is a trademark of The MathWorks, Inc.
6