103,99 €
Mehr erfahren.
- Herausgeber: John Wiley & Sons
- Kategorie: Wissenschaft und neue Technologien
- Serie: Wiley-IEEE Press Book Series on Control Systems Theory and Applications
- Sprache: Englisch
Dynamic System Modeling & Analysis with MATLAB & Python A robust introduction to the advanced programming techniques and skills needed for control engineering In Dynamic System Modeling & Analysis with MATLAB & Python: For Control Engineers, accomplished control engineer Dr. Jongrae Kim delivers an insightful and concise introduction to the advanced programming skills required by control engineers. The book discusses dynamic systems used by satellites, aircraft, autonomous robots, and biomolecular networks. Throughout the text, MATLAB and Python are used to consider various dynamic modeling theories and examples. The author covers a range of control topics, including attitude dynamics, attitude kinematics, autonomous vehicles, systems biology, optimal estimation, robustness analysis, and stochastic system. An accompanying website includes a solutions manual as well as MATLAB and Python example code. Dynamic System Modeling & Analysis with MATLAB & Python: For Control Engineers provides readers with a sound starting point to learning programming in the engineering or biology domains. It also offers: * A thorough introduction to attitude estimation and control, including attitude kinematics and sensors and extended Kalman filters for attitude estimation * Practical discussions of autonomous vehicles mission planning, including unmanned aerial vehicle path planning and moving target tracking * Comprehensive explorations of biological network modeling, including bio-molecular networks and stochastic modeling * In-depth examinations of control algorithms using biomolecular networks, including implementation Dynamic System Modeling & Analysis with MATLAB & Python: For Control Engineers is an indispensable resource for advanced undergraduate and graduate students seeking practical programming instruction for dynamic system modeling and analysis using control theory.
Sie lesen das E-Book in den Legimi-Apps auf:
Seitenzahl: 275
Veröffentlichungsjahr: 2022
Ähnliche
Table of Contents
Cover
Title Page
Copyright
Dedication
Preface
Acknowledgements
Acronyms
About the Companion Website
1 Introduction
1.1 Scope of the Book
1.2 Motivation Examples
1.3 Organization of the Book
Exercises
Bibliography
2 Attitude Estimation and Control
2.1 Attitude Kinematics and Sensors
2.2 Attitude Estimation Algorithm
2.3 Attitude Dynamics and Control
Exercises
Bibliography
Notes
3 Autonomous Vehicle Mission Planning
3.1 Path Planning
3.2 Moving Target Tracking
3.3 Tracking Algorithm Implementation
Exercises
Bibliography
4 Biological System Modelling
4.1 Biomolecular Interactions
4.2 Deterministic Modelling
4.3 Biological Oscillation
Exercises
Bibliography
Notes
5 Biological System Control
5.1 Control Algorithm Implementation
5.2 Robustness Analysis:
‐Analysis
Exercises
Bibliography
6 Further Readings
6.1 Boolean Network
6.2 Network Structure Analysis
6.3 Spatial‐Temporal Dynamics
6.4 Deep Learning Neural Network
6.5 Reinforcement Learning
Bibliography
Appendix A Solutions for Selected Exercises
A.1 Chapter 1
A.2 Chapter 2
A.3 Chapter 3
A.4 Chapter 4
A.5 Chapter 5
Index
Titles in the IEEE Press Series on Control Systems Theory and Applications
End User License Agreement
List of Tables
Chapter 2
Table 2.1 Translational and rotational motions.
Chapter 4
Table 4.1
E. coli
enzyme responses for three different environmental C...
Table 4.2 The nominal values of the model parameters in the tryptophan regul...
Table 4.3 Experimentally known ranges of the parameters from the Supplementa...
Table 4.4 Three optimal parameters changed the most by the genetic algorithm...
Table 4.5 Four optimal parameters changed the most by the differential evolu...
Table 4.6 Laub–Loomis model kinetic parameters.
List of Illustrations
Chapter 1
Figure 1.1 Free‐falling object.
Figure 1.2 Free‐falling object position, velocity, and mass time histories....
Figure 1.3 Ligand–receptor interactions form ligand–receptor complex.
Figure 1.4 (Matlab) EGFR receptor, ligand, and complex time histories.
Figure 1.5 (Python) EGFR receptor, ligand, and complex time histories.
Figure 1.6 The time histories of (a) position (
), (b) velocity (
), and (c)...
Chapter 2
Figure 2.1 Single‐axis rotation about the axis perpendicular to the surface ...
Figure 2.2 The two‐axis rotation is equivalent to the single‐axis rotation a...
Figure 2.3 The object is tumbling in a three‐dimensional space, where the bo...
Figure 2.4 The quaternion time history for
given by (2.9).
Figure 2.5 The time history of the quaternion unit norm error, (2.10), for t...
Figure 2.6 The quaternion time history for
given by (2.9).
Figure 2.7 The time history of the quaternion unit norm error, (2.10), for t...
Figure 2.8 The probability density function of the random number generated b...
Figure 2.9 Five realizations of
whose mean and variance are given by (2.23...
Figure 2.10 Compare the true
and
with the estimated values using 1000 re...
Figure 2.11 The estimated pdf,
, shows the complete picture of the Gaussian...
Figure 2.12 The estimated pdf,
, plot using Python.
Figure 2.13 Bias noise simulation using (a) MATLAB/(b) Python.
Figure 2.14 (MATLAB) Gyroscope measurement simulation.
Figure 2.15 Gyroscope measurement simulation using Python Program 2.14.
Figure 2.16 Identified star in the reference and the bode coordinates.
Figure 2.17 The shortest path rotation.
Figure 2.18 Stochastic simulations using the ODE solvers for
seconds in (a...
Figure 2.19 The Kalman filter for the mass‐spring‐damper system.
Figure 2.20 The Kalman filter states error and 3
bounds for the mass‐spring...
Figure 2.21 Attitude dynamics and kinematics solutions.
Figure 2.22 Quadcopter UAV with the four actuators, where the body frame and...
Figure 2.23 Quadcopter motor force and torques.
Figure 2.24 (a) Attitude stabilization using the quaternion feedback control...
Figure 2.25 (a) Total propeller force and the torque for each direction in t...
Figure 2.26 (a) Settling time distribution over the moment of inertia uncert...
Figure 2.27 CPU loading shows four of them running in 100%.
Figure 2.28 The star sensor frame and the body frame.
Chapter 3
Figure 3.1 Three repulsive potential functions and one attractive potential ...
Figure 3.2 (a) A path converges to the local equilibrium point; (b) the hori...
Figure 3.3 Two equally probably paths shown in (a) and (b), in which the obs...
Figure 3.4 Two equally probably paths shown in (a) and (b), in which the obs...
Figure 3.5 A graph construction for the path planning.
Figure 3.6 (MATLAB) Plot the graph and the shortest path.
Figure 3.7 (Python) Plot the graph and the shortest path.
Figure 3.8 Three repulsive potential functions and one attractive potential ...
Figure 3.9 Random points and Voronoi diagram.
Figure 3.10 Edges inside the obstacles removed and add start and destination...
Figure 3.11 Shortest path in the graph constructed using Voronoi diagram.
Figure 3.12 (a) Original optimal path in the dash‐dotted line and the update...
Figure 3.13 The higher altitude and the closer distance relative to the targ...
Figure 3.14 Global (
–
) and UAV local (
–
) coordinates with the control in...
Figure 3.15 MATLAB
pretty()
function output.
Figure 3.16 UAV target tracking in two‐step prediction.
Figure 3.17 UAV target tracking in two‐step prediction in the mission coordi...
Figure 3.18 Feasible control input space indicated by the painted region.
Figure 3.19 Calculation examples of feasible control input space.
Figure 3.20 Optimal control input for the target tracking.
Figure 3.21 Maximum velocity constraints in the body frame.
Figure 3.22 Maximum velocity arc sampling.
Figure 3.23 Four possible shapes of the feasible control input set: (a) the ...
Figure 3.24 (a–c) Traces of the UAV flight path and of the target for three ...
Chapter 4
Figure 4.1 Transcription from DNA to RNA and translation from RNA to protein...
Figure 4.2
E. coli
tryptophan measurements normalized by the steady st...
Figure 4.3 Comparison phase angles of (2,2)‐Padé approximation and
.
Figure 4.4 Nonlinear optimization problem.
Figure 4.5 (MATLAB) Model fitting results for Experiments A, B, and C, where...
Figure 4.6 (MATLAB) Model fitting results for Experiments A, B, and C, where...
Figure 4.7 (Python) Sensitivity to the decimal points of the optimal
.
Figure 4.8 The cost function,
, variations with respect to the optimal para...
Figure 4.9
,
‐cyclic adenosine monophosphate (cAMP or cyclic AMP, in short)...
Figure 4.10
Dictyostelium
cAMP oscillation network.
Figure 4.11 The concentration oscillations of the internal cAMP molecules an...
Figure 4.12 Experimental and theoretical distributions of
in Gillespie's d...
Figure 4.13 (MATLAB) The concentration oscillations of the internal cAMP mol...
Figure 4.14 (Python) The concentration oscillations of the internal cAMP mol...
Figure 4.15 (Deterministic ODE model) The concentration change time history ...
Figure 4.16 (Deterministic ODE model) The concentration change time history ...
Figure 4.17 (Stochastic simulation) The concentration change time history wi...
Chapter 5
Figure 5.1 Enzyme–substrate reactions as an input (S) and output (P) system....
Figure 5.2 Feedback control structure for the enzyme–substrate system.
Figure 5.3 Comparison of the frequency responses between the true PI control...
Figure 5.4 Comparison of the step and the impulse responses between the true...
Figure 5.5 Dirac Delta function
and its property with integration.
Figure 5.6 Elementary reactions for enzyme–substrate system with feedback co...
Figure 5.7 Concentration changes of the actual P with the PI controller for ...
Figure 5.8 Concentration changes of the actual P with the PI controller for ...
Figure 5.9 Comparison between
and the true difference.
Figure 5.10 Comparison of the terms in the right‐hand side of
.
Figure 5.11 Responses of the closed‐loop system with the one‐sided subtracti...
Figure 5.12 Comparison of frequency response of
and
.
Figure 5.13 Responses of the closed‐loop system with the degradation,
, and...
Figure 5.14
–
block diagram.
Figure 5.15 Maximum
with respect to
.
Figure 5.16
Upper bound for
.
Figure 5.17 Geometric approach for
lower bound calculation.
Figure 5.18 The square box contacts the unstable region shaded in the figure...
Figure 5.19
lower bound obtained by the geometric approach.
Figure 5.20 Geometric approach for
upper bound calculation.
Figure 5.21
upper bound with the lower bound obtained by the geometric app...
Chapter 6
Figure 6.1 Two modules and a grey node in the network.
Guide
Cover Page
Series Page
Title Page
Copyright
Dedication
Preface
Acknowledgements
Acronyms
Table of Contents
Begin Reading
Appendix A Solutions for Selected Exercises
Index
Wiley End User License Agreement
Pages
ii
iii
iv
v
xiii
xv
xvii
xix
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
295
296
297
298
299
300
301
302
303
304
305
307
308
309
310
311
IEEE Press 445 Hoes Lane Piscataway, NJ 08854
IEEE Press Editorial Board
Sarah Spurgeon,
Editor in Chief
Jón Atli Benediktsson
Andreas Molisch
Diomidis Spinellis
Anjan Bose
Saeid Nahavandi
Ahmet Murat Tekalp
Adam Drobot Peter (Yong) Lian
Jeffrey ReedThomas Robertazzi
Dynamic System Modelling and Analysis with MATLAB and Python
For Control Engineers
Jongrae Kim
University of LeedsLeeds, UK
IEEE Press Series on Control Systems Theory and Applications
Maria Domenica Di Benedetto, Series Editor
Copyright © 2023 by The Institute of Electrical and Electronics Engineers, Inc. All rights reserved.
Published by John Wiley & Sons, Inc., Hoboken, New Jersey.
Published simultaneously in Canada.
No part of this publication may be reprodud, stored in a retrieval system, or transmitted in any form or by any means, electronic, mechanical, photocopying, recording, scanning, or otherwise, except as permitted under Section 107 or 108 of the 1976 United States Copyright Act, without either the prior written permission of the Publisher, or authorization through payment of the appropriate per‐copy fee to the Copyright Clearance Center, Inc., 222 Rosewood Drive, Danvers, MA 01923, (978) 750‐8400, fax (978) 750‐4470, or on the web at www.copyright.com. Requests to the Publisher for permission should be addressed to the Permissions Department, John Wiley & Sons, Inc., 111 River Street, Hoboken, NJ 07030, (201) 748‐6011, fax (201) 748‐6008, or online at http://www.wiley.com/go/permission.
Limit of Liability/Disclaimer of Warranty: While the publisher and author have used their best efforts in preparing this book, they make no representations or warranties with respect to the accuracy or completeness of the contents of this book and specifically disclaim any implied warranties of merchantability or fitness for a particular purpose. No warranty may be created or extended by sales representatives or written sales materials. The advice and strategies contained herein may not be suitable for your situation. You should consult with a professional where appropriate. Further, readers should be aware that websites listed in this work may have changed or disappeared between when this work was written and when it is read. Neither the publisher nor author shall be liable for any loss of profit or any other commercial damages, including but not limited to special, incidental, consequential, or other damages.
For general information on our other products and services or for technical support, please contact our Customer Care Department within the United States at (800) 762‐2974, outside the United States at (317) 572‐3993 or fax (317) 572‐4002.
Wiley also publishes its books in a variety of electronic formats. Some content that appears in print may not be available in electronic formats. For more information about Wiley products, visit our web site at www.wiley.com.
Library of Congress Cataloging‐in‐Publication Data Applied for:
Hardback: 9781119801627
Cover Design: Wiley
Cover Images: © Bocskai Istvan/Shutterstock
To Miyoung
Preface
This book is for control engineers to learn dynamic system modelling and simulation and control design and analysis using MATLAB or Python. The readers are assumed to have the undergraduate final‐year level of knowledge on ordinary differential equations, vector calculus, probability, and basic programming.
We have verified all the MATLAB and Python codes in the book using MATLAB R2021a and Python 3.8 in Spyder, the scientific Python development environment. To reduce the confusion in running a particular program, most of the programs are independent on their own. Organizing programming with multiple files is left as an advanced skill for readers to learn after reading this book.
Leeds, West Yorkshire, England, UK30 November 2021
Jongrae Kim
Acknowledgements
I have learned dynamic modelling and simulation through my undergraduate and post‐graduate education and research projects in the past 30 years. Hence, this book will not be possible without having my teachers, supervisors, and collaborators. I thank Dr Jinho Kim, Professor John L. Crassidis, Professor João P. Hespanhna, Professor Declan G. Bates, Dr Daizhan Cheng, Professor Kwang‐Hyun Cho, Professor Frank Pollick, and Dr Rajeev Krishnadas.
Jongrae Kim
Acronyms
DCM
direction cosine matrix
DNA
deoxyribonucleic acid
EKF
extended Kalman filter
KF
Kalman filter
LHS
left‐hand side
LTI
linear time‐invariant
mRNA
messenger RNA
mRNAP
messenger RNA polymerase
N2L
Netwton's second law of motion
ODE
ordinary differential equation
probability density function
PI
proportional integral
QUEST
quaternion estimation algorithm
RHS
right‐hand side
RNA
ribonucleic acid
About the Companion Website
This book is accompanied by a companion website.
www.wiley.com/go/kim/dynamicmodeling
This website includes:
The solutions for the problems listed in the chapters and the program codes used in Python and MATLAB softwares.
1Introduction
1.1 Scope of the Book
This book is for advanced undergraduate students, post‐graduate students, or engineers to acquire programming skills for dynamic system modelling and analysis using control theory. The readers are assumed to have a basic understanding of computer programming, ordinary differential equations (ODE), vector calculus, and probability.
Most engineering curricula at the undergraduate level include only an elementary‐level programming course in the early of the undergraduate years. Only a handful of self‐motivated engineering students acquire advanced level programming skills mainly from self‐study through tedious time‐consuming practices and trivial mistakes. As modern engineering systems such as aircraft, satellite, automobile, or autonomous robots are implemented through inseparable tight integration of hardware systems and software algorithms, the demand for engineers having fluent skills in dynamic system modelling and algorithm design is increasing. In addition, the emergence of interdisciplinary areas merging the experimental domain with mathematical and computational approaches such as systems biology, synthetic biology, or computational neuroscience further increases the necessity of the engineers who understand dynamics and are capable of computational implementations of dynamic models.
This book aims to fill the gap in learning practical dynamic modelling, simulation, and analysis skills in aerospace engineering, robotics, and biology. Learning programming in the engineering or biology domain requires not only domain knowledge but also a robust conceptual understanding of algorithm design and implementation. It is not, of course, the skills to learn in 14 days or less as many online courses claim. To be confident in dynamic system modelling and analysis takes more than several years of practice and dedication. This book provides the starting point of the long journey for the readers to equip and prepare better for real engineering and scientific problems.
1.2 Motivation Examples
1.2.1 Free‐Falling Object
Newton's second law of motion is given by
where is the ‐th external force in Newtons (N) acting on the object characterized by the mass, , in kg, is the time derivative, is the time in seconds, is the velocity in m/s, and is the momentum of the object. Newton's second law states that the sum of all external forces is equal to the momentum change per unit of time.
Consider a free‐falling object shown in Figure 1.1. There exists only one external force, i.e. the gravitational force acting downwards in the figure. Hence, the left‐hand side of (1.1) is simply given by , where is the gravitational force. Introduce the additional assumption that the object is within the reasonable range from the sea level. With the assumption, the gravitational force, , is known to be proportional to the mass, and the proportional constant is the gravitational acceleration constant, , which is equal to 9.81 m/ in the sea level. Therefore, . Replace the left‐hand side of (1.1), i.e. , by provides
where the downward direction is set to the positive direction, which is the opposite of the usual convention. It highlights that establishing a consistent coordinate system at the beginning of modelling is vital in dynamic system simulation.
Figure 1.1 Free‐falling object.
From the kinematic relationship between the velocity, , and the displacement, , we have
where the origin of is at the initial position of the object, , and the positive direction of is downwards in the figure. The right‐hand side of (1.2) becomes
Finally, the leftmost and the rightmost terms are equal to each other as follows:
and it is expanded as follows:
Using the short notations, , , and , and after rearrangements, the governing equation is given by
For purely educational purposes, assume that the mass change rate is given by
We can identify now that there are three independent time‐varying states, which are the position, , the velocity, , and the mass, . All the other time‐varying states, for example, and , can be expressed using the independent state variables. Define the state variables as follows:
Obtain the time derivative of each state expressed in the state variable as follows:
and this is called the state‐space form.
Let the initial conditions be equal to m, m/s, and kg. Equation (1.5) can be written in a compact form using the matrix–vector notations. Define the state vector, , as follows:
and the corresponding state‐space form is written as
The second‐order differential equation, (1.3), and the first‐order differential equation, (1.4), are combined into the first‐order three‐dimensional vector differential equation, (1.6). Any higher order differential equations can be transformed into the first‐order multi‐dimensional vector differential equation, . Numerical integration methods such as Runge–Kutta integration (Press et al., 2007) solves the first‐order ODE. They can solve any high‐order differential equations by transforming them into the corresponding first‐order multi‐dimensional differential equation.
1.2.1.1 First Program in Matlab
We are ready to solve (1.6) with the initial condition equal to , where the superscript is the transpose of the vector. We solve the differential equation from to seconds using Matlab. Matlab includes many numerical functions and libraries to be used for dynamic simulation and analysis. A numerical integrator is one of the functions already implemented in Matlab. Hence, the only task we have to do for solving the differential equation is to learn how to use the existing functions and libraries in Matlab. The complete programme to solve the free‐falling object problem is given in Program 1.1. Producing Figure 1.2 is left as an exercise in Exercise 1.1.
Figure 1.2 Free‐falling object position, velocity, and mass time histories.
Now, we study the first program line by line. The m‐script starts with the command ‘clear’. The clear command removes all variables in the workspace. In the workspace, there would be some variables defined and used in previous activities. They may have the same names but different meanings and values in the current calculation. For example, the gravitational acceleration ‘grv_const’ in the third line is undefined in the current program and uses a variable of the same name used to analyse objects falling on the moon. A falling object program in the Moon was executed earlier, and ‘grv_const’ is still in the workspace. Without the clear command, the incorrect constant is used in the program producing wrong results. Hence, it is recommended to clear the workspace before starting new calculations. We must be careful, however, that the clear command erases all variables in the workspace. Before the clear command, we check if all values, which might be generated from a long computer simulation, were saved.
From line 3 to line 12, several constants are defined. Based on the equations we have seen earlier, it is tempting to write a code as follows:
These seem to look compact and closer to the equations we derived. It is a bad habit to write a program in this way. The list of problems in the above programming style is as follows:
It defines a variable with a single character, ‘g’, ‘x’, ‘v’, etc. Using a single character variable might cause confusion on the meaning of the variable and lead to using them in wrong places with incorrect interpretations.
Numerical numbers are written without units. There is no indication of units of the numerical values, e.g. 9.81, is it m/
or ft/
?
It uses magic numbers. What do the numbers, 0 and 5, mean in defining ‘t’?
Program 1.1 uses a better style. The initial position is defined using the variable name, ‘init_pos’, whose value is 0.0 and the unit is in metres. Appropriately named variables reduce mistakes and confusion in the program. Program 1.1 indicates the corresponding unit for each numerical value, e.g. the ‘init mass’ value 5.0 is in kg. We understand the meaning of each variable by its name. The texts after ‘%’ are the comments, where we could add various information such as the unit of each numerical value.
In line 13, the built‐in Runge–Kutta integrator, ode45(), is used to integrate the differential equation provided by the function, ‘free_falling_obj’, at the end of the m‐script. Frequently, each function is saved as a separate m‐script. It could also be included in the m‐script for the cases that the functions might be used in the specific m‐script only. To include functions in the m‐script, they must be placed at the end of the m‐script as in this example.
Functions in Matlab begin with the keyword function and close with the keyword end. In line 30, ‘dxdt’ is the return variable of the function and ‘free_falling_obj’ is the function name. The function has three input arguments. A function can have any input argument used by the function. This particular function, ‘free_falling_obj’, is not an ordinary function, however. This is the function to describe the ODE. The function is to be passed into the built‐in integrator, ode45. The first two arguments of the function for ode45 must be time and states, i.e. and in (1.6).
In lines 31–33, the variable ‘state’ is assumed to be a three‐dimensional vector, and each element of the vector corresponds to the states, , , and . In line 35, the return variable ‘dxdt’ is initialized as by the built‐in function zeros(3,1). zeros(m,n) creates the matrix filled in zeros. Lines 36, 37, and 38 define the state‐space form ODE, (1.6).
The function works perfectly well without the initialization line for ‘dxdt’, line 35. However, it is not good programming if line 35 is removed. Without the initialization, ‘dxdt’ in line 36 is a one‐dimensional scalar value. In the next lines, it becomes a two‐dimensional value and a three‐dimensional value. Each line, the size of ‘dxdt’ changes, and this requires the computer to find additional memory to store the additional value. This could increase the total computation time longer and could be noticeably longer if this function is called a million times or more. Hence, it is better to acquire all the required memory ahead as in line 35.
Efficiency vs. development cycle: We strive to create efficient programs, but the prototyping phase requires a fast development cycle.
It is vital to have the habit of being conscious of the efficiency of algorithm implementation. On the other hand, try not to overthink the efficiency of the program. Script languages such as Matlab and Python are for rapid implementation and testing. Hence, it needs a proper balance between optimizing codes and saving the development time.
Now, we are ready to solve the differential equation using the built‐in numerical integrator, ode45. ode45 stands for ODE with Runge–Kutta fourth‐ and fifth‐order methods. Details of the Runge–Kutta integration methods can be found in Press et al. (2007).
Recall, the following line from Program 1.1:
When we use ode45, the input argument starts with @ symbol, which is the function handle. The function handle, @, is used when we pass function A, e.g. ‘free_falling_obj’, to function B, e.g. ode45, where function B would call function A multiple times. With the function handle, we can control or construct the function to be passed with some flexibility. ‘@(time,state)’ explicitly indicates that the function to be passed has two arguments, ‘time’ and ‘state’, and they will be passed between ode45 and ‘free_falling_obj’ function in the specific order, i.e. ‘time’ be the first and ‘state’ be the second argument. This order is required by the integrator, ode45.
With the function handle, we can take some freedom to order the function arguments differently in the function definition of ‘free_falling_obj’. For example, we could write the function as follows:
and the integration part is updated to follow the updated function definition as follows:
The program works the same as the ones before the modifications. Also, we notice that we have an additional input argument, ‘grv_const’. Similarly, we could add more input parameters if they are necessary. As long as the first argument, ‘time’, and the second argument, ‘state’, are indicated in the function handle, the function can have any number of input arguments in any order to pass to the integrator, ode45.
Once the integration is completed, the results return to two output variables, ‘tout’ and ‘xout’. Execute the command, whos, in the Matlab command prompt, the following information is displayed:
The first column shows all variables created including the two output results from the integrator. The second column shows the size of each variable: ‘tout’ is 61 rows and 1 column and ‘xout’ is 61 rows and 3 columns. Hence, each row of ‘xout’ corresponds to the time instance of the corresponding row values of ‘tout’. Why is the number of row 61? This is determined by the integrator automatically to adjust the integration accuracy and computation time. We can assign the number of rows or the number of time steps explicitly, and this is covered in the later chapters. The three columns of ‘xout’ correspond to the state, , , and . The first column of ‘xout’ is for , the second column of ‘xout’ is for , and the last column of ‘xout’ is for .
By executing the following line in the Matlab command prompt, we can print out all values of in the command window:
where ‘:’ indicates all rows. If we want to see the values of from the 11th row to the 15th row, then
Similarly, the time history of is xout(:,2) and the time history of is xout(:,3).
The plot command in Matlab plots the results as follows:
Before plotting each figure, open a new figure window using figure(1), figure(2), and figure(3), respectively. The label for each axis is created using the commands xlabel and ylabel for the horizontal and the vertical axes, respectively, where each axis must indicate what quantity and what units are used.
1.2.1.2 First Program in Python
Program 1.3 solves the free‐falling object differential equation. The program is remarkably similar to the Matlab script in Program 1.1. There are, however, many differences between the two languages.
On lines 4 through 14, the constants are defined with the proper naming and the units indicated in the comments. In Python, comments are placed after #.
