Synthesis Lectures on

Engineering, Science, and Technology

Erik Cuevas · Alberto Luque · Héctor Escobar

Computational
Methods
with MATLAB®
Synthesis Lectures on Engineering, Science,
and Technology
The focus of this series is general topics, and applications about, and for, engineers
and scientists on a wide array of applications, methods and advances. Most titles cover
subjects such as professional development, education, and study skills, as well as basic
introductory undergraduate material and other topics appropriate for a broader and less
technical audience.
Erik Cuevas · Alberto Luque · Héctor Escobar

Computational Methods
with MATLAB®
Erik Cuevas Alberto Luque
University of Guadalajara University of Guadalajara
Guadalajara, Jalisco, Mexico Guadalajara, Jalisco, Mexico

Héctor Escobar
University of Guadalajara
Guadalajara, Jalisco, Mexico

ISSN 2690-0300 ISSN 2690-0327 (electronic)

Synthesis Lectures on Engineering, Science, and Technology
ISBN 978-3-031-40477-1 ISBN 978-3-031-40478-8 (eBook)
https://doi.org/10.1007/978-3-031-40478-8

© The Editor(s) (if applicable) and The Author(s), under exclusive license to Springer Nature
Switzerland AG 2024

This work is subject to copyright. All rights are solely and exclusively licensed by the Publisher, whether the whole
or part of the material is concerned, specifically the rights of translation, reprinting, reuse of illustrations, recitation,
broadcasting, reproduction on microfilms or in any other physical way, and transmission or information storage
and retrieval, electronic adaptation, computer software, or by similar or dissimilar methodology now known or
hereafter developed.
The use of general descriptive names, registered names, trademarks, service marks, etc. in this publication does
not imply, even in the absence of a specific statement, that such names are exempt from the relevant protective
laws and regulations and therefore free for general use.
The publisher, the authors, and the editors are safe to assume that the advice and information in this book are
believed to be true and accurate at the date of publication. Neither the publisher nor the authors or the editors give
a warranty, expressed or implied, with respect to the material contained herein or for any errors or omissions that
may have been made. The publisher remains neutral with regard to jurisdictional claims in published maps and
institutional affiliations.

This Springer imprint is published by the registered company Springer Nature Switzerland AG
The registered company address is: Gewerbestrasse 11, 6330 Cham, Switzerland
Preface

The objective of this book is to present a comprehensive and consistent overview of

numerical methods for undergraduate students and professionals in the engineering field.
The book can also serve as a reference for engineering professionals who have the need
to use MATLAB in their applications. This is because many of the MATLAB codes,
presented after introducing the basic ideas of each algorithm, can be easily modified to
solve similar problems even by those who do not know what goes on inside MATLAB
routines and the algorithms they use.
This book can be successfully used by two types of readers. The first type are those
readers who require a detailed explanation of each algorithm and its potential. These
readers are interested in knowing each method in such a way that they can modify the
structure and make adaptations to solve their problems. The second are those readers
who want to use numerical methods generically as subroutines. Just as most users of a
household appliance need only know how to operate it to perform a task, these users
need only know how to formulate their problems that they want to solve using MATLAB
and how to use the corresponding routines to solve them. It is important to clarify that
detailed knowledge of the numerical method is useful for finding a solution to a particular
engineering problem. However, it is only implied that one-time users of any numerical
method can use this book, as well as readers who want to understand the underlying
principle/equations of each algorithm.
This book focuses primarily on helping readers understand the fundamental mathemat-
ical concepts of numerical methods and practice problem-solving skills using MATLAB.
The methodology of the book is to first teach basic concepts so that readers can correctly
formulate problems mathematically, skipping some tedious and unnecessary checks. Then,
readers can directly implement the codes in MATLAB to solve practical problems. All
algorithms presented in this book are followed by a MATLAB code example so that
students can easily modify the code to solve their own problems. This methodology is
grounded in the belief that the majority of students and professionals, especially those with
backgrounds outside of mathematics, recognize the paramount significance of effectively
employing numerical tools to address their pertinent challenges. This practical application
holds greater importance than considering only extensive demonstrations and proofs.
v
vi Preface

The book uses much more code than formal mathematics. This is because we are
convinced that even readers with an excellent mathematical background have trouble
understanding an approach until they see the algorithm implemented in code. This fact is
because the implemented code removes all ambiguities.
For more than two years, we have tried multiple ways to expose this material to dis-
similar audiences. Along the way, we have counted on the invaluable tolerance of our
students, mainly from CUCEI at the University of Guadalajara in Mexico. Special thanks
are due to our fellow professors of the Centro Universitario de Ciencias Exactas e Inge-
nierías. So many collaborations, help and discussions with colleagues would merit an
additional chapter. To all those mentioned, and especially to those omitted, our testimony
of gratitude.

Guadalajara, Mexico Erik Cuevas

Alberto Luque
Héctor Escobar
Contents

1 MATLAB® Fundamentals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.1.1 MATLAB® Environment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.1.2 MATLAB® Languaje . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.1.3 Trigonometric Operators . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.1.4 Special Operators for Complex Numbers . . . . . . . . . . . . . . . . . . . 6
1.1.5 Relational Operators . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.1.6 Logical Operators . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.1.7 Operators with Matrix and Arrays . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.1.8 Input and Output Operators . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.1.9 Operators for Data Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.1.10 MATLAB® Relationship with Matrices and Vectors . . . . . . . . . . 9
1.1.11 Control Flow Structures in MATLAB® . . . . . . . . . . . . . . . . . . . . . 10
1.1.12 MATLAB® Function Declarations and .m Files . . . . . . . . . . . . . 15
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2 Graphing with MATLAB . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.2 Simple Graphing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.3 Contour of Two-Dimensional Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
2.4 Triangular Mesh and Contours . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
2.5 Surface and Mesh Plotting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
2.6 Matlab as a Graphing Tool Implementing the Mapping Toolbox . . . . . . . 46
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
3 Linear Algebra . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
3.2 Matrices and Vectors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
3.3 Matrix and Vector Operations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
3.4 Systems of Linear Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
3.5 Determinant of a Square Matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63

vii
viii Contents

3.6 Ill-Conditioned Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

3.7 Gaussian Elimination . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
3.8 Gauss–Jordan Elimination . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
3.9 LU Decomposition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
3.10 Eigenvalues of Matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
4 Interpolation and Polynomials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
4.2 Polynomial Expressions in Matlab . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
4.3 Linear Interpolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
4.4 Polynomial Interpolation with Power Series . . . . . . . . . . . . . . . . . . . . . . . . 84
4.5 Lagrange Interpolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
4.6 Polynomial Interpolation Error . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
4.7 Differentiation and Integration of the Lagrange Interpolation
Polynomial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
4.8 Chebyshev and Legendre Interpolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
4.9 Cubic Hermite Interpolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
4.10 Two-Dimensional Interpolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
5 Numerical Differentiation Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
5.2 Differentiation of Interpolation Polynomials . . . . . . . . . . . . . . . . . . . . . . . . 104
5.3 Differences Approximation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
5.4 Differentiation by Taylor Expansion Method . . . . . . . . . . . . . . . . . . . . . . . . 107
5.5 Approximation of Differences in Partial Derivatives . . . . . . . . . . . . . . . . . 108
5.6 Numerical Computation of Higher-Order Derivatives . . . . . . . . . . . . . . . . 109
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
6 Numerical Integration Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
6.1.1 Trapezoidal Rule . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114
6.1.2 Simpson’s Integration Rule . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
6.1.3 Integration by Closed Newton–Cotes Formula . . . . . . . . . . . . . . . 122
6.1.4 Numerical Integration in Two-Dimensional Domains . . . . . . . . . 125
6.1.5 MATLAB® Commands for Numerical Integration . . . . . . . . . . . 128
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128
7 Roots of Nonlinear Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129
7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130
7.2 Graphical Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130
7.3 Bisection Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133
Contents ix

7.4 Newton’s Iteration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137

7.5 Secant Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141
7.6 Method of Successive Substitutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143
7.7 Nonlinear Simultaneous Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150
8 Spline Interpolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151
8.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151
8.2 Linear Splines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156
8.3 Quadratic Splines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 160
8.4 Cubic Spline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163
8.4.1 Cubic Splines with Standard Form Polynomials . . . . . . . . . . . . . 164
8.4.2 Cubic Splines Based on Lagrange Polynomials . . . . . . . . . . . . . . 166
8.5 Matlab Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 171
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 177
9 Differential Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 179
9.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 179
9.2 Symbolic Solution of Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 181
9.3 Euler’s Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 183
9.4 Trapezoidal Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 186
9.5 Runge–Kutta Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 187
9.6 Predictor–Corrector Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 191
9.6.1 Functions for Generic Implementation . . . . . . . . . . . . . . . . . . . . . . 193
9.7 System of Differential Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 196
9.8 Second Degree Differential Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 199
9.9 Differential Equations with Boundary Value . . . . . . . . . . . . . . . . . . . . . . . . 204
9.10 Graphical Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 210
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 212
MATLAB® Fundamentals 1

Abstract

In this chapter, the fundamentals of MATLAB® are discussed. MATLAB® is a power-

ful scientific tool due to its ability to generate graphics and animations and its capacity
to develop diverse complex applications. During this chapter, the environment, lan-
guage, its relationship with matrices and vectors, control flow statemets and function
declarations are discussed.

1.1 Introduction

MATLAB® is software originally written by Cleve Moler, with the goal of providing
matrix software in the late 1970s. Over the years, the platform has developed to the
extent that it is now a standard in engineering and computer science.
MATLAB® is an interpreted, simple, yet fast and powerful language, coupled with its
ability to generate graphics, incorporate animations, and its capacity for the development
of diverse complex applications, and it has become a key tool in the advancement of
scientific research.
During this chapter, we will study the basic fundamentals of the environment and
language itself that will allow us to understand, analyze and create codes that are presented
in later chapters of the text.

© The Author(s), under exclusive license to Springer Nature Switzerland AG 2024 1

E. Cuevas et al., Computational Methods with MATLAB® , Synthesis Lectures
on Engineering, Science, and Technology, https://doi.org/10.1007/978-3-031-40478-8_1
2 1 MATLAB® Fundamentals

1.1.1 MATLAB® Environment

The environment is the set of tools that allow you to work as a user or programmer. It
allows importing and processing data, creating and modifying files, generating graphics
and animations, and developing user applications.
The MATLAB® interface is composed of different windows in which instructions can
be executed or necessary information can be obtained from the workspace. Figure 1.1
shows the initial screen when the program is started:
Depending on the need, instructions can be given in real time in the “command
window” shown in Fig. 1.2
To obtain information on the variables with which you are working, you only need to
consult the “workspace” window shown in Fig. 1.3
.

Fig. 1.1 MATLAB® main window

Fig. 1.2 MATLAB® command window

1.1 Introduction 3

Fig. 1.3 Workspace window

If you want to know the directory in which you are working, just check the “current
folder” window shown in Fig. 1.4.
In the case of creating a script in a new m-file, it will be necessary to click on the “ +
” button shown in Fig. 1.5 and select the script option. This will open the editor window

Fig. 1.4 Current folder

window
4 1 MATLAB® Fundamentals

Fig. 1.5 Button used to create

a new script

Fig. 1.6 Function editor

shown in Fig. 1.6, which will be the area where we can edit and create the programs we
want to carry out.

1.1.2 MATLAB® Languaje

1.1.2.1 Variable Assignment

MATLAB® is a scientific programming language that, unlike many others, does not
require a variable declaration to use variables. The only condition is the correct use of the
following syntax:

Variable_Name = Variable_Value

For example:

A = 32
1.1 Introduction 5

Table 1.1 MATLAB® reserved words

Reserved word Function
Eps Machine Accuracy (smallest number that can be expressed in the computer)
Flops Floating point operations counter
Inf Infinite
Pi Pi number π
Nargin Number of input arguments to the function
Nargout Number of Output arguments in a Function
NaN Not a number (occurs in mathematical indefinitions)
Date Machine date

As MATLAB® is a scientifically oriented language, it is also necessary to take

into account certain reserved words and predefined variables in the language, such as
(Table 1.1).

1.1.2.2 MATLAB® Operators

As you will see below, the fundamental units of MATLAB® are the matrices. Below are
some operators that can be applied to both matrices and scalar quantities in the language
(Table 1.2).
In the case of operating with matrices, when is wanted to perform an operation element
by element of the matrix just add a “.” Before the operator uses, for example, the sum
element to element of matrices can be performed with:

C = a.+ b

Table 1.2 MATLAB®

Operator Command
arithmetical operators
Add +
Subtract −
Multiplication *
Division /
Potentitation ^
Square root sqrt()
Residue of a/b division mod(a,b)
6 1 MATLAB® Fundamentals

Table 1.3 MATLAB®

Operator Command
trigonometric operators
Sine Sin()
Cosine Cos()
Tangent Tan()
Cotangent Cot()
Secant Sec()
Cosecant Csc()

Table 1.4 MATLAB ® Special operators for complex numbers

Operador Comando
Magnitude of a complex number Abs(x)
Angle in radians of the complex number Angle(x)
Generation of a complex number x + yi Complex(x,y)
Conjugate of a complex number Conj(x)
Real part of a complex number Real(x)
Imaginary part of a complex number Imag(x)
Returns 1 if the number is real or 0 if the number is complex Isreal()

1.1.3 Trigonometric Operators

See Table 1.3.

1.1.4 Special Operators for Complex Numbers

See Table 1.4.

1.1.5 Relational Operators

They are operators that are generally used to evaluate a condition for which their output
can be a 1 if the statement is true or a 0 if the statement is false (Table 1.5).
1.1 Introduction 7

Table 1.5 MATLAB®

Operador Símbolo
relational operators
Greater than >
Less than <
Equal to ==
Greater than or equal to > =
Less than or equal to < =
Other than ~ =

Table 1.6 MATLAB® logical operators

Operator Command
Logical AND, gives a 1 as long as both conditions are true &
Logical OR, gives 1 one provided that either of the two conditions is true |

1.1.6 Logical Operators

These operators are typically used to evaluate more than one condition in a single
condition sentence and are (Table 1.6).

1.1.7 Operators with Matrix and Arrays

Since the basic structure of MATLAB® is matrices, some special operators for matrix
arrays are shown below (Table 1.7).

Table 1.7 MATLAB®

Operator Command
operators with matrix and
arrays Dimension of an array Size()
Length of a vector, or number of rows in a matrix Length()
Maximum value of an array Max()
Minimum value of an array Min()
Median of the elements of an array Median()
Average of the elements of an array Mean()
Standard Deviation of the elements in an array Std()
8 1 MATLAB® Fundamentals

Table 1.8 MATLAB® input and output operators

Operator Command
Displays the value of a variable on screen or a message disp(variable)or disp(‘mensaje’)
in character format
Used to display text and data program outputs on screen fprintf(variable)or fprintf(“mensaje”)
or to save them to a file
Receives numeric input from user keypad and stores it in Variable = input(prompt)
a variable
Receives keyboard string input and stores it in a variable String = input(prompt,’s’)

1.1.8 Input and Output Operators

On several occasions, it is necessary to perform an interaction with the user of the code
that is being executed, so it is necessary to have input data acquisition or message display
instructions. In MATLAB® , these are the instructions (Table 1.8).

1.1.9 Operators for Data Analysis

MATLAB® can perform statistical analysis on columnar data sets, including the following
functions (Table 1.9).

Table 1.9 MATLAB®

Operador Comando
operators for data analysis
Correlation coefficient Corrcoef(x)
Covariance matrix Cov(x)
Cumulative Column Sum Cumsum(x)
Cumulative Columns Product Cumprod(x)
Histogram Hist(x)
Maximum in each column Max(x)
Sample interquartile range Iqr(x)
Product of column elements Prod(x)
Sum of column elements Sum(x)
Range of each column Range(x)
Variance of the sample Var(x)
Frequency table of the vector Tabulate(x)
Sort data in ascending order Sort(x)
1.1 Introduction 9

Table 1.10 Commands for fast vector building

Command Function
(i:f) Creates a vector with values from i to f in increments of 1
(i: s:f) Creates a vector with values from i to f in increments or decrements of s
linspace(i,f) Creates a linearly spaced vector between the values i and f with 100 elements
linspace(i,f,n) Creates a linearly spaced vector with n elements
logspace(i,f) Creates a logarithmically spaced vector between the values 10i to 10f with 50
elements
logspace(i,f,n) Creates a logarithmically spaced vector between the values 10i to 10f with n
elements

1.1.10 MATLAB® Relationship with Matrices and Vectors

The fundamental unit in MATLAB® for working with and storing values in its variables
is the matrix, hence the use of vectors, which are a special case of matrices.
To create a vector, you enter the desired values separated by spaces or commas
enclosed in square brackets “[]”. On the other hand, if you want to create a matrix,
simply separate each row vector with the symbol “;”.

1.1.10.1 Access to Elements of Vectors and Matrices

The way to access the elements stored in a vector is the use of an index so that “x(n)”
would be the n-th element of the vector “x”. taking into account that unlike other lan-
guages such as C, in MATLAB® , the indexing of a vector starts counting from the value
1. It is also possible to access a set of elements of the vector in an ordered way using the
operator “:” so that if we enter “x(i:f)”, we will access the set of elements comprising
x(i) to x(f) in steps of one in the indexing.
For a matrix, two values separated by a comma must be added in the indexing so
that “x(f,c)” will return the element in row “f” and column “c” of the matrix we are
working with.

1.1.10.2 Commands for Fast Vector Building

In addition to defining a vector manually, it is sometimes necessary to create vectors with
a larger number of elements, and the following commands can help save time in this task
(Table 1.10).

1.1.10.3 Commands for Fast Array Building

As with vectors, it is often necessary to create special matrices automatically; below are
the commands that will help us to perform this action (Table 1.11).
10 1 MATLAB® Fundamentals

Table 1.11 Commands for fast array building

Comando Función
zeros(n) Creates a square matrix of n x n in which all its values are zero
zeros(n,m) Creates a square matrix of n x m in which all its values are zero
ones(n) Creates a square matrix of n x n in which all its values are one
ones(n,m) Creates a square matrix of n x m in which all its values are one
rand(n) Creates a n x n square matrix of uniformly distributed random numbers with values
between 0 and 1
rand(n,m) Creates a n x m square matrix of uniformly distributed random numbers with values
between 0 and 1
randn(n) Creates a n x n square matrix of uniformly distributed random numbers with values
between 0 and 1
randn(n,m) Creates a n x m square matrix of uniformly distributed random numbers with values
between 0 and 1
eye(n) Creates a n x n square matrix with ones on the diagonal and zeros on the remainder
magic(n) Create a square matrix of n x n so that its rows and columns add up to the same
values

1.1.11 Control Flow Structures in MATLAB®

The flow control structures are those that allow us to execute or repeat certain segments
of code depending on a given condition in MATLAB® .

1.1.11.1 Conditional If-Else Structure

In MATLAB® , the conditional structure has the following syntax:

In this structure, when the condition is met, the first block of instructions is executed;
otherwise, the second block is executed. It is not necessary to add the “else” statement
1.1 Introduction 11

in case you do not want to execute a specific block of instructions if the condition is not
met.

Example: Algorithm 1 Algorithm to detect even and odd numbers in MATLAB®

1.1.11.2 Selective SWITCH Structure

This structure is used to replace multiple if-else conditions in which it is analyzed that a
variable can take different values, and its syntax is:
12 1 MATLAB® Fundamentals

In this case, the variable may have different values, and the code segment that fits the
indicated case will be executed. In case it does not fit in any of the previously indicated
cases, the default instruction defined in the statement “otherwise” will be executed.

Example: Algorithm 2 Algorithm that emulates a selection of arithmetic operations in

MATLAB®
1.1 Introduction 13

1.1.11.3 Cyclic Structure While

This structure is used when you want to repeat the same code segment several times. To
be executed, the structure checks a condition and carries out the instructions as long as
this condition is met. Its syntax is:

while (condition)

Instructions

end
14 1 MATLAB® Fundamentals

One of the most commonly used tricks in this statement is to repeat the loop infinitely
by replacing the condition input with a “1”. It is also possible to interrupt the loop and
exit the loop at any time by using the “break” statement.

Example: Algorithm 3 Algorithm that performs the calculation of the factorial of a number
given by the user in MATLAB®

1.1.11.4 Cyclica Structure For

This structure is an alternative to the “while” statement in which, in an analogous way, a
code segment can be repeated a certain number of times and, in addition, the parameters
of cycle start, cycle end and increment in the cycle interval can be defined from its call.
Depending on the way it is programmed, the increment or decrement value in the cycle
interval can be chosen or not. Its syntax is:

for variable = initial_value : final_value

Instructions

end
1.1 Introduction 15

In this way, the cycle will repeat whenever the variable has values from “start” to
“end”, and the increment interval in the variable will be “1”. Another way to use this
structure is:

for variable = initial_value :increment_value : final_value

Instructions

end

In this alternative way, one will be able to define the increment value that one wishes
the variable to be evaluated to take.

Example: Algorithm 4 Algorithm that displays the desired number of terms of the Fibonacci
sequence in MATLAB®

1.1.12 MATLAB® Function Declarations and .m Files

There are two types of “.m” files in MATLAB® , command or script files and function
statement files. The command file simply has a set of instructions to be executed when
the file name is entered on the MATLAB® command line or the name is included in
16 1 MATLAB® Fundamentals

some other “.m” file. On the other hand, the function files allow us to define functions
analogous to those of MATLAB® in which it is possible to include name, input and output
arguments to make use of them later in future scripts.

1.1.12.1 Function Declaration in MATLAB®

When declaring functions, it is of vital importance to save the “.m” file with the same
name given to the declared function. The syntax for declaring functions in MATLAB® is
the following:

Here, it is possible to notice how the first line starts with the “function” statement
followed by the output arguments [out1,out2,…,outN] enclosed in square brackets
(in case of having only one output argument, it can be declared without enclosing it in
square brackets). Once this is done, the desired name of the function to work with is
placed next to its input arguments (ina1,ina2,…,inaN) enclosed in parentheses.
After that, we write the commands that lead us to obtain the output arguments previously
mentioned, and the output of the function is a vector with the values of the output variables
ordered in the way it was declared at the beginning of the function.

References

1. MATLAB¨/Simulink¨ Essentials: MATLAB¨/Simulink¨ for Engineering Problem ... -Sulaymon

L. Eshkabilov—Google Libros. https://books.google.com.mx/books?hl=es&lr=&id=wVh
8DgAAQBAJ&oi=fnd&pg=PA1&dq=MATLAB+and+simulink+for+engineers&ots=o5U
OvPgbSY&sig=a4hxx4GiThsF1hV24JpMLfa8F-Q&redir_esc=y#v=onepage&q=MATLAB%
20and%20simulink%20for%20engineers&f=false (accessed May 30, 2023).
2. H. Moore, “MATLAB for Engineers,” 2012, Accessed: May 30, 2023. Available: https://ds.
amu.edu.et/xmlui/bitstream/handle/123456789/5248/MATLAB%20for%20Engineers%203rd%
20ed.%20-%20H.%20Moore%20%28Pearson%2C%202011%29%20BBS.pdf?sequence=1&
isAllowed=y.
References 17

3. L. Keviczky, R. Bars, J. Hetthéssy, and C. Bányász, “Introduction to MATLAB,” Advanced Text-

books in Control and Signal Processing, pp. 1–27, 2019, doi: https://doi.org/10.1007/978-981-
10-8321-1_1.
4. D. M. Etter and U. Saddle River Boston Columbus San Francisco New York Indianapolis London
Toronto Sydney Singapore Tokyo Montreal Dubai Madrid Hong Kong Mexico City Munich Paris
Amsterdam Cape Town, “Introduction to MATLAB,” 2002, Accessed: May 30, 2023. [Online].
Available: https://www.pearsonhighered.com/assets/preface/0/1/3/6/0136081231.pdf.
5. E. A. Sobie, “An introduction to MATLAB,” Sci Signal, vol. 4, no. 191, Sep. 2011, doi: https://
doi.org/10.1126/SCISIGNAL.2001984.
6. T. Michalowski, “Applications of MATLAB in Science and Engineering,” 2011, Accessed: May
30, 2023. [Online]. Available: https://books.google.com/books?hl=es&lr=&id=qn-ZDwAAQ
BAJ&oi=fnd&pg=PR11&dq=MATLAB+and+simulink+for+engineers&ots=KqUKAFMF1h&
sig=vGDChrUonRY17vL8kR4wOK78A0k.
7. D. Xue and Y. Chen, “System simulation techniques with MATLAB and Simulink,” 2013,
Accessed: May 30, 2023. [Online]. Available: https://books.google.com/books?hl=es&lr=&id=
6d7iAAAAQBAJ&oi=fnd&pg=PT6&dq=MATLAB+and+simulink+for+engineers&ots=w6P
KsMDIBQ&sig=crq_xM3YFZOY8k7FZaPFjTEt57k.
8. “Applications of MATLAB in Science and Engineering - Google Libros.” https://books.google.
com.mx/books?hl=es&lr=&id=qn-ZDwAAQBAJ&oi=fnd&pg=PR11&dq=MATLAB+and+sim
ulink+for+engineers&ots=KqUKAFMF1h&sig=vGDChrUonRY17vL8kR4wOK78A0k&redir_
esc=y#v=onepage&q=MATLAB%20and%20simulink%20for%20engineers&f=false (accessed
May 30, 2023).
Graphing with MATLAB
2

Abstract

The relationship between calculations and graphics in computing was nearly nonex-
istent in the past, leaving computational users to rely solely on numerical listings to
evaluate research outcomes. However, with significant technological advancements,
researchers, engineers, and professionals from various domains can now easily inte-
grate graphs into their work and research. The importance of incorporating graphs has
been recognized throughout our educational journey, from basic to advanced levels,
and this significance extends to various fields of knowledge, including biology, applied
mathematics, and evolutionary computation. Graphs find direct application in scientific
and technological works, aiding in analyzing research findings, mathematical evalua-
tions, and presentations. Graphs have become an essential tool in research, as they
provide a different perspective and facilitate evaluating research results. MATLAB
stands out among the widely used software for generating high-quality and precise
graphs due to its ability to produce visually striking graphs with just a few simple com-
mands. In Sect. 2.1, a short introduction to the chapter is given. Section 2.2 describes
the essential elements of graphing with MATLAB software. Section 2.3 presents the
implementation of different commands to plot the contour of two-dimensional func-
tions. Likewise, Sect. 2.4 reviews the most important features of MATLAB to create
a triangular mesh. Section 2.5 gives an overview of the combined plots, and Sect. 2.6
The Mapping Toolbox is implemented in creating maps, manipulating satellite images,
and creating a 3D globe.

© The Author(s), under exclusive license to Springer Nature Switzerland AG 2024 19

E. Cuevas et al., Computational Methods with MATLAB® , Synthesis Lectures
on Engineering, Science, and Technology, https://doi.org/10.1007/978-3-031-40478-8_2
20 2 Graphing with MATLAB

2.1 Introduction

Visualizing complicated mathematical functions graphically in up to 3 dimensions and not

having a tool to develop such graphs complicates understanding them [1]. The importance
of implementing graphs has been supported from our basic education to more advanced
stages. This same importance is present in various areas of knowledge, such as biology,
applied mathematics, and evolutionary computation. We can identify a direct application
of graphs by associating them to works carried out in any area of science and technology,
some mathematical analysis developed or simply in some presentation within our work
environment.
In past decades, the relationship between calculations and graphics in computing was
almost nonexistent, forcing computational users to evaluate results from research devel-
oped only in numerical listings. Nevertheless, with the boom that has occurred to date,
technological development has changed radically, allowing researchers, engineers or any-
one who develops a different work to implement graphics in their work and research
easily. Currently, graphs in the research area are of obligatory use since they provide a
different perspective at the moment of the evaluation of research results. One of the most
common software programs in this area is the MATLAB computer program since striking
and precise graphing can be performed with only a few commands.
Understanding and developing mathematical equations through graphs is a simple and
pleasant way to learn mathematics. Therefore, it is recommended that any mathematical
function or the result of an analysis found in this book be graphed to obtain a better
understanding of it. Finally, it is recommended that once the graph analysis is finished,
it should be eliminated in the same way as the implemented variables since MATLAB
could present problems when implementing certain commands. If the problems persist, it
is recommended to restart the program and try to execute the instruction again.

2.2 Simple Graphing

In the MATLAB graphing environment, to make a plot, we use the plot command, which
requires two variables, x and y, as arguments. These variables represent the position of
a set of data, so these variables must constitute an identical arrangement of rows or
columns of the same length. For example, suppose that we want to plot the function
f (x) = cos(x)ex p(−0.4x), 0 ≤ x ≤ 10, which is shown in Fig. 2.1 and plotted with
code 2.1.
2.2 Simple Graphing 21

Fig. 2.1 Graphic of code 2.1 or 2.2

Code 2.1
x=(0:0.05:10);
y=cos(x).*exp(-0.4*x);
plot(x,y)
xlabel('x'); ylabel('y')

Column vectors can also be implemented and established in the plot command
arguments, as shown in code 2.2.

Code 2.2
x=(0:0.05:10)';
y=cos(x).*exp(-0.4*x);
plot(x,y)
xlabel('x'); ylabel('y')
22 2 Graphing with MATLAB

Fig. 2.2 Graphic of code 2.3

The two previous codes provide the same plot presented in Fig. 2.1. The axis labels
are established under the xlabel and ylabel commands, which will be explained in more
detail later. Figure 2.2 shows the graph from code 2.3, constituted by a series of points
connected in a complex plane.

Code 2.3
x=(0:0.05:8*pi);
y=(sin(-x)+i*sin(-2*x)).*exp(-0.05*x)+0.01*x;
plot(real(y), imag(y))

xlabel('Re(y)'); ylabel('Im(y)')

Among the possibilities available at the time of plotting, it is possible to implement

a plot with different symbols without being connected by lines [2]. Implementing this
function is required to give one more argument to the plot function, which is set after the
coordinates. The available types of marks are ‘.’, ‘+’, ‘*’, ‘o’ and ‘x’. These symbols are
dots, plus, stars, circles and x marks, respectively (Fig. 2.3).
2.2 Simple Graphing 23

Fig. 2.3 Graphic of code 2.4 with circle symbols only

Code 2.4
x=(0:0.3:10)';
y=cos(x).*exp(-0.4*x);
plot(x,y,'o')
xlabel('x'); ylabel('y' )

To be able to plot a graph with both lines and symbols, it is only necessary to plot
a second time. This behavior can be implemented in the same line of code in the plot
command, so the code in code 2.4 would look like this: plot(x,y,x,y,‘o’). The default line
type set is continuous, but MATLAB allows plotting on different line types and the option
to specify a color for the line type.
24 2 Graphing with MATLAB

LINE SYMBOL
TYPE OF LINE SYMBOL COLOR

SOLID
- RED r

YELLOW y
DASHE
-- MAGENTA m

DOTTED
:
TURQUOISE c
DASHES- DOTS
-. GREEN g

BLUE b

WHITE w

BLACK k

After the coordinates plot(x,y,‘–’), the desired line type must be added to stipulate
with the command plot command a different line from the default one. In the same way,
when you want to implement a different line color, you should stipulate it in the plot
command followed by the coordinates plot(x,y,‘r’). Combining both modifications, the
symbols and colors are also possible. Specifying it in the plot command is necessary after
the coordinates plot(x,y,‘or’); this line of code will return a plot with red symbols.
An alternative to the plot command is the fplot command, which is implemented as
follows: fplot(‘fun name’,[xmin, xmax]). where fun name refers to the name of the
MATLAB file or the name of the function to be plotted, xmin refers to the minimum
limit of the plot, and xmax refers to the maximum limit of the plot. These limits are
considered the actual minimum and maximum of the function, but this can be adjusted
with the command axis. The command clf removes all existing information in the plot
windows, while the cla command removes the previously plotted curves and redraws
the axes. On the other hand, there are functions in the implicit state, for example, y 3 +
exp(y) = tanh(x). This prevents a correct plotting since it is impossible to express y as a
function of x or x as a function of y. However, it is possible to plot the function thanks
to the contour command. The plotting values, such as the scale marks, the minimum
and maximum coordinates and the values of the coordinates in the scale marks, are set
automatically. However, we can modify the minimum and maximum coordinates and the
frame’s shape thanks to the axis command. Therefore, we obtain a square-shaped graph
with the axis(‘square’) command, as shown in Fig. 2.4.
Another alternative is the axis command; this instruction removes the scale marks and
coordinates axes by specifying the word off, axis(‘off’). T his effect can be reversed at any
time by implementing the counterpart of the previous command, axis(‘on’). The minimum
2.2 Simple Graphing 25

Fig. 2.4 Plot with axis(‘square’) command

and maximum coordinates within the plot can be implemented in a specific way using the
axis command ([xmin, xmax, ymin, ymax]). This code will cut the lines that could go out
of the established limits. This command is implemented after the plot command to modify
the display of the plot as many times as desired. For example, the reader is encouraged
to implement the axis([−1.5, 4, −0.9, 0.9]) command in code 2.4. MATLAB can add a
grid to the graph by implementing the grid on command (as shown in Fig. 2.5). While
the grid off command removes the grid, the implementation of the grid command alone
toggles between displaying or not displaying the grid.

Code 2.5
x=(0:0.1:10);
y=cos(x).*exp(-0.4*x);
plot(x,y,'k')
grid on
xlabel('x'); ylabel('y')
26 2 Graphing with MATLAB

Fig. 2.5 Graphic with command grid on

Code 2.6
x=(0:0.01:pi+1);
y=sin(7*x).*exp(-0.1*x);
polar(x,y,'k')
grid on

The polar command allows graphing a function in polar coordinates; an example is

code 2.6, which graphs Fig. 2.6.
The graphs can also be displayed on a logarithmic and semilogarithmic scale with the
loglog command, as illustrated in code 2.7 and Fig. 2.7.
2.2 Simple Graphing 27

Fig. 2.6 Graphic with polar command

Fig. 2.7 Graphic implementing loglog command

28 2 Graphing with MATLAB

Code 2.7
t=(0.01:0.1:5);
x=exp(t);
y=exp(t.*sinh(t));
loglog(x,y)
grid on
xlabel('x'); ylabel('y')

Code 2.8 results in a semilogarithmic plot with the y-axis on the logarithmic scale.

Code 2.8
t=(.1:.1:3);
semilogy(t, exp(t.*t));
grid on
xlabel('t'); ylabel('t.*t')

In the same way, code 2.9 produces a semilogarithmic plot with the x-axis on the
logarithmic scale.

Code 2.9
t=(.1:.1:3);
semilogx(t, exp(t.*t));
grid on
xlabel('t'); ylabel('t.*t')

All coordinate pairs are written in the command to plot different curves in a single
plot command. Additionally, a different color or mark should be assigned to each curve
for practical purposes. This change is indicated after the pair of coordinates, as shown in
Fig. 2.8, a product of code 2.10.
2.2 Simple Graphing 29

Fig. 2.8 Graphic implementing the plot command with two curves

Code 2.10
x = 0:0.05:5;
y = sin(x);
z = cos(x);
plot(x,y,'k-',x,z,'k*')

We have implemented the plot command to plot two different curves, but sometimes
it is necessary to plot an additional curve on top of another previously plotted curve. This
addition is performed using the hold-on command. Once the command is stipulated in the
code, the plot remains on the screen even if another code block is executed. Therefore,
executing the hold-off command is necessary when the desired procedure is finished. The
above can be visualized in code 2.11, shown graphically in Fig. 2.9.
30 2 Graphing with MATLAB

Fig. 2.9 Graphic implementing hold on and hold off commands

Code 2.11
x = 0:0.05:5;
y = sin(x);
plot(x,y,'k-')
hold on
z = cos(x);
plot(x,z,'kx')
hold off

At the moment of having different curves, thanks to the hold-on command, it is nec-
essary to stipulate the different minimum and maximum of the coordinates in the graph.
This property, as we saw before, is codified with the axis command; if the minimum and
maximum are not stipulated, MATLAB will determine them by default and based on the
first graph shown, which can bring several problems at the moment of visualizing the
second curve.
The hold-on command also has a very important function when preparing a graphic
that takes a long time to be displayed. When commands are implemented to redesign the
2.2 Simple Graphing 31

parameters of the figure to be displayed, such as the perspective angles, the color map or
the color axis, among others, to make these modifications, it is necessary to redraw the
entire figure. To save time, all the modification commands are issued first, and then the
plot command is implemented.
We can also stipulate different titles or texts to the coordinates of the figure using the
xlabel, ylabel and title commands. In addition, placing text inside the graph with the text
or gtext command is possible. The first command requires three arguments for its correct
operation, and the first two values are the coordinates where we want the text to start.
The last argument refers to the text string we want to display in the graph, which can be
a predefined string or a variable. When we want to include different titles or texts in our
graph, we can implement the num2str and int2str commands. These commands allow us
to convert numbers to character strings to display in titles or texts. For example, string =
num2str(pi) allows us to obtain a string containing 3.1416.
On the other hand, string = int2str(pi) will return a string containing only the integer
part of the number, in this case, 3. Once we have the string, we can include it in the title
command as follows: title [‘number =’, string]. This example can be clarified with code
2.12 and Fig. 2.10.

Fig. 2.10 Graphic implementing the commands title, text, and gtext
32 2 Graphing with MATLAB

Code 2.12
x = (0:0.01:1)';k = 1.4;
y = (1+(k-1)/2+x.^2).^(k/(k-1));
plot(x,y)
xlabel('x-axis'), ylabel('y-axis')
title('Title set for the graph')

When it is necessary to specify the coordinates in a more guided way, it would be

appropriate to implement the gtext command since it allows you, once the graphic is fin-
ished, to use a pointer with which it is possible to display the text indicated in the position
where it is clicked. The text command can be implemented in three-dimensional axes.
Otherwise, the gtext command can only be used in two-dimensional axes. We can modify
the font and color by implementing the following code text (2, 0, ‘abcdef’, ‘FontName’,
‘Color’).
The color of the text can vary between red, yellow, green, cyan, blue and magenta,
which can be implemented in the code as r, y, g, c, b, and m, respectively, concerning the
font size change depending on the MATLAB edition. The text command also offers us the
possibility to convert the stipulated letters into Greek letters. This change is implemented
with the argument‘symbol’ followed by the argument ‘FontName’. The size of the axes
can be modified by implementing the set command, and the only drawback is that if the
text is placed in the wrong place, there is no way to delete it. The figure will have to be
redrawn after correcting the code. MATLAB allows us to implement subplots with the
subplot command. Using this command, we can display m by n plots in a single figure,
where m, n and k are integers, m and n refer to an array of plots, and k refers to the
sequential number of the plot. The above is illustrated in Fig. 2.11 and code 2.13.

Code 2.13
x = 0:0.3:30;
subplot(2,2,1),plot(x, sin(x),'k'), title('Subplot 2,2,1')
xlabel('x'), ylabel('sin(x)')
subplot(2,2,2),plot(x,x.*sin(x),'k'), title('Subgraph 2,2,2')
xlabel('x'), ylabel('x.*sin(x)')
subplot(2,2,3),plot(x,x.*sin(x).^2,'k'), title('Subplot 2,2,3')
xlabel('x'), ylabel('x.*sin(x).^2')
subplot(2,2,4),plot(x,x.^2.*sin(x).^2,'k'), title('Subplot 2,2,4')
xlabel('x'), ylabel('x.^2*sin(x).^2')
2.2 Simple Graphing 33

Fig. 2.11 Graphic where subgraphs are implemented

The plot3 command allows us to plot in three dimensions. All the rules we have seen
for the plot command can be implemented with plot3. In Fig. 2.12, we can observe the
spiral movement of a particle from an initial point to an endpoint implemented with the
plot3 command. The viewing angle can be modified with the view command or the axis
command to define the space limits. The command plot3 is implemented in code 2.14.

Code 2.14
t = 0:0.1:20;
r = exp(-0.2*t);
th = pi*t*0.5;
z = t;
x = r.*cos(th);
y = r.*sin(th);
plot3(x,y,z,'k')
hold on
plot3([1,1], [-0.5,0], [0,0],'k')
text(1, -0.7, 0 ,'Start')
n = length(x);
text(x(n), y(n), z(n)+2, 'Final')
xlabel('X'); ylabel('Y'); zlabel('Z');
34 2 Graphing with MATLAB

Fig. 2.12 Graphic implementing the plot3 command

2.3 Contour of Two-Dimensional Functions

In code 2.15, we implement a MATLAB function using the meshgrid command, which
allows us to create two-dimensional arrays of the form (x,y), where x is an array of the
x-coordinates of the grid and y is an array of the y-coordinates of the grid. The arrays
mentioned above are necessary to calculate a two-dimensional z-array (Fig. 2.13).

Code 2.15
xa = -2:0.2:2;
ya = -2:0.2:2;
[x,y] = meshgrid(xa, ya);
z = x.*exp(-x.^2-y.^2);
mesh(x,y,z)
title('This is a 3-D graph of z = x * exp(-x^2 - y^2)')
xlabel('x'); ylabel('y'); zlabel('z');
2.3 Contour of Two-Dimensional Functions 35

Fig. 2.13 Graphic implementing the mesh command

Among the tools offered by MATLAB, we have the contour command, which allows
us to plot the contour of a function contained in a two-dimensional array. The tool is
implemented with the command contour(x, y, z, level). Within the command, the argument
z represents the two-dimensional array of the function, and the x and y values constitute
the x and y coordinates in the array.
Finally, the level parameter refers to the contour levels contained in a vector. The x and
y parameters can also be one-dimensional arrays, so the first index of z rotates or changes
direction relative to the y value. Otherwise, the second index rotates or changes direction
relative to the x value. If the mesh is equispaced, then for convenience, it is possible to
implement only the contour(z) function. It is also possible to replace the level argument
with an integer c, which is the number of levels. These levels are obtained by dividing
the maximum and minimum z values into m − 1 intervals.
36 2 Graphing with MATLAB

Code 2.16
axis('square')
xm=-2:.2:2;
ym=-2:.2:2;
[x,y]=meshgrid(xm, ym);
z=x.*exp(-x.^2-y.^2);
zmax=max(max(z));
zmin=min(min(z));
dz=(zmax-zmin)/10;
level = zmin + 0.5*dz:dz:zmax;
h=contour(x,y,z, level,'edgecolor','k');
clabel(h,'manual')
title('Contour plot made with contour(x,y,z, level)')
xlabel('x'); ylabel('y')

The result of code 2.16 can be seen in Fig. 2.14, in which the function implemented for
its graphing is defined in a two-dimensional function. The contours of the plotted figure
are drawn with the command clavel(h,‘manual’). This function enables the programmer
to specify the position of the numbers with the pointer. The contour levels can be drawn
automatically with clavel(h). This behavior can be identified in Fig. 2.15.
2.3 Contour of Two-Dimensional Functions 37

Fig. 2.14 Graphic implementing contour command with level location

Fig. 2.15 Graphic implementing contour command with automatic levels

38 2 Graphing with MATLAB

2.4 Triangular Mesh and Contours

A triangular mesh has geometrically triangular elements, and its simplest implementation
is in finite element analysis [3]. To implement a triangular mesh, two data sets must be
used, one with the numbering of the elements and the second with the numbering of the
points. These files will be structured as follows:

In addition, a function will be integrated that will implement the data files mentioned
above. This function will generate the triangular mesh, and it is presented as follows:
2.4 Triangular Mesh and Contours 39

Triangular_mesh.m
function triangular_mesh(data_elem, coord_point, scale)
hold off
[n_tr,n]=size(data_elem);
[n_pt,n]=size(coord_point);
nmax=data_elem(1,1);
x=coord_point(:,2);
y=coord_point(:,3)*scale;
element_number=input('Do you want to number the elements? 1
yes/0 no:');
points_number=input('Do you want to number the elements? 1 yes/0
no:');
xmin=min(x); xmax=max(x); x_cen=0.5*(xmin+xmax);
ymin=min(y); ymax=max(y); y_cen=0.5*(ymin+ymax);
Dx=xmax-xmin; Dy=ymax-ymin;
if Dx<Dy, xmin=x_cen-Dy/2; xmax=x_cen+Dy/2; end
if Dx>Dy, ymin=y_cen-Dx/2; ymax=y_cen+Dx/2; end
clf; hold off; clc;
axis([xmin, xmax,ymin, ymax])
xlabel('Triangular mesh plot'); hold on
del_x=0.1; del_y=0.1;
for k=1:n_tr
for l=1:3
p=data_elem(k,l+1);
xx(l)= x(p);
yy(l)=y(p);
end
xx(4)=xx(1);
yy(4)=yy(1);
plot(xx, yy)
x_cen=sum(xx(1:3))/3;
y_cen=sum(yy(1:3))/3;
if element_number == 1
text(x_cen - del_x, y_cen-del_y, int2str(k))
end
end
if points_number == 1
for n=1:n_pt
text(x(n), y(n), ['(',int2str(n),')'])
end
end
axis('off')
40 2 Graphing with MATLAB

A main file described in code 2.17 is created to integrate the function described above
and the data defined in the point and element files (Figs. 2.16 and 2.17).

Fig. 2.16 Triangular mesh highlighting numbered elements

Fig. 2.17 Triangular mesh numbering only the points

2.5 Surface and Mesh Plotting 41

Code 2.17
load elements.m
load points.m
triangular_mesh(elements, points, 1.8)

2.5 Surface and Mesh Plotting

Three-dimensional graphics will be discussed and analyzed in this section. This type of
graphics has been improved over time due to the new versions of MATLAB that have
been released. Previously, we discussed the mesh command. However, we will discuss it
in more detail in the following pages. When referring to a mesh plot of a matrix, it is
necessary to remember that one of the most common applications of three-dimensional
graphics is precisely the plotting of a matrix, considering an m by n matrix x. This plot
can be plotted with code 2.18 and visualized in Fig. 2.18.

Fig. 2.18 Mesh graph of a matrix

42 2 Graphing with MATLAB

Fig. 2.19 Implementation of the view command with different perspectives

Code 2.18
for i=1:8
for j=1:7
x(j,i)=sqrt(i^2+j^2);
end
end
mesh(x,'edgecolor','k')
xlabel('i')
ylabel('j')
zlabel('z')

Within the three-dimensional plot, the limits implemented in the three-dimensional

space of the plot are automatically weighted, but depending on the application to which
you want to orient this can be modified with the axis command [xmin, xmax, ymin,
2.5 Surface and Mesh Plotting 43

ymax, ymax, zmin, zmax], this command structure is logically inherited from the two-
dimensional plots. Likewise, the x-, y-, and z-axis legends are added in the usual way
with xlabel, ylabel, and zlabel, respectively.
A very important aspect of three-dimensional graphics is the perspective the user can
take to visualize the graphic. This perspective can be modified with the command view
([azimuthal, elevation]) or view ([x, y, z]). Within the above structure, azimuthal represents
the azimuthal angle, and the elevation argument represents the elevation angle. It should
be noted that when the azimuthal and elevation values are 0, the user’s perspective is
at the reference perspective angle, which is located along the negative y-axis. From this
perspective, the three-dimensional graph becomes a two-dimensional plane, with the z-
axis vertically and the x-axis extending horizontally to the right.
Modifying the azimuthal parameter allows us to move the perspective position in a
counterclockwise azimuthal degree around the z-axis from the reference angle. To com-
plete the command analysis, the elevation parameter will enable us to raise the perspective
angle elevation degrees from the x–y plane. By not specifying the azimuthal and elevation
values, the default values are −37.5° and 30°, respectively.
With the implementation of the view command by giving it three-dimensional coor-
dinates as parameters such as view([x, y, z]), it is assumed that the user’s perspective
is along the vector [x, y, z] along the origin. view([0, −1, 0]) is equivalent to view([0,
0]) expressed in angles. This implementation can be seen in code 2.19 and visualized
graphically in Fig. 2.19.

Code 2.19
yp=1:5;
xp=1:4;
[x,y]=meshgrid(xp, yp);
z=sqrt(x.^2+y.^2);

%subplot 1
subplot(221)
mesh(x,y,z,'edgecolor','k')
44 2 Graphing with MATLAB

axis([0,5,0,5,5,0,10])
title('Default perspective')
xlabel('X-axis')
ylabel('Y axis')
zlabel('Z-axis')

%subplot 2
subplot(222)
mesh(x,y,z,'edgecolor','k')
axis([0,5,0,5,5,0,10])
title('Perspective [40,25]')
view([40,25])
xlabel('X-axis')
ylabel('Y axis')
zlabel('Z-axis')

%subplot 3
subplot(223)
mesh(x,y,z,'edgecolor','k')
axis([0,5,0,5,5,0,10])
title('Perspective [70,15]')
view([70,15])
xlabel('X-axis')
ylabel('Y axis')
zlabel('Z-axis')

%subplot 4
subplot(224)
mesh(x,y,z,'edgecolor','k')
axis([0,5,0,5,5,0,10])
title('Perspective [0,0]')
view([0,0])
xlabel('X-axis')
ylabel('Y axis')
zlabel('Z axis')zlabel('z')

Another element that needs to be addressed is the contour with mesh. Implementing
the meshc command allows plotting the contour of z in the x and y planes in addition to
the mesh(z) plots. The above is referenced in code 2.20 and Fig. 2.20.
2.5 Surface and Mesh Plotting 45

Fig. 2.20 Implementation of the meshc command

Code 2.20
dth=pi/20;
j=1:21;
i=1:10;
x=log(i);
y=log(j);
[x,y]=meshgrid(x,y);
z=sqrt(0.1*((x-log(5)).^2+(y-log(5)).^2))+1;
meshc(x,y,z,'edgecolor','k')
view([41,27])
xlabel('Axis X')
ylabel('Y-axis')
zlabel('Axis z')

An important utility of the hold on command is found within the implementation of

graphics that take time to be displayed because when the properties are modified, either
of the color map, the color axis or even the same axis, in addition to other parameters that
can be implemented after having made the figure, but if another command that performs
46 2 Graphing with MATLAB

some of the modifications mentioned above is used after plotting the figure, the figure
is redrawn. Therefore, it saves considerable time to indicate the necessary configurations
with their respective commands before plotting and is retained with the hold on.

2.6 Matlab as a Graphing Tool Implementing the Mapping

Toolbox

Within the possibilities offered by MATLAB, we address an embedded feature within

the mapping toolbox [4]. This toolbox will allow us to plot a set of structured data and
visualize them in the form of city maps or world maps, depending on the data we have.
In addition, we can implement satellite images to point out through legends or signals
specific points to locate important cities inside a country or the outstanding colonies of
certain towns. Additionally, these graphing maps will be manipulated to orient them to an
application.
In addition, the mapping toolbox allows the generation of animated globes. This exer-
cise will be discussed in more detail in its respective section. The type of file needed to
plot is the shapefile format. These file types contain geospatial data vectors for geographic
information software and are identified by their ‘.shp’ ending. Shapefiles specify a vector
of spatial features such as points, lines and polygons, which are implemented to represent
lakes, rivers and others. Each element usually has attributes that describe it, such as a
name or its temperature.
This type of file can be obtained from the following web page https://la.mathworks.
com/help/map/finding-vector-geodata.html. However, it is not subject to only being able
to implement such data. In code 2.21, the mapshow() function is implemented, which
displays the coordinates of the x and y vectors as lines.

Code 2.21
streets = shaperead('concord_roads.shp');
figure
mapshow(streets);
xlabel('East in meters')
ylabel('North in meters')

In code 2.21, the shaperead() function is used to read the shapefile and assign it to
a variable named ‘roads’, which will be the argument of the mapshow function. Finally,
the names of the x and y-axes are specified, which can be identified in Fig. 2.21. As in
the two-dimensional or three-dimensional graphics, we can also modify the style of the
2.6 Matlab as a Graphing Tool Implementing the Mapping Toolbox 47

Fig. 2.21 Graphing implementing the mapshow function

lines of our map by specifying the word LineStyle followed by its parameter within the
mapshow function, as we can see in code 2.22 and appreciate graphically in Fig. 2.22.

Code 2.22
figure
mapshow('concord_roads.shp','LineStyle',':');
xlabel('East in meters')
ylabel('North in meters')

Another way to graphically display maps with the help of mapping toolbox is the
inclusion of geospatial images, which allow us to visualize the satellite image of a specific
area and display labels indicating areas of interest. Code 2.23 shows the implementation
of GeoTIFF files, which can store and transfer satellite images, elevation models and
scanned maps. These files have embedded data in tags that can be implemented when
locating an area of interest (Fig. 2.23).
48 2 Graphing with MATLAB

Fig. 2.22 Modifying the line style with the mapshow function

Fig. 2.23 Data plot boston.tif

2.6 Matlab as a Graphing Tool Implementing the Mapping Toolbox 49

Code 2.23
figure
mapshow boston.tif
axis image manual off

S = shaperead('boston_placenames.shp');
surveyFeetPerMeter = unitsratio('sf','meter');
for k = 1:numel(S)
S(k).X = surveyFeetPerMeter * S(k).X;
S(k).Y = surveyFeetPerMeter * S(k).Y;
end

text([S.X], [S.Y], {S.NAME}, 'Color', [0 0 0], 'BackgroundColor',[0.9

0.9 0],'Clipping','on');

A feature of interest is to visualize in three-dimensional form a globe. For this, we will

establish the perspective that it will have. We use some commands such as axesm, which
allows us to create the access map and has a property called ‘globe’; gridm provides the
access map with various features such as line style and color, among others. The spin
function is used to give the globe the ability to rotate. The spin function is shown below.

SPIN.M
for i=90:-5:-270
view(i,23.5);
drawnow
end

The next objective will be to provide the globe with a solid color for the base, then
proceed to plot the globe to finally add the lighting touches and some perspective adjust-
ments, to the procedure described above can be followed with code 2.24 and observe the
result graphically in Fig. 2.24.
50 2 Graphing with MATLAB

Fig. 2.24 Graphing of a globe

Code 2.2
figure
axesm('globe');
gridm('GLineStyle','-','Gcolor',[.7 .8 .9],'Grid','on')

%% Rotate the globe once only

set(gca,'Box','off','Projection','perspective')
spin

%% Adding color to the base and implementing a softer twist

base = zeros(180,360);
baseR = georefcells([-90 90],[0 360],size(base));
copperColor = [0.62 0.38 0.24];
hs = geoshow(base, baseR,'FaceColor',copperColor);
setm(gca,'Galtitude',0.025);
axis vis3d
spin
2.6 Matlab as a Graphing Tool Implementing the Mapping Toolbox 51

%% Show globe
clmo(hs)
load mole
topo = topo/(earthRadius('km')* 20);
hs = meshm(topo, topolegend,size(topo), topo);
demcmap(topo)
set(gcf,'color','black');
axis off;
spin
camlight right
lighting Gouraud;
material ([.7, .9, .8])

The need to be able to visualize graphically in up to 3 dimensions various com-

plicated problems developed in research and not having a tool to establish graphs to
help us understand the proposed solution complicates the fact of being able to under-
stand the immediate impact of this issue. The importance of implementing graphics
has been supported from our basic training to more advanced stages of our research.
This same importance is present in various areas of knowledge, such as biology, applied
mathematics, and evolutionary computation.
We can identify a direct application of graphs by associating them with work done
in any area of science and technology, some mathematical analysis developed or simply
in some presentation within our work environment. Understanding and displaying our
solutions through graphs supports a simple and pleasant way to make known the results
of our research. Therefore, it is recommended that any function resulting from an analysis
found in relevant research be graphed to understand it better. Finally, it is recommended
that once the analysis of a graph is finished, an exhaustive investigation is carried out
to find different ways of graphing to see an innovative and dynamic visualization, which
could present better parameters for understanding in case the techniques presented in this
section do not meet the requirements of the project. We recommend that the reader use
the most known and coined graphing techniques in the literature, such as those seen in
the first sections of this chapter.
In various investigations that may involve a global vision or a broader view of some
important feature, it is possible to implement certain innovative techniques to ensure
an understanding of what you want to explain. Through innovative designs, it is pos-
sible to arouse the user’s curiosity. Moreover, correct graphing allows us to visualize the
problem and the proposed solution more clearly. Unfortunately, the recent pandemic has
taught us that our area of analysis should not only focus on our country, city or colony,
but we must expand the boundaries when addressing a problem. Therefore, the correct
52 2 Graphing with MATLAB

global visualization through MATLAB graphing is an important tool among the program’s
offerings.
When approaching a global problem, have a correctly sampled space graph. For this,
we resort to not-so-common techniques such as those offered by MATLAB toolbox map-
ping. As we developed pages ago, displaying a moving globe makes it possible to attract
the public’s attention to our experiments and consolidate, as described above, the knowl-
edge related to our research. Another way of graphing our planet is through a world map.
This time, to take it a little further, our priority will be to display a world map where
you can identify each of the continents, as well as the most important cities of the most
significant countries. However, the intention is to make the world map described above
by adding an interactive component in which the user can select a part of the map with
the mouse and the closest city to the selected point, and its location in coordinates will
be displayed on the screen.
To perform the above described, the reader will be guided step by step to develop and
understand how to make an interactive world map and manipulate each section to adapt
effectively to any field or area you want to apply. The world map will be displayed on
the screen, and a series of red dots will be observed on it, indicating a city of interest, so
that the user can intuitively move the mouse to that point to receive the information of
the selected area, which will be graphed and displayed as Fig. 2.25.
First, we will build the map with each region and terrain. For this, we will need to
configure the axes of the map object and render the grid elevations; in addition to building
the axes mentioned above, this is done with the following code block, which begins with

Fig. 2.25 Graph of a World Map

2.6 Matlab as a Graphing Tool Implementing the Mapping Toolbox 53

basic instructions for cleaning variables and output console. Then, the axesm function
is implemented to indicate the projection implemented, which will be the Briesemeister,
allowing us to observe a 2d figure with more emphasis on the continents than the oceans.

clc
clear
close all
figure
axesm bries
text(.8, -1.8, 'Briesemeister projection')
framem('FLineWidth',1)

After configuring the above, we can proceed to load the elevations to be displayed
by loading the topo variable and implementing the geoshow command, which has as a
parameter the display shape and texture of the map. Then, we will plot the most significant
cities within the map, which can be identified with a red circular marker. For this, it is
necessary to load and read a shapefile that contains the names and coordinates in latitude
and longitude of the cities mentioned above. Implementing the shaperead command,
as we have done previously, it is possible to obtain such information, which will be
stored in multiple variables for proper implementation. The name of the cities will be
held in “cities”, as well as the latitude and longitude in the variables “lats” and “lons”,
respectively. Then, the geoshow command will be used again, adding as parameters the
latitude and longitude of the cities and specifying the attributes of the markers. As the
last parameter, one of the titles that will contain the graph is stipulated.

cities = shaperead('worldcities', 'UseGeoCoords', true);

lats = extractfield(cities,'Lat');
lons = extractfield(cities,'Lon');
geoshow(lats, lons,...
'DisplayType', 'point',...
'Marker', 'o',...
'MarkerEdgeColor', 'r',...
'MarkerFaceColor', 'r',...
'MarkerSize', 3)
text(-2.8, -1.8,'Major cities of the world')
54 2 Graphing with MATLAB

With the above, we will only have a 2d map; however, as a last resource, we will add
the interactive aspect where the user can select the cities or a region close to it. It will
display the name of the same and its coordinates in latitude and longitude. For this, it is
necessary to stipulate the use of the whole cycle, which will allow iterating infinitely until
a stop criterion is found; this enables the user to identify a single city and discover the
rest. It is important to note which instructions are used more simply than they are given
that at the moment of developing something, it is possible to leave aside instructions that
seem obvious. However, it is important to remember to be as explicit as possible when
developing a program such as this one. Therefore, we stipulate that it is necessary to give
a click with the mouse to a city or red point to obtain information about the same one.
Additionally, it is specified that for the execution of the program, it is necessary to press
the key entry of the keyboard.

The elevation 1-to-1 is loaded and displayed on the grid.

load mole
geoshow(topo, topolegend, 'DisplayType', 'texturemap')
The terrain to be deployed is improved. A map color suitable for the
elevation is obtained.
demcmap(topo)
Map brightness is increased
brighten(.5)
The coordinates of the coasts are loaded.
load coastlines
%We generalize coastlines to 0.25- degrees of tolerance
[rlat, rlon] = reducem(coastlat, coastlon, 0.25);
The coordinates are plotted in brown color.
geoshow(rlat, rlon, 'Color', [.6 .5 .2], 'LineWidth', 1.5)
disp('End of input.')
end

Finally, we obtain the information of the selected city and display it to the user, or
previously coded, obtaining the values of the variables chosen by implementing the find
function. The geoshow function is executed again, as well as its parameters. However,
this time, we stipulate the display of the name of the city by assigning to h2.string that
attribute. We display the coordinates in latitude and longitude, assigning to h3.string the
value of the necessary variables. It is appropriate to note that it is necessary to convert
the numeric value of the variables city.lat and city.lon with the MATLAB native function
num2string. With the above, it is possible to visualize an interactive world map, as shown
in Fig. 2.26.
2.6 Matlab as a Graphing Tool Implementing the Mapping Toolbox 55

Fig. 2.26 Interactive world map, City of Guadalajara

runCitySelectionLoop = true;
if(runCitySelectionLoop)
h1 = text(-2.8, 1.9, 'Click on the red dots to get the city
information. Press enter to exit');
h2 = text(-2.8, 1.7, '');
h3 = text(-2.8, 1.5, 'City coordinates.');
while true
[selected_lat, selected_lon] = inputm(1);
if isempty(selected_lat)
break
end
d = distance(lats, lons, selected_lat, selected_lon);
k = find(d == min(d(:))),1);
city = cities(k);
geoshow(city.Lat, city.Lon, ...
'DisplayType','point',...
'Marker', 'o', ...
'MarkerEdgeColor','k', ...
'MarkerFaceColor','y', ...
'MarkerSize', 3)
h2.String = city.Name;
h3.String = num2str([city.Lat, city.Lon],'%10.2f');
end
56 2 Graphing with MATLAB

References

1. McMahon, D. (2007). MATLAB demystified. New York, NY: McGraw-Hill.

2. Etter, D. M., Kuncicky, D. C., & Hull, D. W. (2002). Introduction to MATLAB. Hoboken, NJ,
USA: Prentice Hall.
3. Menke, W., & Menke, J. (2016). Environmental data analysis with MATLAB. Academic Press.
4. Mapping toolbox. Mapping Toolbox Documentation - MathWorks Latin America (n.d.).
Retrieved June 8, 2022, from https://la.mathworks.com/help/map/index.html.
Linear Algebra
3

Abstract

The study of linear algebra is fundamental for understanding and working effectively
with numerical methods in MATLAB® . This branch of mathematics forms the basis
of MATLAB® ’s usage due to its powerful matrix and vector operations capabilities,
enabling users to leverage the software’s features fully. This chapter will delve into
the essentials of linear algebra, including solving linear equations and manipulating
matrices and vectors. We will explore the proper representation of linear equations
and familiarize ourselves with the basic operations performed using matrix and vector
notation. Additionally, we will review MATLAB® syntax for working with linear equa-
tions in matrix form. We will also address challenging and unsolvable problems and
supplementary topics that enhance the practical application of algebraic concepts. In
Sect. 3.1, a short introduction to the chapter is given. Section 3.2 describes the essential
elements of matrices and vectors. Section 3.3 presents the implementation of different
operations of matrices and vectors. Additionally, Sect. 3.4 reviews the most important
features of systems of linear equations. Section 3.5 gives an overview of the determi-
nant of a square matrix, and in Sect. 3.6, ill-conditioned problems are addressed. In
Sects. 3.7 and 3.8, the implementation of Gauss elimination and Gauss–Jordan elimi-
nation are reviewed, respectively. Section 3.9 LU decomposition is applied, and finally,
in Sect. 3.10, the eigenvalues of matrices are discussed.

3.1 Introduction

One of the most fundamental tools for studying numerical methods and the main basis
for working in MATLAB® is the study of linear algebra. The main basis for working
in MATLAB® is the study of linear algebra, given its capabilities for the operation of

© The Author(s), under exclusive license to Springer Nature Switzerland AG 2024 57

E. Cuevas et al., Computational Methods with MATLAB® , Synthesis Lectures
on Engineering, Science, and Technology, https://doi.org/10.1007/978-3-031-40478-8_3
58 3 Linear Algebra

matrices and vectors, which means that the knowledge of this tool allows us to make
greater use of the study software.
In this chapter, we will study the fundamentals of linear algebra, linear equation solv-
ing, and the resolution of linear equations. We will examine the correct way to express
linear equations and the basic operations in matrix and vector notation. The syntax of
MATLAB® in linear equations with matrix notation will be reviewed. The difficult and
unsolvable problems will be analyzed, as well as additional topics that will help to make
a greater application of algebra [1, 2].

3.2 Matrices and Vectors

Since in MATLAB® , the terms array and matrix are synonymous. It is necessary to clarify
whether an array is used as a mathematical element or a simple variable. In mathematical
notation, matrices are usually enclosed in square brackets and follow certain rules. This
subtopic will review these rules and discuss matrices and vectors in their strict mathemat-
ical sense. In mathematics, a matrix is a rectangular array of numbers enclosed in square
brackets and can be represented as follows:
⎡ ⎤
x1,1 x1,2 . . . x1,n
⎢ ⎥
⎢ x2,1 x2,2 . . . x2,n ⎥
⎢ ⎥ (3.1)
⎢ .. .. .. . ⎥
⎣ . . . .. ⎦
xm,1 xm,2 . . . xm,n

This array is a matrix of dimension m × n. It can be observed how the subscripts m

change in the vertical direction and the subscripts n change in the horizontal direction.
The matrices are represented by a symbol, for example:
⎡ ⎤
x1,1 x1,2 . . . x1,n
⎢ ⎥
⎢ x2,1 x2,2 . . . x2,n ⎥
M =⎢
⎢ .. .. .. .
⎥
⎥ (3.2)
⎣ . . . .. ⎦
xm,1 xm,2 . . . xm,n

Once the matrix M has been determined, it can be represented in a mathematical

equation in terms of M without having to write the whole matrix, abbreviated as follows:
[ ]
M = xi, j (3.3)

On the other hand, a vector is a special case of a matrix with dimensions either m × 1
or 1 × n, thus having a column vector or row vector, respectively, and can be represented
as follows:
3.2 Matrices and Vectors 59

⎡ ⎤
x1,1
⎢ ⎥
⎢ x2,1 ⎥
vc = ⎢
⎢ .. ⎥
⎥
⎣ . ⎦
xm,1
vr = [x1,1 , x1,2 , . . . , x1,n ] (3.4)

Since, in both cases, one of the indices remains constant, it is omitted from the notation,
leaving the vectors represented in (3.4) as follows:
⎡ ⎤
x1
⎢ ⎥
⎢ x2 ⎥
vc = ⎢⎢ .. ⎥
⎥
⎣ . ⎦
xm
vr = [x1 , x2 , . . . , xn ] (3.5)

In MATLAB® , it is possible to define vectors and matrices by separating column values

with spaces “ ” or commas “,” and separating rows using a semicolon “;”.
For example:

c"?"]3"4="6"5="2"4_

Result:

c"?
3"4
6"5
2"4

The same notation is used for the use of vectors. Other special cases of matrices are:
Scalar matrix: Matrix of dimension 1 × 1 being M = [x].
Square matrix: is a matrix of dimensions m = n.
60 3 Linear Algebra

Identity Matrix: is that square matrix for which its elements in the main diagonal have
the value of 1 and the rest of the elements are 0. An example of the case of a matrix of
3 × 3 is as follows:
⎡ ⎤
100
⎢ ⎥
I = ⎣0 1 0⎦ (3.6)
001

The command in MATLAB® to define an identity matrix is I = eye(n), where I

is assigned an identity matrix of dimensions n × n.
Null Matrix: a matrix in which all its elements have zero value. In MATLAB® , null
matrices are defined as N = zeros(m,n), resulting in a null matrix of dimensions
m × n in the case of using only one argument as N = zeros(n), returning a null
matrix of size n × n.

3.3 Matrix and Vector Operations

In this section, we will study the basic operations with matrices and how to perform them
in MATLAB® [3].
Transpose of a Matrix: The transpose of a matrix is obtained by changing the rows of
the matrix by columns of the matrix.
⎡ ⎤
x1,1 . . . x1,n
⎢ . . . ⎥
para M = ⎢
⎣ .. . . .. ⎦
⎥
xm,1 · · · xm,n
⎡ ⎤
x1,1 . . . xm,1
⎢ . . .. ⎥
MT = ⎢
⎣ .
. . . . ⎦
⎥ (3.7)
x1,n · · · xn,m

For example:
[ ]
38
para M =
46
[ ]
34
M =
T
o
86
para M = [ 7 2 ]
3.3 Matrix and Vector Operations 61

[ ]
7
MT = (3.8)
2

In MATLAB® , this operation is performed with the apostrophe operator being, for
example:

O"?"]9.4_=
OV"?"O)=

Assigning to MT the transpose of the matrix M results in the second example shown
in (3.8).
Addition and Subtraction of Matrices: It is possible to add and subtract matrices or
vectors as long as both have the same dimension. Having for:
[ ]
A = ai, j and B = [bi j ] (3.9)

where A and B are matrices of the same dimension. We define addition and subtraction
as:

C = A±B (3.10)

where C = [ci j ] for:

ci j = ai j ± bi j (3.11)

In MATLAB® , it is only necessary to define two matrices of equal dimensions and use
the addition “+” or subtraction “-” operator depending on the case.
Multiplication of Matrices: Having a matrix of dimensions of A of dimensions m 1 × n 1
and a matrix of B dimension matrix m 2 × n 2 these can be multiplied only if m 1 = n 2 .
Otherwise, the matrices will be incompatible under the multiplication criterion. Then,
since matrices A = [ai j ] and B[bi j ] satisfy the above condition, the product of the two
matrices is given by:

C = AB (3.12)

For C = [ci j ] where:

ci j = ai ,1 b1, j + ai,2 b2, j + · · · + ai ,n bn, j (3.13)

62 3 Linear Algebra

In MATLAB® , it is only necessary to define the matrices and use the multiplication
operator “*” between them. The system will throw an error if the compatibility criterion
is not met under multiplication.
The Inverse of a Matrix or Inverse Matrix: The inverse of a matrix M written as M −1
is that matrix such that when calculating:

M M −1 = I (3.14)

where I is the identity matrix mentioned previously in Sect. 3.1. It is worth noting that
to obtain the inverse of a matrix, it must be a square matrix; that is, its dimensions are
m × m.
In MATLAB® , the inverse of a matrix can be calculated using the command
“inv(M)”, where M is the matrix whose inverse is to be known.

3.4 Systems of Linear Equations

In linear algebra, a system of linear equations is a set of linear equations defined on a

commutative ring. Consider a set of m equations with n unknowns given by:

a1,1 x1 + a1,2 x2 + · · · + a1,n xn = s1

a2,1 x1 + a2,2 x2 + · · · + a2n xn = s2
..
.
am,1 x1 + am,2 x2 + · · · + am,n xn = sm (3.15)

where ai , j are known coefficients, xi are the unknowns of the equation and si are
nonhomogeneous terms. The system can be expressed in matrix notation as follows:
First, the system of linear equations is compacted as follows:

Ax = s (3.16)

where A is the known coefficient matrix defined as:

⎡ ⎤
a1,1 . . . a1,n
⎢ . . .. ⎥
A=⎢ ⎣ .
. . . . ⎦
⎥ (3.17)
am,1 · · · am,n

x is the column vector of unknowns given by:

3.5 Determinant of a Square Matrix 63

⎡ ⎤
x1
⎢ . ⎥
x =⎢ ⎥
⎣ .. ⎦ (3.18)
xn

s is the column vector of nonhomogeneous terms:

⎡ ⎤
s1
⎢ . ⎥
s=⎢ ⎣ .. ⎦
⎥ (3.19)
sm

Once the system is compacted and cleared for x, it follows that the solution to the
system of equations is given by:

x = A−1 s (3.20)

This can be achieved in MATLAB® with the command “x = inv(A)* s” by stor-

ing in “x” the vector of solutions to the system of equations y, where “A” and “s” are the
coefficient matrix and the vector of nonhomogeneous terms, respectively. Unfortunately,
systems of linear equations do not always have a numerical solution because the solution
of the system is geometrically located at the point of intersection of the straight lines
defined by the linear equations of the system. This issue leads us to the following cases:
Illustration 3.1 shows the different solution scenarios for a system of linear equations:
A.—The lines defined in the system have only one intersection and therefore have a
unique solution.
B.—The lines defined in the system are superimposed one over the other, causing
an infinite number of solutions. This scenario makes it impossible to give a numerical
answer.
C.—There are two parallel lines in this case, so the system has no intersection. There-
fore, it does not have a solution or three independent equations for two unknowns, which
causes the system to never be satisfied simultaneously, so it does not have a solution
either.

3.5 Determinant of a Square Matrix

One of the quantities of greatest interest is related to the square matrices since it allows
us to know the unique solution in the systems of linear equations. When the value of the
determinant of a matrix is equal to zero, the system of equations that conform to it does
not have a unique solution. In addition to playing an important role in calculating these
eigenvalues is the determinant.
64 3 Linear Algebra

Illustration 3.1 Possible solutions for a system of linear equations

The determinant of a matrix is denoted by |M| depending on the dimension of the

matrix, and it has different ways of being calculated; in the case of matrices of 2 × 2, it
can be calculated as:
|[ ]|
| |
| x1,1 x1,2 |
|M| = | | = x1,1 · x2,2 − x1,2 · x1,2 (3.21)
| x2,1 x2,2 |

For a matrix of 3 × 3:
3.6 Ill-Conditioned Problems 65

|⎡ ⎤|
| x x x |
| 1,1 1,2 1,3 |
|⎢ ⎥|
|M| = |⎣ x2,1 x2,2 x2,3 ⎦|
| |
| x3,1 x3,2 x3,3 |
= x1,1 · x2,2 · x3,3 + x1,3 · x2,1 · x3,2 + x1,2 · x2,3 · x3,1
− x1,3 · x2,2 · x3,1 − x1,2 · x2,1 · x3,3 − x1,1 · x2,3 · x3,2 (3.22)

For higher dimensional matrices, the determinant can be calculated as follows:

Σ
|M| = (±)xi ,1 · x j,2 · xk,3 . . . kr ,n (3.23)
i , j,k,...r

where the summation implies all possible permutations of the first subscript of x, y (±)
is “+” for the case of an even permutation and “−” for the case of an odd permutation.
To calculate the determinant of a matrix in MATLAB® , it is only necessary to use the
command “det(m)”, where “m” is the square matrix from which we want to obtain the
determinant.

3.6 Ill-Conditioned Problems

Several problems involving many linear equations that are solved tend to be inaccurate
values because serious mistakes occur when implementing the rounding technique. When
this happens, it is called an ill-conditioned problem. The small errors that arise when
solving linear equations by rounding can lead to erroneous results when solving an ill-
conditioned problem.

0.12032x + 0.98755y = 2.00555 (3.24)

0.12065x + 0.98775y = 2.01045 (3.25)

Rounding errors are mostly seen in the example illustrated in Eqs. 3.1 and 3.2. When
analyzing the instance, it is evident that the two equations are very close. The solution
is found and denoted by x1 = 14.7403 y y1 = 0.23942. 0.001 is artificially increased
in Eq. 3.2 to make the coefficient error even more evident in its nonhomogeneous term.
Therefore, the resulting system is the following:

0.12032x + 0.98755y = 2.00555 (3.26)

0.12065x + 0.98775y = 2.01145 (3.27)

The solution to the above system will be identified as x2 = 17.9756 and y2 =

−0.15928. Analyzing the above, it is possible to understand that the differences between
66 3 Linear Algebra

the first solution set (x1 , y1 ) and the second solution set (x2 , y2 ) represent a significant
disparity. This disparity is more important when considering the small increase applied
and the alteration obtained. Using such changes at the time of rounding can cause similar
phenomena when solving the system of equations.
The previously mentioned effect becomes increasingly significant as the number of
equations increases. In the following sections, we will analyze two techniques with which
we will be able to validate whether a specific matrix is ill-conditioned or not. The first
technique consists of implementing the condition number; however, it is necessary to
remember certain fundamental aspects, such as rule 2 of a specific matrix, before applying
this technique. This rule is defined below:
┌⎛ ⎞
┃
┃ Σ
┃
||A|| = √⎝ |Ai , j |2 ⎠ (3.28)
i, j

The value of the condition of matrix A is identified by Cond(A) and is defined by

Eq. 3.29, and the value of the condition always satisfies Cond(A) > 1. The value of the
condition usually increases concerning whether the matrix is ill-conditioned, but it should
be noted that this value does not measure the error.
| |
Cond(A) = ||A|| · || A−1 || (3.29)

In MATLAB, it is possible to obtain the condition value by applying the line of code
“cond(A)”, using the above coding to solve the problem of Eq. 3.1. We have “cond
([0.12065, 0.98775; 0.12032, 0.98755])”. To evaluate a poorly conditioned problem more
accurately, it is also important to have a higher degree of accuracy in the calculations. The
inverse and the determinant of an ill-conditioned matrix are not exact. It is possible to
evaluate if the mentioned errors exist by validating certain points, such as validating that
the determinant of A multiplied by the determinant of the inverse of A is different from
one, validating that the determinant of A multiplied by the determinant of the inverse of
( )
A is different from one: det(A)det A−1 / = 1, validating that the inverse of the inverse
−1
of A is different from A: (A−1 ) / = A, validating that multiplying A times the inverse
of A is different from the identity matrix: A A−1 / = I (A), validating that multiplying the
inverse of A by the inverse of the inverse of A is significantly different from the identity
( )−1
matrix of A, even more than the above validation: A−1 A−1 .

Exercise A
Commonly, Hilbert matrices have the problem that they are ill conditioned. Hilbert
matrices can be defined as A = [ai, j ], where

1
ai, j = (3.30)
i + j −1
3.7 Gaussian Elimination 67

With the above, it is necessary to obtain the condition value and the multiplication of
( )
the determinant of A by the determinant of the inverse of A: det(A)det A−1 for Hilbert
matrices from 6 by 6 to 11 by 11. The solution to the above exercise can be obtained
with the following coding:

engct"cnn
hqt"p?8<33
hqt"k?3<p
hqt"l?3<p
C*k.l+?31*k-l/3+=
gpf
gpf

It is possible to observe that the condition value is high during iteration six by ana-
lyzing the obtained results; however, this value will continue to increase faster as the
( )
iterations increase. The values obtained in det(A)det A−1 are susceptible to the round-
ing techniques implemented by the computer. Therefore, it is possible that the results
shown below are not the same:

iteration 6, conditional value of (A) = 1.495106e + 07, det(A) ∗ det(A∧ − 1) = 1.000000e + 00

iteration 7, conditional value of (A) = 4.753674e + 08, det(A) ∗ det(A∧ − 1) = 1.000000e + 00
iteration 8, conditional value of (A) = 1.525758e + 10, det(A) ∗ det(A∧ − 1) = 1.000000e + 00
iteration 9, conditional value of (A) = 4.931534e + 11, det(A) ∗ det(A∧ − 1) = 1.000001e + 00
iteration 10, conditional value of (A) = 1.602503e + 13, det(A) ∗ det(A∧ 1) = 1.000014e + 00
iteration 11, conditional value of (A) = 5.220207e + 14, det(A) ∗ det(A∧ 1) = 1.000323e + 00

3.7 Gaussian Elimination

Gaussian elimination can be understood in two main factors, forward elimination and
backward elimination. We will take as a basis Eq. 3.11, defining that m = n [4].

a1,1 x1 + a1,2 x2 + a1,3 x3 + · · · + a1,n xn = y1

a2,1 x1 + a2,2 x2 + a2,3 x3 + · · · + a2,n xn = y2
...
am,1 x1 + am,2 x2 + am,3 x3 + · · · + am,n xn = ym (3.31)
68 3 Linear Algebra

We will proceed to implement the first factor of Gaussian elimination by multiplying

the first equation by a2,1 /a1,1 and subtracting it from the second equation to eliminate the
first term of the second equation. Likewise, we apply the same methodology with the first
term of the following equations. Subsequently, i > 2 is eliminated by multiplying the first
equation by ai ,1 /a1,1 and subtracting. Implementing the above, we obtain the following
equations:

a1,1 x1 + a1,2 x2 + a1,3 x3 + · · · + a1,n xn = y1

a ' 2,2 x2 + a ' 2,3 x3 + · · · + a ' 2,n xn = y ' 2
...
'
a n,2 x 2 + a ' n,3 x3 + · · · + a ' n,n xn = y ' n (3.32)

where
( )
' ai,1
a i, j = ai , j − ai , j
a1,1
( )
ai ,1
= yi − y1 (3.33)
a1,1

It is important to note that the first equation is not altered; once this is done, we
eliminate all the initial terms of the following equations by applying the methodology
mentioned above. Once the process is finished, the system of equations will have the
following appearance:

a1,1 x1 + a1,2 x2 + a1,3 x3 + · · · + a1,n xn = y1

a ' 2,2 x2 + a ' 2,3 x3 + · · · + a ' 2,n xn = y ' 2
a '' 3,3 x3 + · · · + a '' 3,n xn = y '' 3
...
(n−1)
an,n xn = yn(n−1) (3.34)

The initial terms of Eq. 3.12 are called a pivot, and it is possible to normalize the
equations by applying a methodology that divides each equation by its main coefficient;
however, the Gaussian elimination technique does not support normalization since the
execution time increases significantly and better performance is obtained by avoiding
normalization. The second factor of Gaussian elimination is backward elimination, which
(n−1) (n−1)
obtains the solution. xn = yn /an,n thus obtaining
( )
(n−2) (n−2) (n−2)
xn−1 = yn1 − an−1,n xn /an−1,n−1
...
3.7 Gaussian Elimination 69

( )
Σ
n
x1 = y1 − a1,i xi /a1,1 (3.35)
i=2

This equation concludes the Gaussian elimination methodology. Pivoting is a technique

that consists of permuting the order of the equations, allowing the coefficient of the pivot
to have a greater value than any other that is located below it in the same column and,
therefore, to be eliminated. ai,i to have a higher value than any other coefficient located
below it in the same column and, therefore, to be eliminated. The second pivot can be
implemented before applying the backward elimination factor, as found in Eq. 3.14. The
second | pivot | technique compares whether the second pivot coefficient is the same as the
first. |a ' 2,2 | > |a ' i,2 , |, If the above is not met, the equation permutes the pivot to be the
one with the highest value before starting the elimination, which is implemented at each
pivot until the forward elimination process is completed.

Exercise B
Solve the equations by implementing double precision based on the following system
of equations, composed of four equations and unknowns. It is prudent to note that the
exact answer is unity by applying double precision and the pivoting technique since each
nonhomogeneous term equals the sum of the coefficients in the same row.

1.334 × 10 − 4 4.123 × 10 + 1 7.912 × 10 + 2 −1.544 × 10 + 3 −711.5698662

1.777 2.367 × 10 − 5 2.070 × 10 + 1 −9.035 × 10 + 1 −67.87297633
9.188 0 −1.015 × 10 + 1 1.988 × 10 − 4 −0.961801200
1.002 × 10 + 2 1.4442 × 10 + 4 −7.014 × 10 + 2 5.321 13824.12100

Applying the following coding

engct"cnn
C?]305562g/26"603452g-23"90;342g-24"/
307662g-25=
309992g-22"405892g/27"402922g-23"/
;02572g-23=
;03::2g-22"2"/302372g-23"30;::2g/26=
302242g-24"306664g-6"/902362g-24"
705432g-22_=
{"?uwo*C)+)=
hqtocv"nqpi"g
70 3 Linear Algebra

We obtain the following result.

x=
9.999999999999570e − 01
1.000000000000000e + 00
9.999999999999980e − 01
9.999999999999991e − 01

3.8 Gauss–Jordan Elimination

The Gauss–Jordan elimination methodology is a variant of the Gaussian elimination tech-

nique where the numbers above and below a pivot are eliminated without differentiating
forward elimination in contrast to backward substitution. However, it is necessary to add
the pivoting technique.

Exercise C
Solve the following equation by Gauss–Jordan elimination defining the augmented matrix
where the first three columns are the coefficients, and the last one is the right member.
This augmented matrix is described below:

C?"]/2026"2026"2034"5=
2078"/3078"2054"3=
/2046"3046"/24:"2_=

Implementing the first pivot (exchange of rows 1 and 2) in addition to normalizing the
first row by dividing it by the pivot, we then proceed to eliminate all elements below the
pivot by applying additions or subtractions of a multiple of the first row with the second
and third rows.
We obtain the second pivot by dividing the row by the element’s value [2, 2] to elimi-
nate the top and bottom elements by applying the same methodology. For the third pivot,
we perform the division and eliminate the top elements obtaining the solution set. The
following code performs the process described above:
3.8 Gauss–Jordan Elimination 71

engct"cnn
C?]/2026"2026"2034"5=
2078"/3078"2054"3=
/2046"3046"/204:"2_=
z?]2.2.2_)=
hktuv"rkxqv"ku"korngogpvgf
v?C*4.<+=
C*4.<+?C*3.<+=
C*3.<+?v=
C*3.<+?C*3.<+1C*3.3+=
Tgoqxcn"qh"gngogpvu"dgnqy"vjg"hktuv"rkxqv
C*4.<+?C*4.<+/C*3.<+,C*4.3+1C*3.3+=
C*5.<+?C*5.<+/C*3.<+,C*5.3+1C*3.3+=
ugeqpf"rkxqv"ku"korngogpvgf
v?C*5.<+=
C*5.<+?C*4.<+=
C*4.<+?v=
C*4.<+?C*4.<+1C*4.4+=
C
Gnkokpcvkqp"qh"vjg"gngogpvu"cdqxg"cpf"dgnqy"
vjg"ugeqpf"rkxqv
C*3.<+?C*3.<+/C*4.<+,C*3.4+1C*4.4+=
C*5.<+?C*5.<+/C*4.<+,C*5.4+1C*4.4+=
C
vjktf"tqy"urnkv"korngogpvgf"dgvyggp"kvu"qyp"
rkxqv
C*5.<+?C*5.<+1C*5.5+=
C
'Gngogpv"cdqxg"vjktf"rkxqv"gnkokpcvgf
C*3.<+?C*3.<+/C*5.<+,C*3.5+1C*5.5+=
C*4.<+?C*4.<+/C*5.<+,C*4.5+1C*5.5+=
C
Qdvckpgf"cu"c"ugv"qh"uqnwvkqpu
z*5+?C*5.6+1C*5.5+=
z*4+?*C*4.6+/C*4.5+,z*5++1C*4.4+=
z*3+?*C*3.6+/C*3.4<5+,z*4<5++1C*3.3+=

At the end of the process, we obtain an identity matrix where the last column is the
solutions.
72 3 Linear Algebra

1.000000000000e + 00 0 0 7.000000000000003e + 00
0 100000000000e + 00 0 7.000000000000003e + 00
0 0 1.00000000000e + 00 2.500000000000000e + 01

3.9 LU Decomposition

The LU decomposition methodology is based on modifying matrix A into a multiplication

of two matrices: matrix L, a lower triangular matrix, and matrix U, an upper triangular
matrix [5]. Assuming that A = LU, we can understand that Ax = y can be rewritten as
LU x = y. The above equation can be solved as Ux = z. Finally, we obtain the following
equation:

Lz = y (3.36)

Analyzing the previous equation, we can conclude that obtaining z does not entail
any difficulty, given the triangular nature of L. The resolution of a system of equations
through the LU decomposition technique performs better than Gaussian elimination in
cases where several systems of linear equations with the same matrix of coefficients, but
different nonhomogeneous terms, must be solved.

Exercise D
Solve the linear equation by LU decomposition, Ax = y, where
⎡ ⎤ ⎡ ⎤
1 2 1 8
⎢ ⎥ ⎢ ⎥
A = ⎣ −1 3 −2 ⎦, y = ⎣ 1 ⎦
3 4 −7 10

Y A = LU with
⎡ ⎤ ⎡ ⎤
1 0 0 1 2 1
⎢ ⎥ ⎢ ⎥
L = ⎣ −1 1 0 ⎦ U = ⎣ 0 5 −1 ⎦
3 −0.4 1 0 0 −10.4

First, solve Lz = y for z

⎡ ⎤⎡ ⎤ ⎡ ⎤
1 0 0 z1 8
⎢ ⎥⎢ ⎥ ⎢ ⎥
⎣ −1 1 0 ⎦⎣ z 2 ⎦ = ⎣ 1 ⎦
3 −0.4 1 z3 10

Analyzing the above, we find the solution: z 1 = 8, z 2 = 1 − (−8)=9,z 3 = 10 − 3(8) −

9(−0.4) = −10.4. Subsequently, we operate Ux = z
3.10 Eigenvalues of Matrices 73

⎡ ⎤⎡ ⎤ ⎡ ⎤
1 2 1 x1 8
⎢ ⎥⎢ ⎥ ⎢ ⎥
⎣ 0 5 −1 ⎦⎣ x2 ⎦ = ⎣ 9 ⎦
0 0 −10.4 x3 −10.4

Therefore, the resulting solution is

x3 = −10.4/ − 10.4 = 1

x2 = (9 + x3 )/5 = 2

x1 = 8 − 2x2 − x3 /1 = 3

3.10 Eigenvalues of Matrices

The characteristic polynomial of a matrix A of n × n is defined by the function.

f (λ) = det[A − λI ] (3.37)

The above function is considered a polynomial of order n. Consider a matrix A

[ ]
21
A=
13

Taking Eq. 3.17, we obtain

det[A − λI ] = (2 − λ)(3 − λ) − 1
= λ2 − 5λ + 5 (3.38)

Therefore, Eq. 3.17 is defined for this example:

f (λ) = λ2 − 5λ + 5 (3.39)

Equating f (λ) = 0, and solving this equality, it is possible to obtain the characteristic
values equal to the eigenvalues of matrix A. To obtain the characteristic values’ coeffi-
cients in MATLAB, the command “c = poly(A)” is implemented, where A is a matrix and
c is an array of coefficients. The command “roots(c)” obtains the characteristic values.
However, we can directly obtain the eigenvalues of matrix A with “eig(A)”; it is necessary
to point out again that the answers of “roots(c)” and “eig(A)” are the same.
74 3 Linear Algebra

Exercise E
Obtain the eigenvalues directly with the command “eig”, expand A to its characteristic
polynomial and obtain the roots of the polynomial. Matrix A is defined as
⎡ ⎤
3 12
⎢ ⎥
A = ⎣1 1 4 ⎦
2 42

To solve the above problem, it is necessary to define the matrix A in MATLAB as

follows: “A = [3 1 2; 1 1 4; 2 4 2]”. Since A is defined, we can obtain the eigenvalues
with “eig(A)” as follows:

−2.5967

1.8179

6.7788

In addition, with the command “poly(A)”, it is possible to obtain the coefficients of

the characteristic polynomial

1.0000 −6.0000 −10.0000 32.0000

The above can be interpreted by implementing Eq. 3.17, obtaining

f (λ) = λ3 − 6λ2 − 10λ + 32

Finally, we implement the command “roots(A)” to obtain the roots of the polynomial;
therefore, we store the coefficients of the characteristic polynomial in a variable “c =
poly(A)” and use it with the command roots “roots(c)”. The result is the same as the one
obtained by implementing the command “eig(A)”.

6.7788
−2.5967
1.8179
References 75

References

1. Greub, W. H. (2012). Linear algebra (Vol. 23). Springer Science & Business Media.
2. Wilkinson, J. H., Bauer, F. L., & Reinsch, C. (2013). Linear algebra (Vol. 2). Springer.
3. Ford, W. (2014). Numerical linear algebra with applications: Using MATLAB. Academic Press.
4. Hartfiel, D. J. (2017). Matrix Theory and Applications with MATLAB® . Crc Press.
5. Sewell, G. (2014). Computational methods of linear algebra. World Scientific Publishing Com-
pany.
Interpolation and Polynomials
4

Abstract

In this chapter, we will discuss the polynomials and interpolation topics, which are
very useful tools to describe and evaluate different data sets. During the chapter, the
expressions in MATLAB® to properly declare a polynomial and their operators will
be discussed; on the other hand, the purpose of interpolation is to obtain discrete data
from known data sets and thus to be able to estimate functional values between these
points. During the chapter, how polynomial expressions are managed in MATLAB® ,
the numerical methods to realize linear interpolation, interpolation with power series,
Lagrange interpolation, how to calculate the error in the interpolation methods, the
integration and differentiation of the Lagrange interpolation polynomial, the Cheby-
shev Legendre interpolation polynomial, Hermite interpolation, and two-dimensional
interpolation methods are discussed.

4.1 Introduction

The purpose of interpolation is to obtain discrete data from known data sets and thus to
be able to estimate functional values between these points. In addition, this methodology
is exploited and derived in other numerical methods, for example, numerical methods
of integration and infinite difference approximations that are derived from interpola-
tion polynomials. Therefore, it is of utmost importance to analyze the mathematics of
such polynomials, their accuracy, and their effects on the selection of data points. Since
there are different ways of expressing interpolation polynomials, this chapter will focus
on studying the forms of the power series and Lagrange interpretation, differentiation,
and integration of interpolation polynomials. Polynomials using nonequidistant points

© The Author(s), under exclusive license to Springer Nature Switzerland AG 2024 77

E. Cuevas et al., Computational Methods with MATLAB® , Synthesis Lectures
on Engineering, Science, and Technology, https://doi.org/10.1007/978-3-031-40478-8_4
78 4 Interpolation and Polynomials

along with Chebyshev points, Lagrange interpolation and two-dimensional transfinite

interpolation will be presented.

4.2 Polynomial Expressions in Matlab

A polynomial is an algebraic expression consisting of a finite set of variables and constants

(coefficients) related through various arithmetic operations, which can be expressed in the
following ways:

P(x) = Cn x n + Cn−1 x n−1 + · · · + C1 x + C0 x 0 (4.1)

where n is the order and C I are the coefficients of the polynomial. In grouped form, it can
be expressed as:

P(x) = ((. . . ((Cn x + Cn−1 )x + Cn−2 )x · · · + Cn )x + C0 ) (4.2)

Finally, in its factored form, it is expressed:

P(x) = Cn (x − r1 )(x − r2 ) . . . (x − rn ) (4.3)

where r i are the roots of the polynomial.

Considering a polynomial of order n, it can be said that it has n roots, which can be of
real or complex values. If all the coefficients in the polynomial are real, all its complex
roots will be in conjugate pairs.
In MATLAB® , polynomials are represented through a row vector by adding the
coefficients in order of descending powers. For example:

P(x) = 5x 3 + 3x 2 − 2x + 8 (4.4)

It is expressed as follows:

p = [5 3 -2 8]

and its roots can be calculated using the "roots()" command so that:

roots(p)p
4.2 Polynomial Expressions in Matlab 79

-1.5373 + 0.0000i

0.4686 + 0.9062i

0.4686 - 0.9062i

That being assigned to a variable is stored as a column vector.

Similarly, knowing the value of the roots, the coefficients can be calculated, considering
that the vector generated from the roots is arbitrary without considering any constant
multiplier, considering the following as roots of the polynomial:

r = [4,3,5]

The derived polynomial is obtained using the command:

P = poly(r)

resulting in the row vector:

P = 1 -12 47 -60

which is equivalent to the polynomial:

P(x) = x 3 − 12x 2 + 47x − 60 (4.5)

On the other hand, to evaluate a polynomial, it is only necessary to use the command:

polyval(Px,Val)

where Px is the polynomial in question and Val is the value at which the polynomial
is going to be evaluated. In the case where Val is a vector of desired values, the output
will be a vector with the answers with the same length as the Val vector.
Furthermore, a polynomial of order n is uniquely determined when n + 1 points
are provided, i.e., a polynomial of order n fitted to n + 1 data points, (xi , yi ), i =
1, 2, . . . , n + 1 is unique.
The coefficients of the polynomial are determined using the function:
80 4 Interpolation and Polynomials

polyfit(x,y,n)

whose input arguments are the sampling point x, the value of the function at the sampling
point y, and the order of the polynomial to be fitted n.
Differentiation and integration: to integrate a polynomial given by the expression:

c1 n+1 c2 n cn
Y = yd x = x + x + · · · + x 2 + cn+1 x + cn+2 (4.6)
n+1 2 2
where cn+2 is the constant of integration.
Having the coefficients of Eq. 4.6 arranged in a row vector, the coefficients of Y are
calculated using the algorithm with the function:

Algorithm 4.1.1 Algorithm to obtain the coefficients of polynomial integration in

MATLAB® .

function pi = poly_itg(c)
n = length(c)
pi = [c.*[n:-1:1].^(-1),0].
end

where c is the vector of coefficients of the polynomial resulting in the vector of

coefficients after integration given by:
 
c1 c2
, , . . . , cn+1 (4.7)
n+1 n

Note: the constant of integration cn+2 is not included in the expression.

On the other hand, by deriving Eq. 4.6, we obtain:

y  = nc1 x n−1 + (n − 1)c2 x n−2 + · · · + cn (4.8)

And its coefficients are calculated with the expression:

Der = polyder (c)

where c is the vector of coefficients, and Der is the vector of coefficients given by:

[nc1 (n − 1)c2 , . . . , cn ] (4.9)

4.2 Polynomial Expressions in Matlab 81

Operations with polynomials: To perform operations on two polynomials, it is necessary

to define two polynomials of the form:

pa = a1 x m + a2 x m−1 + · · · + am x + am+1
pb = b1 x n + b2 x n−1 + · · · + bn x + bn+1 (4.10)

To perform a sum, you can use the algorithm with the function:

Algorithm 4.1.2 Algorithm for Polynomial Addition in MATLAB® .

function pr = polyadd(p1,p2)
n1 = length(p1);
n2 = length(p2);
if n1 == n2
pr = p1 +p2;
end
if n1>n2
p = p1 + [zeros(1,n1-n2),p2];
end
if n1<n2
p = [zeros(1,n2-n1),p1]+p2;
end
end

where pr is the resulting polynomial. To perform a subtraction, it is only necessary to

change the sign of the second argument, leaving the expression.

pr = polyadd(pa, -pb)

The product of polynomials of order mxn results in a polynomial of order m + n given

by the equation:

pr = papb = c1 x d + c2 x d−1 + · · · + cd x + cd+1 (4.11)

In MATLAB® , pr is obtained with the expression:

pr = conv(a,b)

On the other hand, the division of the two polynomials satisfies the equation:

pa = pqpb + pr (4.12)
82 4 Interpolation and Polynomials

where pq is the coefficient and pr is the remainder of the division. These polynomials
( pq y pr ) are obtained with the command:

[pq, pr] = deconv(pa, pb)

4.3 Linear Interpolation

One of the bases for a great variety of existing numerical methods, for example, when
deriving the trapezoidal integration method or the use of the gradient as an approximation
of the first derivative of a function, is linear interpolation (Fig. 4.1).
This operation consists of fitting a line to two or more data points and is given by the
equation:

f (b) − f (a)
g(x) = (x − a) + f (a) (4.13)
b−a

where f (a) and f (b) are known values of f (x) at x = a and x = b, respectively.

Fig. 4.1 Example Linear Interpolation, where the points are the data values to be adjusted and the
line is the result of the process
4.3 Linear Interpolation 83

This adjustment presents an error that is expressed as:

e(x) = 0.5(x − a)(x − b) f  (ξ ), a ≤ x ≤ b, a ≤ ξ ≤ b (4.14)

where ξ is a value dependent on x, which lies in an interval point between a and b. Since
this error cannot be evaluated exactly, this equation is difficult to handle. Even so, it can
be dimensioned as follows |e(x)|:
 
|e(x)| ≤ 0.5|(x − a)(x − b)| max  f  (ξ ) (4.15)
a≤x≤b

Noting that the error is a function of x, which is equal to 0 at x = a and x = b,

the maximum error occurs at approximately the midpoint. If it is assumed that f  (x) is
almost constant over the interval, then f  (ξ ) can be approximated as f  (xm ).
Linear interpolation can be performed in MATLAB® with the following function:

interp1(x,y, xi,method)

where ’x’ is the abscissa of the points to be interpolated expressed as a row vector,
’y’ is the ordinate of the points to be interpolated expressed as a row vector, and ’xi’
are the abscissae to construct the interpolation function expressed as a row vector. The
methods can be ’linear’ for a linear interpolation,’spline’ for a cubic spline interpolation,
and ’cubic’ for a cubic interpolation.

Example 4.1: Assume a function y = y(x) function given in tabular form:

x y
0 1.3558
0.25 0.8002
0.5 0.7255
0.75 0.4359
1 0.2289

where y(x) is a decreasing function of x. Find the values of x that satisfy y =

1, 0.9, 0.75, 0.5 y 0.3 using MATLAB® .
84 4 Interpolation and Polynomials

Solution:

Algorithm 4.2.1 Example of the use of interp1() in MATLAB® .

is considered to be a function of
clear all; close all;
%vector generation and
x = [0,0.25,0.5,0.75,1];
y = [1.3558,0.8002,0.7255,0.4359,0.2289];
%vector of values to interpolate
yi = [1,0.9,0.75,0.5,0.3];
%linear interpolation
xi = interp1(y,x,yi,'linear');
%tabular acomode
[yi',xi']

Resulting in:

%result
ans =
1.0000 0.1601
0.9000 0.2051
0.7500 0.4180
0.5000 0.6947
0.3000 0.9141

4.4 Polynomial Interpolation with Power Series

It is possible to express polynomial interpolation in different ways that can be transformed

into each other, such as Lagrange interpolation, Newton forward and backward interpo-
lation, and power series. Regardless of the expression equation, all forms of polynomial
interpolation are equivalent.
4.4 Polynomial Interpolation with Power Series 85

Assuming data of the form n + 1 data of the form:

x1 x2 . . . xn+1
(4.16)
y1 y2 . . . yn+1

where x1 , x2 , . . . , xn+1 are the abscissae of the data points in increasing order. Since the
increment between consecutive points of x is arbitrary, the polynomial that passes through
n + 1 data points can be written as a power series of the form:

p(x) = c1 x n + c2 x n−1 + · · · + cn+1 (4.17)

where ci are coefficients of the polynomial. Doing p(xi ) = yi for each of the n + 1 data
points, we obtain n + 1 linear equations expressed in matrix notation of the form:

Ac = y (4.18)

For:
⎡ ⎤ ⎡ ⎤ ⎡ ⎤
x1n x1n−1 . . . x1 1 c1 y1
⎢ n n−1 ⎥ ⎢ ⎥ ⎢ ⎥
⎢ x2 x2 . . . x2 1 ⎥ ⎢ c2 ⎥ ⎢ y2 ⎥
A=⎢
⎢ .. .. .. ..
⎥, c = ⎢
⎥ ⎢ .. ⎥, y = ⎢
⎥ ⎢ .. ⎥
⎥ (4.19)
⎣ . . . .1 ⎦ ⎣ . ⎦ ⎣ . ⎦
n−1
n
xn+1 xn+1 . . . xn+1 1 cn+1 yn+1

By solving the system of equations given in Eq. 4.18, the coefficients can be deter-
mined. cn These can also be determined with the command polyfit(x,y,n) as previously
explained.

Example 4.3 Given the following data set:

x y
1 3.552
2.6 4.576
4.2 5.232
5.5 2.852

Find the coefficients of the polynomial fitted to the data set by solving for (4.18)
and determine the value of y for x = 3.5y5.2 using the interpolation formula. Plot the
polynomial along with the data point.
86 4 Interpolation and Polynomials

Solution:

Algorithm 4.3.1 Example of polynomial interpolation in MATLAB® .

clear all; close all; hold off

%generation of the x and y vectors
x = [1,2.6,4.2,5.5];
y = [3.552,4.576,5.232,2.852];
%generation of vectors for data clouds
n = length(x)-1;
%determination of coefficients
coefficient = polyfit(x,y,n);

%Vector to interpolate
xi = [3.5,5.2];

%Polynomial interpolation

yi = polyval(coefficient,xi)
%Graphication
xgraf = 1:0.05:6;
ygraf = polyval(coefficient,xgraf);
plot(xgraf,ygraf,x,y,'x')
xlabel('x')
ylabel('p(x):-,data points:x')

Resulting in:

coefficient =

-0.1557 1.1429 -1.8609 4.4257

yi =

5.2356 3.7545

and the graph (Fig. 4.2).

4.5 Lagrange Interpolation 87

Fig. 4.2 Polynomial Interpolation Graph, solution to Example 4.3.1

4.5 Lagrange Interpolation

One of the alternatives to power series interpolation is Lagrange polynomial interpola-

tion. This method has two interesting advantages: first, it eliminates the need to solve
a system of linear equations and allows interpolation when the interpolation data values
are unknown by simply expressing them as symbols. This increases the computational
efficiency of interpolation over power series and reduces susceptibility to rounding
error effects. In addition, it makes it possible to express the polynomial in terms of
undetermined ordinates of the data points.
Consider a polynomial interpolation expression passing through the points:

x1 x2 . . . xn+1
(4.20)
y1 y2 . . . yn+1

To respect the principle of Lagrange’s method, the product of factors is considered:

u(x) = (x − x2 )(x − x3 ) . . . (x − xn+1 ) (4.21)

88 4 Interpolation and Polynomials

Related to the n + 1 data points shown in (4.20). u(x) is a polynomial of order n at x

such that:

0, x2 , x3 , . . . , xn+1
u(x) = (4.22)
 = 0, x1

Dividing u(x) by u(x1 ) obtains:

u(x) (x − x2 )(x − x3 ) . . . (x − xn+1 )

u 1 (x) = = (4.23)
u(x1 ) (x1 − x2 )(x1 − x3 ) . . . (x1 − xn+1 )
Satisfying: u 1 (x1 ) = 1 y u 1 (xi ) = 0 for i = 2, 3, . . . , n + 1. Similarly:

u(x) (x − x1 )(x − x3 ) . . . (x − xn+1 )

u 2 (x) = = (4.24)
u(x2 ) (x2 − x1 )(x2 − x3 ) . . . (x2 − xn+1 )
Satisfying: u 2 (x2 ) = 1 y u 2 (xi ) = 0 ∀i  = 2. Being able to generalize the expression
of the form:
n+1
u(x) x − xj
u i (x) = = (4.25)
u(xi ) xi − x j
j=1, j=i

where u i (x) is a polynomial of coefficients, also known as a shape function.

In addition, multiplying u 1(x) , u 2 (x), . . . , u n+1 (x) by y1 , y2 , . . . , yn+1 and adding the
results, the result is in turn a polynomial of ordern, which is equal toyi ∀x = xi . The
Lagrangian interpolation method of order n can be expressed by the equation:


n+1
i L (x) = u i (x)yi (4.26)
i=1

In MATLAB® , there is no specific function to perform this action, so you have to

create a code that performs the previously seen formulas, as in the following example:

Example 4.4.1 Given the following data set:

x y
0 3.245
0.5 2.789
1 2.350
1.5 4.532
2 5.783
(continued)
4.5 Lagrange Interpolation 89

(continued)
x y
2.5 3.789
3 2.245

Find the coefficients of the polynomial fitted to the data set and determine the value
of y for x = 0.835 y 2.355 using the Lagrange interpolation method. Plot the polynomial
along with the data point.

Solution:

Algorithm 4.4.1 Example of Lagrange interpolation in MATLAB® .

function [coefficients,yi] = Lagrange(x,y,xi)

%polynomial size
n = length(x);
%Symbolic variable for finding the coefficients
syms t;
polynomial = 0;

plot(x,y,'x')%polynomial calculation
for i =1:n
L = 1;
%development of the multiplication of ui in formula 4.25
for j = 1:n
if i~=j
L=L*(t-x(j))/(x(i)-x(j));
end
end
%summatory to find iL in formula 4.26
polynomial = polynomial +L*y(i);

end
coefficients = expand(polynomial);
%Polynomial graph
hold on;
ezplot(polynomial,[0,max(x)])
if nargin == 3
%Polynomial interpolation evaluation
t=xi;
yi = eval(polynomial)
end
end

Finally, the instructions are executed to provide a solution to Example 4.4.1.

90 4 Interpolation and Polynomials

Algorithm 4.4.2 Solution to Example 4.4.1 in MATLAB® .

clear all, close all

x = [0,0.5,1,1.5,2,2.5,3];
y = [3.245,2.789,2.350,4.532,5.783,3.789,2.245];
xi = [0.835,2.355];
[coefficients,yi] = Lagrange(x,y,xi)

Resulting in:

coefficients =

- (2623*t^6)/11250 + (9301*t^5)/2500 - (34051*t^4)/1800 +

(47513*t^3)/1200 - (2949851*t^2)/90000 + (57877*t)/7500 +
649/200

yi =

2.1224 4.6305

and the graph (Fig. 4.3).

4.6 Polynomial Interpolation Error

To analyze the error produced in the interpolations, it is necessary to express it in a

systemic form, i.e., the error of the polynomial interpolation is given by:

e(x) = f (x) − i L (x) = L(x) f n+1 (ξ ) (4.27)

x1 = a ≤ ξ ≤ b 0 = xn+1 (4.28)

where f n+1 is the n + 1-th derivative of f (x), n + 1 is the number of data points and
L(x) is given by:
4.7 Differentiation and Integration of the Lagrange Interpolation Polynomial 91

Fig. 4.3 Lagrangian interpolation plot resulting from Example 4.4.1

(x − x1 )(x − x2 ) . . . (x − xn )(x − xn+1 )

L(x) = (4.29)
(n + 1)!

ξ is dependent on x and is in the interval of a y b. When f (x) is a polynomial of order

n or smaller, its n + 1-disappears, meaning an error equal to zero. On the other hand, if
the error does not disappear, it can be deduced that:
 
|e(x)| ≤ |L(x)| max  f n+1 (ξ ) (4.30)
a≤ξ ≤b

Thus, having a maximum error estimate for any value of x depends on L(x).

4.7 Differentiation and Integration of the Lagrange

Interpolation Polynomial

Since a function can be approximated by an interpolation polynomial, it is also possible

to approximate the derivative and the integral of a function by deriving and integrating
the interpolation polynomial, respectively. This is the basic principle for the derivation of
92 4 Interpolation and Polynomials

numerical methods of differentiation and integration. In this section, we will study how
to evaluate the derivative and integral of the Lagrangian interpolation polynomial.
By deriving Eq. 4.26, one can express the first derivative of the Lagrangian interpola-
tion in the form:


n+1
i L (x) = u i (x)yi (4.31)
i=1

To evaluate the derivative of the function u i given in Eq. 4.25 in factored form, it must
be expressed in the form of a power series. Considering that u 1 is a polynomial of order
n fitted to a set of data points:

x = [x1 , x2 , . . . , xn+1 ]
(4.32)
y = [1, 0, . . . , 0n+1 ]

Similarly, for u 2 :

x = [x1 , x2 , . . . , xn+1 ]
(4.33)
y = [0, 1, . . . , 0n+1 ]

Having in general for u i

x = [x j ]
0, y j (4.34)
y=
1, yi

Therefore, the polynomial u i can be expressed in the form of a power series by fitting
the polynomial of order n to the data.
The calculation of the coefficients of the power series for all i of u i (x) in MATLAB®
can be performed with the following algorithm:

Algorithm 4.6.1 Algorithm for power series shape fitting in MATLAB® .

function matcoef = serpot(x)

np = length(x);
%calculation of coefficients
for j=1:np
y = zeros(1,np);
y(j) = 1;
matcoef(j,:) = polyfit(x,y,np-1);
end
end
4.7 Differentiation and Integration of the Lagrange Interpolation Polynomial 93

where x is the vector of abscissae and matcoef is a matrix in which i represents the
coefficients of the powers for u i (x).
Now, to calculate the first derivative of a Lagrangian interpolation polymony, it is suf-
ficient to convert each row of matcoef into the vector of coefficients of the first derivative
with the polyder function.

Example 4.6.1 Given the following set of points from Example 4.3, find the vector of
coefficients of the first derivative at the abscissa points.

Solution:

Algorithm 4.6.1 Algorithm for calculating the evaluation of the first derivative at given
points in MATLAB® .

clear all
close all
Generation of the x and y vectors
x = [1,2.6,4.2,5.5];
y = [3.552,4.576,5.232,2.852];
%Vector to interpolate
xi = [3.5,5.2];
n = length(x);
matcoef = serpot(x);
Calculation of the coefficients of the derivative
for i =1:n
cd(i,:) = polyder(matcoef(i,::));
end
%evaluation of the coefficients of the derivative

for inp =1 : length(xi)

for i = 1: n
vd(i)= polyval(cd(i,:),xi(inp));
end
yi(inp) = vd*y';
end
yi
94 4 Interpolation and Polynomials

resulting in:

yi =

0.4159 -2.6084

where matcoef(i,:) is the row of the power coefficients of u i (x), cd(i,:) the
coefficients of powers of xi (x). In turn, yi is the value of the first derivative of the
Lagrange interpolation polynomial evaluated for xi.

4.8 Chebyshev and Legendre Interpolation

One of the usual problems when approximating a function by means of interpolation is

to reduce the maximum error. One of the common measures when performing Lagrange
interpolation is to make use of equispaced points, minimizing the maximum error at the
midpoints of the function. The problem is that the error increases at the extremes of the
function.
An optimal solution is to use the points determined by a Chebyshev polynomial since
the distribution of L(x) is more uniform. Another advantage in this respect is that the
error in the data is not amplified at the extremes, unlike in the case of equispaced points.
In addition, another similar solution is the Legendre point distribution, which is mainly
used for numerical integration.

Polynomials and Chebyshev Points:

The Chebyshev polynomials are given by the expression:

Ch n (x) = 2xCh n−1 (x) − Ch n−2 (x) (4.35)

In MATLAB® , the coefficients of the Chebyshev polynomial in the series can be

calculated through the following algorithm:
4.8 Chebyshev and Legendre Interpolation 95

Algorithm 4.7.1 Algorithm for obtaining coefficients in the power series of the Chebyshev
polynomial in MATLAB® .

function coef = Chebyshev(order)

pbb = [1];
if order == 0
coef = pbb;
return;
end
pb = [1 -1];
if order == 1
coef = pb;
return;
end
for i = 2:order
coef = 2*[pb,0] - [ 0,0,pbb];
pbb = pb;
pb = coef;
end
end

where coef is a row vector with the coefficients of the polynomial. The order is the
order of the desired polynomial.
The roots of the polynomial can be calculated in ascending order using the function:

sort(roots(Chebyshev(order)))

This polynomial of order n can also be expressed as:

 
Ch n (x) = cos ncos−1 (x) , −1 ≤ x ≤ 1 (4.36)

where Ch n (x) has n roots in an interval of [−1, 1] and can be calculated with the
following expression:
 
n + 0.5 − i
xi = cos π , i = 1, 2, . . . , n (4.37)
n

If the interpolation interval is [a, b], the roots given by Eq. 4.37 fit the interval given
by:
   
1 n + 0.5 − i
xi = (b − a)cos π + a + b , i = 1, 2, . . . , n (4.38)
2 n
96 4 Interpolation and Polynomials

Polynomials and Legendre Points:

It is a polynomial used only for interpolation purposes, especially in the Gauss–Legendre
numerical integration method, based on the integration of the interpolation polynomial
using the Legendre points.
The Legendre polynomial is given by the expression.
1
Legn (x) = [(2n − 1)x Legn−1 (x) − (n − 1)Legn−2 (x)] (4.39)
n
Thus, the algorithm for calculating the coefficients of the polynomial can be as follows:

Algorithm 4.7.2 Algorithm for obtaining coefficients in the power series of the Legendre
polynomial in MATLAB® 0.8.

function coef = Legendre(order)

pbb = [1];
if order == 0
coef = pbb;
return;
end
pb = [1,0];
if order == 1
coef = pb;
return;
end
for i = 2:order
coef = ((2*i-1)*[pb,0] - (i-1)*[0,0,pbb])/i;
pbb = pb;
pb = coef;
end
end

4.9 Cubic Hermite Interpolation

It is possible to fit polynomials to derivatives in addition to functional values. These poly-

nomials are called Hermite Interpolation Polynomials or Oscillating Polynomials. Consider
a cubic Hermite polynomial:

H er ( p) = c1 p 3 + c2 p 2 + c3 p + c4 (4.40)
4.9 Cubic Hermite Interpolation 97

Fitted to two functional values and two derivatives. The interval between the points p1
y p2 and assuming that the functional values y of the first derivatives are specified at the
points 1 y 2. This results in the following four equations:

H er ( p1 ) = c1 p13 + c2 p12 + c3 p1 + c4 = f 1
H er  ( p1 ) = 3c1 p12 + c2 p1 + c3 = f 1
H er ( p2 ) = c1 p23 + c2 p22 + c3 p2 + c4 = f 2
H er  ( p2 ) = 3c1 p22 + c2 p2 + c3 = f 2 (4.41)

Equation 4.41 is a system of four linear equations with coefficients to be determined.

The implementation of the cubic Hermite interpolation form is simplified by approx-
imating the boundary conditions of the derivative through a finite difference formula
because the Hermite interpolation method can be approximated using the Lagrangian
method. This approximation strategy has a negligible margin of error. The approximation
of the boundary conditions for 4.41 are given by:

H er ( p1 ) = f 1
H er ( p1 + ζ ) = f 1 + ζ f 1
(4.42)
H er ( p2 ) = f 2
H er  ( p2 − ζ ) = f 2 + ζ f 2

where ζ is a random parameter satisfying ζ  x2 − x1 . Therefore, it is possible to deter-

mine a third-degree polynomial fitted to the conditions of (4.41) through the Lagrange
interpolation method.

Example 4.8.1 Determine a curve passing through the point A and the point B passing
through the coordinates x − y that satisfies the following conditions:
dy
A : x = 1, y = 1 =0
dx
dy
A : x = 2, y = 2 =0
dx

Solution:
Using the criteria seen in (4.41), the data points are established as follows:

p = 0 : x = 1, y = 1,
p = ζ : x = 1 + aζ, y = 0
p = 1 − bζ : x = 4, y = 2 − bζ
98 4 Interpolation and Polynomials

p = 1 : x = 4, y = 2 (4.43)

Adhering to the criterion ζ  x2 − x1 the value of ζ is taken as a very small positive

quantity. However, if the quantity is too small, rounding errors will occur when deter-
mining the coefficients, so the value of ζ = 0.01 is selected, and the values a and b are
arbitrary parameters.
The coefficients of the cubic polynomials are determined by the following algorithm:

Algorithm 4.8.1 Algorithm for obtaining Cubic Hermite coefficients in MATLAB® .

clear
hold off

z = 0.01; %zeta
a = 3; %parameter a
b = 3; %parameter b
%%equations 4.40
p(1) = 0;
x(1) = 1;
y (1) = 1;
p(2) = z;
x(2) = 1+z*a;
y (2) = 1;
p(3) = 1-z;
x(3) = 4;
and (3) = 2-z*b;
p(4) = 1;
x(4) = 4;
y(4) = 2;

c = polyfit(p,x,length(p)-1);
d = polyfit(p,y,length(p)-1);
pp = 0:0.05:1;
px = polyval(c,pp);
py = polyval(d,pp);
plot(px,py)
fprintf("The result is: \n");
c
d
4.10 Two-Dimensional Interpolation 99

Fig. 4.4 Hermite curve produced by Algorithm 4.8.1 in MATLAB®

Throwing:

The result is:

c =
-3.0921 3.1231 2.9691 1.0000
d =
1.0307 -0.0309 0.0002 1.0000

and the curve (Fig. 4.4).

4.10 Two-Dimensional Interpolation

Bilinear Interpolation
An acceptable way to interpolate data from a two-dimensional function table is to use
linear interpolation twice (Fig. 4.5).
100 4 Interpolation and Polynomials

Fig. 4.5 Bilinear interpolation

in a two-dimensional domain

A two-dimensional function table is an array of values in a rectangular grid. f i , j =

f (xi , y j ) in a rectangular grid (xi , y j ). Assume that it is necessary to estimate the function
value at a given point in a rectangular domain defined by xi−1 ≤ x ≤ x_i and y j−1 ≤
y ≤ y j , as shown in Illustration 5. Using a linear interpolation in the direction of y, it is
obtained that the values in E and F are:
yj − y y − y j−1
E= f i−1, j−1 + f i−1, j
y j − y j−1 y j − y j−1
yj − y y − y j−1
F= f i−1, j−1 + fi , j (4.44)
y j − y j−1 y j − y j−1

Therefore, the linear interpolation of E y F results in:

xi − x x − xi−1
p(x, y) = E+ F (4.45)
xi − xi−1 xi − xi−1
Combining both steps in one equation, bilinear interpolation is written:
1
p(x, y) =
(xi − xi−1 )(y j − y j−1
    
× (xi − x) y j − y f i−1, j−1 + (xi − x) y − y j−1 f i−1, j
    
+(x − xi−1 ) y j − y f i, j−1 + (x − xi−1 ) y − y j − 1 f i, j (4.46)
References 101

Double Lagrange Interpolation

It consists of applying the Lagrange interpolation method in two dimensions, thus being
able to apply the interpolation formula for all the data points of the table. Assume that
the table of the function has N rows and M columns. The coordinates of the points
are denoted by (xm , yn ) and the function values by F(xm , yn ). Thus, having a double
Lagrange Interpolation given by:


M 
N
L 2D (x, y) = φm (x)n (y)F(xm , yn ) (4.47)
m=1 n=1

where φm and n are given by:

M
x − xk
φm (x) =
xm − xk
(k=1,k=m)
N
y − yk
ψn (y) = (4.48)
yn − yk
(k=1,k=m)

It is observed that φm (x)n (y) takes a value of zero at all data points except for
(xm , yn ). For M = N = 2 Eq. 4.47 reduces to the bilinear interpolation Eq. 4.46.
Chebyshev or Lobatto points can be used to select xm y yn .

References

1. G. Phillips, “Interpolation and approximation by polynomials,” 2003, Accessed: May 30, 2023.
[Online]. Available: https://books.google.com/books?hl=es&lr=&id=87vciTxMcF8C&oi=fnd&
pg=PR7&dq=interpolation+and+polynomial+approximation&ots=bG-Nc80Cgw&sig=v-VA7
aZa7z89V_c_PQnmlIr7mSM
2. R. Burden and J. Faires, “Interpolation & Polynomial Approximation Cubic Spline Interpola-
tion III,” 2016, Accessed: May 30, 2023. [Online]. Available: http://www.personal.psu.edu/jjb23/
web/htmls/sl455SP12/ch3/CH03_5C.pdf
3. D. Lee, R. Cheung, … W. L.-I. T. on, and undefined 2008, “Hardware implementation tradeoffs
of polynomial approximations and interpolations,” ieeexplore.ieee.org, Accessed: May 30, 2023.
[Online]. Available: https://ieeexplore.ieee.org/abstract/document/4384474/
4. M. P.-T. C. Journal and undefined 1967, “On the maximum errors of polynomial approximations
defined by interpolation and by least squares criteria,” academic.oup.com, Accessed: May 30,
2023. [Online]. Available: https://academic.oup.com/comjnl/article abstract/9/4/404/390372
Numerical Differentiation Methods
5

Abstract

During this chapter, the topic of numerical differentiation methods is discussed.

Obtaining the derivative of a function is undoubtedly one of the most widely used
mathematical tools to solve various engineering problems, and calculating it with a
computer can be an interesting problem. During the rest of the chapter, the themes
of differentiation of interpolation polynomials, differences approximations, and dif-
ferentiation by Taylor expansion methods are treated; furthermore, the topics of
approximation of differences in partial derivatives and high-order derivatives are
discussed.

5.1 Introduction

One of the most widely used mathematical tools to solve various engineering problems
is undoubtedly obtaining the derivative of a function. A common problem when working
with functions through a computer is that the data set of the functions is obtained in a
discrete way.
Numerical methods of differentiation or difference approximation are useful for eval-
uating the value of the derivative of a function using discrete data points. Assuming that
the values of the discrete data function are known, then it is possible to express the
approximation of the function by an interpolation polynomial.

© The Author(s), under exclusive license to Springer Nature Switzerland AG 2024 103
E. Cuevas et al., Computational Methods with MATLAB® , Synthesis Lectures
on Engineering, Science, and Technology, https://doi.org/10.1007/978-3-031-40478-8_5
104 5 Numerical Differentiation Methods

5.2 Differentiation of Interpolation Polynomials

Considering a data set:

(xi , yi ), i = 1, 2, . . . , n + 1 (5.1)

It is possible to fit it to an interpolation polynomial through the power series

interpolation method, as explained in the previous chapter, expressed as follows:

p(x) = c1 x n + c2 x n−1 + · · · + cn x + xn+1 (5.2)

To find the first derivative of p(x) for x = 0, we differentiate p(x) and take x = 0.
Giving as a result p(0) = cn . In turn, the second derivative is p '' (x) = 2Cn−1 , and the
general method for the i-th derivative is:

p (i ) (0) = cn+1−i i!, i = 1, 2, . . . , n (5.3)

Eliminating ci from Eq. 5.2 by Eq. 5.3 and inverting the order of the terms, we obtain:

p''(0) 2 p'''(0) 3 p (n) (0) n

p(x) = p(0) + p'(0)x+ x + x + ··· + x (5.4)
2 3! n!
This is the Taylor series expansion form of Eq. 5.1 around x = 0, also known as
Mclaurin’s series.
These relationships can be applied to differential calculus for any value of x. Assume
that we wish to obtain the derivative of p(x) for x = a. Then, using the coordinate
transformation z = x − a, it is possible to express p(x) in terms of z as follows:

p(x) = p(z) = b1 z n + b2 z n−1 + · · · + bn z + bn+1 (5.5)

Therefore, the derivatives p (i) (a):

[( ) [
d i
p(z) = i! bn+1−i, i = 1, 2, . . . , n (5.6)
dz
z=0

An advantage is that it is possible to obtain the value of the coefficients from bk from
ck by fitting p(z) to the data set (xk − a, yk ), k = 1, 2, . . . , n + 1 in MATLAB® with the
following instruction:

b = polyfit(x-a, y, lenght(xd)-1)
5.2 Differentiation of Interpolation Polynomials 105

where x is a vector of xi , i = 1, 2, . . . , n + 1 and y is a vector of yi , i = 1, 2, . . . , n + 1.

In addition, it is possible to compute all derivatives of an interpolation polynomial fitted
to a data set using the following algorithm:
Algorithm 5.1 Algorithm to obtain the derivatives of a polynomial interpolation in
MATLAB® .

function der = derivapoly(xd, yd,a)

m = length(xd)-1;
b = polyfit(xd-a, yd,m);%fit of points to the data set
c =b(m:-1:1);
fact(1)=1;
for i =2:m
fact(i) = i*fact(i-1); %factorial calculation
end
der = c.*fact; %obtaining the coefficients of the derivative
end

where xd are the abscissae of the points fitted by the interpolation, yd are the ordinates
of the points fitted by the interpolation, and a is the value at which the derivative will be
evaluated.
Example 5.1: Taking into account a given data set:

xd yd
0 0.2838
0.25 0.4583
0.50 0.5893
0.75 0.6852
1 0.7625
1.25 0.8964

Estimate the derivatives for a = 0.2.

Solution:

xd = [0 ,0.25, 0.5,0.75,1,1.25];
yd = [0.2838,0.4583,0.5893,0.6852,0.7625,0.8964];
a = 0.2;
der = derivapoly(xd, yd,a)
106 5 Numerical Differentiation Methods

Resulting in:

right =

0.6315 -0.7084 1.6968 -13.4707 51.8144

5.3 Differences Approximation

In the previous section, we worked out a simple way to evaluate all the derivatives of
an interpolation polynomial fitted to a data set. On the other hand, in practice, it is not
usually necessary to obtain all the derivatives; instead, one or two low-order derivatives
are usually needed. The method of obtaining the values for these derivatives is called
difference approximation by polynomial interpolation.
To facilitate the deduction of the difference approximation method, we consider a
function f (x), as shown in Illustration 5.1, and assume that we wish to obtain the first
derivative of f (x) at x = 0. Knowing the values of f −1 = f (x0 − h), f 0 = f (x0 ) and
f 1 = f (x0 + h), where h is an interval between two consecutive points on the axis x axis,
it is possible to approximate f 0 '(x) = f '(x0 ) using the gradient of the linear interpolation
A, B or C, illustrated in Fig. 5.1.

(a) (b) (c)

Fig. 5.1 Different difference approximations for f'(x0 ). a Forward difference approximation,
b backward difference approximation, c central difference approximation
5.4 Differentiation by Taylor Expansion Method 107

Such approximations using the gradient of lines A, B and C are called the for-
ward difference approximation, backward difference approximation and central difference
approximation, respectively, and are mathematically expressed as follows:
Forward Difference Approximation (A):
f1 − f0
f0 ' = (5.7)
h
Backward Difference Approximation (B):
f 0 − f −1
f0 ' = (5.8)
h
Central difference approximation (C):
f 1 − f −1
f0 ' = (5.9)
2h
These approximations are closely related to the interpolation polynomials. Consider
n+1 abscissa points, x1 , x2 , . . . , xn+1 and their corresponding ordinates f 1 , f 2 , . . . , f n+1 .
It is possible to express a Lagrange interpolation polynomial fitted to the set of data points.
This results in a polynomial of order n fitted to n + 1 points so that the derivatives can
be evaluated up to their n-th order.

5.4 Differentiation by Taylor Expansion Method

An alternative to obtain the difference approximation is the Taylor expansion method,

which, in addition to deriving the difference expressions, also obtains the error terms in
them.
Assuming that it is desired to obtain a derivative of order n requires at least a number
of data points to derive n + 1 data points to derive a difference approximation. That is, to
obtain a first-order derivative of a function, at least two points are needed.
Consider the deduction by approximation of differences for f i ' = (xi ) in terms of
f i = ( f (xi ) and f i+1 (xi+1) . The form of the Taylor expansion for f i+1 around xi is
given by:

h 2 '' h 3 ''' h 4 ''''

f i+1 = f i + h f i ' + f + f + f + ... (5.10)
2 i 3 i 24 i
Solving 5.9 for f i ', we obtain:

f i+1 − f i h h 2 '''
fi ' = − f i'' + f − ... (5.11)
h 2 6 i
108 5 Numerical Differentiation Methods

Ignoring all the terms on the right-hand side of the equation except for the first one, we
obtain the difference approximation, which turns out to be the same as Eq. 5.6. Truncating
the whole expression 5.10 leads to a computational error represented in its largest weight
2
by − h2 f i'' This is because the remaining terms tend to be negligible when reducing the
value of h. Thus, the forward difference equation can be expressed as follows:
f i+1 − f i
fi ' = +E (5.12)
h
Being:
h ''
E =− f (5.13)
2 i
This tells us that an approximation of the error E is proportional to the interval given
by h and the second derivative f '' .
In a similar manner, the backward difference approximation is obtained using the
Taylor expansion for f i−1 and f i , leaving the expression as:
f i − f i−1
fi ' = +E (5.14)
h
where the calculation of the error E is the same as in Eq. 5.12.
Finally, for the calculation of the central difference approximation, the following are
used f i+1 y f i−1 which by subtracting the second term from the first term and performing
its Taylor expansion and calculating its derivative we obtain:
f i+1 − f i−1 1
fi ' = − h 2 f i'' + . . . (5.15)
2h 6
Thus, leaving the deduction of the approximation including the term E as follows:
f i+1 − f i−1
fi ' = +E (5.16)
2h
To:
1
E = − h 2 f i'' (5.17)
6

5.5 Approximation of Differences in Partial Derivatives

A partial derivative of a multivariable function is nothing more than the derivative of the
function with respect to each of the variables keeping the rest as a constant in the function.
The difference approximation method for partial derivatives is in essence the same as the
methods of numerical differentiation of one-dimensional functions.
5.6 Numerical Computation of Higher-Order Derivatives 109

Consider a two-dimensional function f (x, y). The calculation of the partial derivative
with respect to x by means of the difference approximation method is deduced by taking
y as a constant value y0 and considering f (x, y0 ) as a one-dimensional function. The
approximations previously seen (A, B, C) are expressed as follows:

f (x0 + ∆x, y0 ) − f (x0 , y0 )

fx =
∆x
f (x0 , y0 ) − f (x0 , −∆x, y0 )
fx =
∆x
f (x0 + ∆x, y0 ) − f (x0 − ∆x, y0 )
fx = (5.18)
2∆x
where f x is the partial derivative of f with respect to x.

5.6 Numerical Computation of Higher-Order Derivatives

The usefulness of the methods discussed in the previous sections is greatly reduced as
the order of the derivative to be evaluated increases. This is because the accuracy of the
difference approximations is quickly lost due to the truncation and rounding errors of the
function.
An alternative to deal with higher-order derivatives is the Cauchy integral method,
since it has very good accuracy for the calculation of high- and low-order derivatives.
This integral is mathematically expressed as follows:
∫
k! f (z)
f (k) (z 0 ) = dz (5.19)
2πi C (z − z 0 )k+1

where z is a complex variable. The integral along a closed curve C in the complex length
in which f (z) is an analytic function and z 0 is contained. The integration is performed
counterclockwise. The shape of the curve C is arbitrary, making it convenient that it is a
circle.
The circle centered at z 0 with radius r in the complex plane is given by

z = r ei π θ + z 0 , 0 ≤ θ ≤ 2π (5.20)

where θ is the angle traveled. Substituting this in 5.17, we obtain:

∫ 2π ( iπθ )
(k) k! f re + z0
f (z 0 ) = dθ (5.21)
2πi 0 ei (k+1)π θ
Once the trapezoidal rule of numerical integration has been applied (which will be
explained in more detail in Chap. 6), Eq. 5.18 can be expressed as follows:
110 5 Numerical Differentiation Methods

( )
k! ∑
N
f r ei πθ + z 0
(k)
f (z 0 ) = ∆θ (5.22)
2πi ei(k+1)πθn
n=1

where:
2π
∆θ = (5.23)
N
Y N is the number of intervals in the numerical integration along the circle.
The following algorithm presents a MATLAB® implementation of Eq. 5.20:
Algorithm 5.2 Implementation of the Cauchy Integral for numerical computation of
higher order derivatives in MATLAB®

function drift_f = CauchyI(name,z0,k)

N = 2480 + k*10;%Intervals of the Cauchy Integral
r = 1; %circle radius
delta = 2 * pi/N; %delta
teta = 0:delta:2*pi-delta; %values of angle teta
z = r * exp(1i*teta) + z0; %analytic function z
kf = 1;
for m=1:k
kf = kf*m;
end
sum_f = delta * sum(feval(name,z)./exp(1i*k*teta));%summatory
of the function

derivative_f = kf/(2*pi)*sum_f/r^k;%derivative of function

end

where name is the name of the function to be evaluated enclosed in quotation marks, x
is the abscissa of the function and k is the order of the derivative.
Example 5.2: Find the first ten derivatives of the function sen(x) for x = 0 using the
Cauchy integral in MATLAB® .
Solution:

for k = 1 : 10
drift_f = CauchyI('sin',0,k);
fprintf('k=%2d, real(drift_f)=%12.5e,
imag(drift_f)=%12.5e',....
k, real(drift_f), imag(drift_f))
end
References 111

Throwing as result:

k= 1, real(drift_f)= 1.00000e+00, imag(drift_f)=-6.42862e-18

k= 2, real(drift_f)= 6.10179e-17, imag(drift_f)=-3.33625e-17
k= 3, real(drift_f)=-1.00000e+00, imag(drift_f)= 4.62093e-17
k= 4, real(drift_f)=-2.79142e-16, imag(drift_f)=-2.71302e-16
k= 5, real(drift_f)= 1.00000e+00, imag(drift_f)=-5.88215e-16
k= 6, real(drift_f)=-1.33374e-13, imag(drift_f)=-1.70385e-14
k= 7, real(drift_f)=-1.00000e+00, imag(drift_f)= 3.47972e-14
k= 8, real(drift_f)=-7.18675e-13, imag(drift_f)=-4.19637e-13
k= 9, real(drift_f)= 1.00000e+00, imag(drift_f)= 4.99286e-12
k=10, real(drift_f)=-3.86013e-10, imag(drift_f)=-2.85811e-11

Here, real(drift_f) and imag(drift_f) are the real and imaginary parts,
respectively, resulting from the Cauchy integral calculation. It is desirable that the imag-
inary part of the calculation is equal to zero since this comes from a calculation error
either of truncation or rounding. It can also be seen that as the order of the derivative
increases, this error increases, although it is still somewhat negligible.
The accuracy of this method will strongly depend on the number of points to be
evaluated and the radius chosen for the curve.
It is also necessary to consider that the result of the Cauchy integral will be invalid in
the case of any singularity inside the circle.

References

1. M. D. preprint math/0306092 and undefined 2003, “Formulae of numerical differentiation,”

arxiv.org, 2006, Accessed: May 30, 2023. Available: https://arxiv.org/abs/math/0306092.
2. B. Fornberg, “Numerical Differentiation of Analytic Functions,” ACM Transactions on Mathe-
matical Software (TOMS), vol. 7, no. 4, pp. 512–526, Dec. 1981, doi: https://doi.org/10.1145/
355972.355979.
3. J. Cheng, X. Jia, Y. W.-I. problems in S. and, and undefined 2007, “Numerical differentiation and
its applications,” Taylor & Francis, vol. 15, no. 4, pp. 339–357, Jan. 2007, doi: https://doi.org/10.
1080/17415970600839093.
4. K. Owolabi and A. Atangana, “Numerical methods for fractional differentiation,” 2019,
Accessed: May 30, 2023. [Online]. Available: https://link.springer.com/content/pdf/https://doi.
org/10.1007/978-981-15-0098-5.pdf.
5. M. A. Abdelhakem et al., “A numerical method based on Legendre differentiation matrices for
higher order ODEs,” researchgate.net, 2020, doi: https://doi.org/10.18576/isl/090303.
6. Numerical methods differentiation—Google Académico. https://scholar.google.com/scholar?
start=10&q=numerical+methods+differentiation&hl=es&as_sdt=0,5 (accessed May 30, 2023).
Numerical Integration Methods
6

Abstract

During this chapter, the topic of numerical integration methods is discussed. As the
inverse operation of a derivative, the integration methods are crucial complements
when solving engineering problems. As in the previous sections, addressing the prob-
lem from a computational point of view involves the use of numerical methods. During
the rest of the chapter, the trapezoidal, Simpson’s 1/8 and 3/8 rules, integration using
the closed Newton–Cotes formula, numerical integration in a two-dimensional domain,
and MATLAB® integration commands are discussed.

6.1 Introduction

In the previous section, we observed in detail the numerical methods to obtain an accept-
able approximation of the calculation of the derivative of a function. On the other hand,
a mathematical tool that, in addition to being the inverse operation of the derivative, is a
crucial complement when solving engineering problems is the calculation of the integral
of a function.
Analogous to the previous sections, addressing this problem from a computational
point of view involves the use of numerical integration methods. These methods allow us
to integrate functions defined analytically or given in tabular form. The solving princi-
ple of these methods consists of fitting a polynomial to data points of the function and
then integrating it. In this way, it is possible to derive different integration methods by
modifying the distribution of the abscissae of data points.

© The Author(s), under exclusive license to Springer Nature Switzerland AG 2024 113
E. Cuevas et al., Computational Methods with MATLAB® , Synthesis Lectures
on Engineering, Science, and Technology, https://doi.org/10.1007/978-3-031-40478-8_6
114 6 Numerical Integration Methods

Throughout this chapter, we will analyze the methods of the trapezoidal rule (used
to solve the Cauchy integral in the previous chapter), Simpson’s rule and the general
deduction of the Newton–Cotes formulas. Finally, we will explain how to obtain two-
dimensional integrals numerically, and we will show a MATLAB® command to calculate
them.

6.1.1 Trapezoidal Rule

One of the most widely used methods for calculating numerical approximations of definite
integrals is the trapezium rule. It is the first of the Newton–Cotes quadratures based on
the integration of the linear integration formula used when the interpolation polynomial
is a polynomial of the first degree.
The integral of a function f (x) is evaluated as:
 b
I ( fx ) = f (x)d x (6.1)
a

It is possible to approximate f (x) by linear interpolation:

b−x x −a
p(x) = f1 + f2 (6.2)
b−a b−a
For f 1 = f (a) y f 2 = f (b).
Substituting Eq. (6.1), we can express:
 b  b h
I ( fx ) = f (x)d x ≈ p(x)d x = ( f1 + f2 ) (6.3)
a a 2
For h = b − a.
Thus, Eq. (6.3) is left as follows:
 b h
I ( fx ) = f (x)d x = ( f1 + f2 ) + E (6.4)
a 2
where E is a truncation error. Figure 6.1 graphically shows the trapezium rule. where
the area under the linear interpolation curve p(x) is the integral calculated using the
trapezoidal rule, while the area under the linear interpolation curve is the exact value of
this integral. y = f (x) is the exact value of the trapezoidal rule. Consequently, the error
given in Eq. (6.1) is equal to the area between p(x) and f (x) and is approximately:

1 3
E≈ h f  (6.5)
12
6.1 Introduction 115

Fig. 6.1 Trapezoidal rule

Furthermore, Eq. (6.4) can be expanded in multiple intervals. This is done by express-
ing the function to be integrated through n +1 data points with equidistant abscissa points.
This allows Eq. (6.4) to be applied repeatedly to each interval. Leaving the extended
equation of the trapezoidal rule as:
 b h
I ( fx ) = f (x)d x = ( f 1 + 2 f 2 + · · · + 2 f n + f n+1 ) + E (6.6)
a 2

where: h = b−a
n , x i = a + (i − 1)h, and f i = f (x i ) for i = 1, 2, . . . , n + 1 (Fig. 6.2).
The error of the trapezoidal rule extension is given by:
b − a 2 
E ≈− h f (6.7)
12
Or equivalently:

(b − a)3 
E ≈− f (6.8)
12n 2

where f  is the mean of f (x) for a < x < b.

Assume that f is an array of f i points of equispaced abscissae in an interval h. It is
possible to obtain the calculation of the extended trapezoidal rule in MATLAB® using
the following algorithm.
116 6 Numerical Integration Methods

Fig. 6.2 Trapezoidal rule extension

Algorithm 6.1 Algorithm to obtain the extended trapezoidal rule integral in MATLAB®

function I = InTrap('fun_name',a,b,n)
h = (b-a) /n;
x = a+(0:n)*h;
f = feval(f_name,x);
I = h * (sum(f)-(f(1)f(length(f)))/2);
end

where fun_name is the name of the function to be integrated, a and b are the limits
of integration and n is the number of equispaced intervals used in the extension of the
reglatrapezoidal.
6.1 Introduction 117

6.1.2 Simpson’s Integration Rule

Simpson’s integration rule is a numerical method for the calculation of definite integrals
in which, unlike the trapezoidal rule, the approximation of the intervals is through second-
degree polynomial interpolation.
There are two complementary methods of Simpson’s rule: Simpson’s 1/3 and 3/8
rules.
Simpson’s Integration Rule of 1/3
From Eq. (6.1), it is possible to obtain a second-order interpolation polynomial with three
data points at x1 = a, x2 = (a + b)/2 and x3 = b. Once the points have been evaluated in
the functions with f 1 , f 2, and f 3 . Simpson’s rule of 1/3 can be deduced by substituting
the quadratic polynomial by f (x) of Eq. (6.1), and it is expressed as (Fig. 6.3):

h
I S1/3 ≈ ( f1 + 4 f2 + f3 ) + E (6.9)
3

where
b−a
h= (6.10)
2
and E is the error given by:

Fig. 6.3 Simpson’s 1/3 rule

118 6 Numerical Integration Methods

h5
E= f  (6.11)
90
The extension of Simpson’s 1/3 rule consists of repeatedly applying Eq. (6.9) over a
space divided into an even number of intervals n. Thus, the expression of the expansion
of Simpson’s rule of 1/3 is as follows:
 b h
I E S1/3 = f (x)d x ≈ ( f 1 + 4 f 2 + 2 f 3 + 4 f 4 + · · · + 2 f n−1 + 4 f n + f n+1 ) + E
a 3
(6.12)

where

f i = f (a + (i − 1)h) (6.13)

And
b−a
h= (6.14)
n
Finally, the calculation of the error is given by:

h 4 
E = −(b − a) f (6.15)
180

where f  is the mean of f  for a < x < b.

Example 6.1
Evaluate the following integral, using the extension of Simpson’s rule from 1/3 to n = 2, 4
and 8 knowing that the exact result is 4.006994.
 2√
I = 1 + ex d x
0
6.1 Introduction 119

Solution

Algorithm 6.2 Algorithm to obtain the integral by extension of Simpson’s 1/3 rule in
MATLAB®

clear
clc
Ireal = 4.006994;
a = 0;
b = 2;
fprintf('Simpson's Rule Extension of 1/3')
fprintf(' n IS13 E)
n = 1;
for k = 1:3
n = 2*n;
h = (b-a)/n;
i = 1:n+1;
x = a + (i-1)*h;
fun = sqrt(1+exp(x));
IS13 = (h/3)* (fun(1)+4*sum(fun(2:2:n))+fun(n+1));
if n > 2
IS13 = IS13 + (h/3)*2*sum(fun(3:2:n));
end
fprintf('%3.0f %10.5f %10.5f %10.5f',n,IS13,Ireal-IS13);
end

Resulting in:

%Extension of Simpson's 1/3 Rule

n IS13 E
2 4.00791 -0.00092
4 4.00705 -0.00006
8 4.00700 -0.00000

Simpson’s 3/8 Rule

One of the main problems when working with the extension of Simpson’s rule of 1/8 is
when working with equispaced intervals, specifically when working with an odd number
of intervals. In this case, the process to follow is to apply Simpson’s rule of 3/8 to either
the first three or the last three intervals and then apply the extension of Simpson’s rule of
120 6 Numerical Integration Methods

1/8 to the rest. This is possible because the order of the error of both rules is the same
and there is no loss of accuracy in the integral calculation.
Similar to Simpson’s 1/8 rule, Simpson’s 3/8 rule is based on the fit of an interpolation
polynomial only, this time, a third-degree polynomial expressed as:
3h  
I S3/8 ≈ f1 + 3 f2 + 3 f 3 + f4 + E (6.16)
8
with
b−a
h= (6.17)
3

f i = f (a(i − 1)h)

and a mistake:
3h 5
E= f  (6.18)
80
The following algorithm can be used to perform the calculation of an integral using
Simpson’s quadrature by checking whether the number of intervals is even or odd, in
which case first apply the 3/8 rule and then the 1/3 rule to the remainder, as explained
in the previous examples.

Algorithm 6.3 Algorithm to obtain the integral by Simpson’s method for different intervals
in MATLAB®
6.1 Introduction 121

function int = Simpson('funname',a,b,n)

h = (b-a)/n;
x = a +(0:n)*h;
fun = feval(funname,x);
if n ==1
fprintf("The data has only one interval");
return;
end
if n ==2
int = h/3*(fun(1) + 4*fun(2) + fun(3));
return;
end
if n == 3
int = 3/8 * h * (fun(1) + 3 * fun(2)....
+ 3*fun(3) + fun(4));
return;
end
if 2* floor(n/2)~= n
int = 3/8* h * (fun(n-2) + 3*fun(n-1)...
+ 3*fun(n)+fun(n+1));
m = n-3;
else
m = n;
end
int = int + (h/3) * (fun(1) + 4*sum(fun(2:2:m))...
+ fun(m +1));
if m > 2
int = int + (h/3)*2*sum(f(3:2:m));
end
end

Another algorithm in which the number of intervals is calculated automatically is

presented as follows.

Algorithm 6.4 Algorithm to obtain the integral by Simpson’s method for different intervals
in MATLAB®
122 6 Numerical Integration Methods

function int = simpsv2(f,h)

n = length(f)-1;
if n ==1
fprintf("The data has only one interval");
return;
end
if n ==2
int = h/3*(fun(1) + 4*fun(2) + fun(3));
return;
end
if n == 3
int = 3/8 * h * (fun(1) + 3 * fun(2)....
+ 3*fun(3) + fun(4));
return;
end
int = 0;
if 2* floor(n/2)~= n

int = 3/8* h * (fun(n-2) + 3*fun(n-1)...

+ 3*fun(n)+fun(n+1));
m = n-3;
else
m = n;
end
int = int + (h/3) * (fun(1) + 4*sum(fun(2:2:m))...
+ fun(m +1));
if m > 2
int = int + (h/3)*2*sum(f(3:2:m));
end
end

6.1.3 Integration by Closed Newton–Cotes Formula

In addition to the trapezoidal rule and Simpson’s rules, there are several ways to perform
the numerical calculation of an integral from the integration of the Lagrange polynomial.
One of them is the use of the closed Newton–Cotes formula, which uses n + 1 data points
with equispaced intervals. The abscissae of these points are given by:
b−a
xi = (i − 1) + a, i = 1, 2, . . . , n + 1 (6.19)
n
Its order is given by f i = f (xi ). In general, the closed Newton–Cotes formula is
expressed as follows:
 b 
n+1
IN C = f (x)d x ≈ h wi f i (6.20)
a i=1
6.1 Introduction 123

where n + 1 is the total number of points, in addition:

b−a
h= (6.21)
n
and wi are weighting factors given by the following Table 6.1.
When n takes values n =1, 2 and 3, the closed Newton–Cotes formula becomes the
rules: Trapezoidal, Simpson’s 1/3 and Simpson’s 3/8.
The following algorithm is used to calculate the integral of an analytical function.

Algorithm 6.5 Algorithm to obtain the integral by the Newton–Cotes closed formula
method in MATLAB®

function int = NewtonCot(funname,a,b,n)

ninter = n+1;
if ninter <0
return;
end
en =ninter:-1:1;
x = 1:ninter;
%% Factor weighting calculation for the formula
for i = 1: ninter
pot2(i) = ninter^en(i);
pot1(i) = 1^en(i);
end
for j = 1:ninter
z = zeros(1,ninter);
z(j) = 1;
a1 = polyfit(x,z,ninter-1);
w(j) = sum(a1.*(pot2-pot1)./en);
end
%% Integral Calculus
x = a :(b-a)/(ninter-1):b;
y = feval(funname,x);
h = (b-a)/(ninter-1);
int = h * sum(w.*y);
%fprintf("x, y w");
%for j =1:ninter
% fprintf("%e%e%e%e",x(j),y(j),w(j));
%end
end

where funname is the name of a function file “.m” that defines the integrand, a and b are
the limits of integration and n is the number of equispaced intervals between the limits of
integration.
 124

Table 6.1 Weighting factors for the Newton–Cotes closed Newton–Cotes formula
n i
0 1 2 3 4 5 6 7 8
1 w 0.5 0.5
2 w 0.333 1.333 0.333
3 w 0.375 1.125 1.125 0.375
4 w 0.3111 1.422 0.533 1.422 0.311
5 w 0.3298611 1.3020833 0.0868055 0.868055 1.3020833 0.3298611
6 w 0.29285714 1.54285714 0.19285714 1.94285714 0.19285714 1.54285714 0.29285714
7 w 0.30422453 1.44901620 0.53593749 1.21082175 1.21082175 0.53593749 1.44901620 0.30422453
8 w 0.27908289 1.66151675 −0.26186948 2.96183421 −1.28112874 2.96183421 -0.26186948 1.66151675 0.27908289
6 Numerical Integration Methods
6.1 Introduction 125

Example 6.2
Evaluate the following integral for n = 2, 4, 6:
 π
sen(x)d x
0

Solution

Algorithm 6.6 Example solution 6.2 in MATLAB®

n = [2,4,6];
a = 0;
b = pi;
for i = 1:length(n)
fprintf("n = %d",n(i))
int = NewtonCot('sin',a,b,n(i));
fprintf("int = %2.4f",int);
end

Resulting in:

n = 2
int = 2.0944
n = 4
int = 1.9986
n = 6
int = 2.0000

6.1.4 Numerical Integration in Two-Dimensional Domains

In the previous sections, the calculation of numerical integrals for functions of a one-
dimensional domain was performed; however, in real world, functions usually belong to
several n-dimensional domains.
Throughout this section, we will consider functions with a domain, as shown in
Fig. 6.4, where the limits of integration on the left and right are vertical lines, and the
126 6 Numerical Integration Methods

Fig. 6.4 Two-dimensional

integration domain

upper and lower limits are given by analytic functions. Therefore, the integral of this
function f (x, y) is presented as follows:

 b d(x)
Ib = f (x, y)d xd y (6.22)
a c(x)

Or
¨
Ib = f (x, y)d xd y (6.23)
A

To achieve this two-dimensional numerical integration, it is sufficient to reduce the

problem to a combination of one-dimensional problems by reducing Eq. (6.22) as follows:
Defining:
 d(x)
I1 (x) = f (x, y)dy (6.24)
c(x)

Substituting in Eq. (6.22), we obtain:

 b
Ib = I1 (x)d x (6.25)
a
6.1 Introduction 127

With this, it is possible to apply any of the numerical methods of integration previously
seen. In addition, it is possible to write a numerical approximation of Eq. (6.25) as:


n+1
Ib ≈ wi I1 (xi ) (6.26)
i=1

where xi are the points of the method used and wi are weighting factors. By doing x = xi
Eq. (6.24), the result is
 d(x)
I1 (x) = f (xi , y)dy (6.27)
c(x)

Transforming the problem into a one-dimensional problem since the integration vari-
able is y makes it possible to evaluate this equation with a numerical integration
method.
The following algorithm allows us to perform a numerical integral in a two-
dimensional domain.

Algorithm 6.7 Algorithm to obtain the numerical integral in a two-dimensional domain in

MATLAB® .

function dint = intbid(funname,cx,dx,a,b,m,n)

if m<2 | n < 2
fprintf("number of invalid intervals");
return
end
''pwodgt"qh"kpvgtxcnu
interm = m+1;
inten = n+1 ;
hx = (b-a)/m;
x = a+(0:m)*hx;
for i = 1:interm
yinf = feval(cx,x(i));
ysup = feval(dx,x(i));
hy =(ysup-yinf)/n;
y(i,:) = yinf +(0:n)*hy;
f(i,:) = feval(funname,x(i),y(i,::));
Int1(i) = simpsv2(f(i,:),hy);
end
dint = simpsv(Int1,hx);
end
128 6 Numerical Integration Methods

6.1.5 MATLAB® Commands for Numerical Integration

MATLAB® presents an alternative to the algorithms previously seen for the numerical
calculation of the integral of a function, using the reserved word quad. This function uses
Simpson’s recursive rule, and its syntax is as follows:

quad('funname',a,b)

quad('funname'a,b,tolerance)

quad('funname'a,b,tolerance,trace)

Depending on the number of input arguments to the function, it can obtain the value of
Simpson’s integral with a tolerance of 0.001, obtain the recursive quadrature calculation
until the tolerance value is satisfied and even plot a graph on the screen showing the
progress of the integration as long as the plot value is nonzero.

References

1. E. Hairer, M. Hochbruck, A. Iserles, C. L.-O. Reports, and undefined 2006, “Geometric numer-
ical integration,” ems.press, Accessed: May 30, 2023. [Online]. Available: https://ems.press/jou
rnals/owr/articles/1196
2. P. Davis and P. Rabinowitz, “Methods of numerical integration,” 2007, Accessed: May 30, 2023.
[Online]. Available: https://books.google.com/books?hl=es&lr=&id=gGCKdqka0HAC&oi=
fnd&pg=PP1&dq=numerical+methods+integration&ots=NDBzvvCOaV&sig=fNauXbgjM316
NC9zU6IPHcjv0Lw
3. P. Theocaris, N. I.-Q. of A. Mathematics, and undefined 1977, “Numerical integration meth-
ods for the solution of singular integral equations,” ams.org, vol. 4, no. 1, 1992, Accessed: May
30, 2023. [Online]. Available: https://www.ams.org/qam/1977-35-01/S0033-569X-1977-044587
3-X/
4. A. Schwartz, “Theory and implementation of numerical methods based on Runge–Kutta integra-
tion for solving optimal control problems,” 1996, Accessed: May 30, 2023. [Online]. Available:
https://search.proquest.com/openview/e34a6c8f6769f31ffb595e22ce97269f/1?pq-origsite=gsc
holar&cbl=18750&diss=y&casa_token=wNCxNK82B8AAAAAA:uVwyRZlB2h3U6L9Yiuxr8
tYfWg4CaIjssCjhuHpn1KnI7VHyNTSp0DpKAHBh8Sl4EBHjSJSgfzA
Roots of Nonlinear Equations
7

Abstract

When we encounter zeros or roots of a function f(x), it indicates that we are referring
to finding the solutions of a scalar equation f(x) = 0. In this chapter, we explore sev-
eral widely used methods in mathematics to find the real roots of nonlinear equations,
including the bisection method, the secant method, the method of successive substi-
tution, Newton’s iteration method, and the graphical method. The bisection method
involves repeatedly dividing an interval until a root is found, and the secant method
approximates the root using linear interpolation between two points. In the successive
substitution method, an initial guess is iteratively refined until the equation is satis-
fied. This method can be particularly useful when dealing with equations that do not
have a closed-form solution. Newton’s iteration method is widely used for finding the
roots of nonlinear equations. It involves iteratively improving an initial guess by using
the tangent line of the function at that point. The graphical method plots the process
and visually identifies the points where it intersects the x-axis, indicating the roots.
While MATLAB provides the roots command for finding the roots of polynomials,
the additional methods discussed in this chapter offer more flexibility and precision
when dealing with nonlinear equations. Understanding and applying these methods
to find the roots of equations can greatly enhance problem-solving capabilities in
various mathematics, engineering, and scientific research fields. Section 7.1, a short
introduction to the chapter, is given. Section 7.2 describes the essential elements of the
graphical method. Section 7.3 presents the implementation of the bisection method.
Additionally, Sect. 7.4 reviews the most important features of Newton’s iteration
method. Section 7.5 gives an overview of the secant method, and Sect. 7.6 addresses
the method of successive substitutions. Finally, in Sect. 7.7, nonlinear simultaneous
equations are discussed.

© The Author(s), under exclusive license to Springer Nature Switzerland AG 2024 129
E. Cuevas et al., Computational Methods with MATLAB® , Synthesis Lectures
on Engineering, Science, and Technology, https://doi.org/10.1007/978-3-031-40478-8_7
130 7 Roots of Nonlinear Equations

7.1 Introduction

When dealing with zeros or the roots of a function f (x) we automatically identify that
we are talking about the solution of a scalar equation f (x) = 0. In this chapter, we will
discuss some of the most implemented methods within the mathematical field to obtain
the real roots of nonlinear equations; these are the bisection method, the secant method,
the method of successive substitution, the method of Newton’s iteration or the graphical
method.
In MATLAB, the roots command can be implemented, which returns the roots of a
polynomial represented as a column vector. This vector contains the coefficients of the
polynomial; however, the methods discussed in this chapter are very useful if a higher
degree of accuracy is desired.

7.2 Graphical Method

If we wish to identify the positive root of

f (x) = 0 (7.1)

where
( )
1
f (x) = xsen 3 − 0.9exp(−x) (7.2)
x

If we questioned ourselves about the resolution of the positive root of the above equa-
tion, we would possibly start with the inspection of the equation; after a few moments, we
would
( identify
) that when x approaches the value of 0 from the positive side of x, the term
sen x13 oscillates exponentially with increasing frequency and becomes singular when
x takes
( the
) value of 0. When we have understood the above, we will begin to glimpse
xsen x13 and 0.9exp(−x). To further clarify the problem, perhaps we would proceed to
make a graph, either in some computer software or we could even do it by hand, once the
graph is done, it would be possible to identify that there is only one root and it is near
0.75.
It may be thought that identifying the most accurate solution is the most important
when implementing computational software. However, another very significant section
is the graphica4l method. Knowing the behavior of the function we are evaluating, its
approximate value or the number of roots, it would be feasible to implement them in
software to easily solve the problem. ( )
It can be displayed graphically in a simple way xsen x13 and 0.9exp(−x),, as we can
see in Fig. 7.1; however, we can also plot the combined function, as shown in Fig. 7.2.
7.2 Graphical Method 131

Fig. 7.1 Graphing two functions

Inspecting the corresponding graphs, it is possible to identify that the root is approxi-
mately 0.73. Additionally, if we implement some tool that expands a selected area, we
can observe a more accurate value. Nevertheless, enforcing any of the methods addressed
in subsequent sections would be more effective.

Exercise A In this exercise, we will apply graphical methodology to solve real problems;
therefore, we will start by using an equation that models the natural vibration frequencies
of a uniform beam clamped at one end.

cos(x)cosh(x) + 1 = 0 (7.3)

where

x = pω2 L/E I (7.4)

In Eq. (7.4), L is the value of the beam length expressed in meters, ω is the frequency
defined in cycles per second, E.I. represents the bending stiffness expressed in pascals,
and p refers to the density of the beam material expressed in kilograms per cubic meter.
132 7 Roots of Nonlinear Equations

Fig. 7.2 Graphing the combined function

Code 7.1
z?"2<203<42=
{?"equ*z+0,equj*z+-3=
rnqv*z.{.z."|gtqu*ngpivj*z++.)eqnqt).)dncem)+=
zncdgn*)z)+={ncdgn*){)+=
vgzv*32.72222222.)h*z+"?"equ*z+0,equj*z+-3)+

To obtain the values of the three positive roots using the graphical methodology, we
first take Eq. (7.3); as we do not have a wide knowledge of the function, we can first
graph it without the y-value when the values of x range from 0 to 20. The appreciation
we obtain in Fig. 7.3 indicates that one of the three roots is found approximately when
x takes the value of 17.5. However, we can also observe that other roots are possibly
located in the x interval between the values 0 and 15.
By using the Axis command, we modify the edges of the graph shown in Fig. 7.3 to
identify the other three smaller positive roots. The line of code implemented is axis ([0
20 −20 20]), which allows us to identify that the other three roots are at x values of 1.8,
4.6 and 7.8, as shown in Fig. 7.4.
7.3 Bisection Method 133

Fig. 7.3 Plot of Eq. (7.3)

In the previous problem, we can see that the graphical methodology has no difficulty
or drawback; however, it is common that the visual method yields poor graphs. This
issue commonly occurs when we approach a singularity; therefore, when implementing
the graphical methodology, we must be very careful when plotting the graph to avoid
losing important features or confusing singularities with roots. When there is a possibility
of encountering a singularity, it is necessary to observe the amplification of the graph to
ensure that it exists.

7.3 Bisection Method

When we know that a root is in an interval, it is possible to implement a very versatile,

effective, and simple method such as the bisection method. The bisection method has the
ability to be implemented in nonanalytic functions [1]. However, the only requirement is
that this methodology is used after a graphical analysis.
Let us take as a reference that there is a root of f (x) = 0 in an interval that is
defined between a value an of x and a value c of x; we can rewrite this mathematically as
a ≤ x ≥ c or denote it simply as a closed set [a, c]. The bisection methodology is based
134 7 Roots of Nonlinear Equations

Fig. 7.4 Plotting Eq. (7.3) with axis command

on the statement that we take as a reference. When we have an interval that has a root,
then the sign of f (x) is different at each of the extremes; this is graphically explicit in
Fig. 7.5.
The bisection method operates first by dividing into two halves the interval where the
root is located, followed by verifying the signs of the now three functions since before;
they were two functions implementing only the ends of the interval. Nevertheless, by
dividing the interval, we now consider half as another function. Therefore, we verify both
the sign of the left end (f(a)) and half (f(b)) and half (f(b)) and right end (f(c)). In this way,
it is possible to identify in which section the root is located. Once the place where the root
is located is identified, we perform the same bisection so that the space where the root
is situated is bounded. It is necessary to indicate that in each bisection, the midpoint is
taken, in addition, when the new interval is smaller than a tolerance, which is commonly
established in 10−6 .

Exercise B Find the root of the above equation

f (x) = (1 − x ∗ cos(x)) ∗ x, 2 < x < (7.5)

7.3 Bisection Method 135

Fig. 7.5 Implementation of the bisection methodology

To graph the solution to this problem, it is possible to implement Code 7.2 to observe
graphically that the root is found approximately when the value of x is 4.7.
136 7 Roots of Nonlinear Equations

Code 7.2
hwpevkqp"dkuaitcrj*hwpevkqp.c.e."nkoak|s."
nkoafgt."rqkpv+
enh."jqnf"qhh
engct"[ac."engct"[ae
ykfaz"?"nkoafgt/nkoak|s="fz?*nkoafgt/
nkoak|s+1rqkpv=
zr?nkoak|s<fz<nkoafgt="{r?hgxcn*hwpevkqp."zr+=
rnqv*zr."{r.)eqnqt).)dncem)+="zncdgn*)z)+="
{ncdgn*)h*z+)+=
jqnf qp
{okp?okp*{r+="{ocz?ocz*{r+="ykfa{?{ocz/{okp=
{r?20,zr="rnqv*zr."{r.)eqnqt).)dncem)+
vqngtcpeg?20222223="kvankokv?52=
kv"?"2=
{ac?hgxcn*hwpevkqp.c+="{ae?hgxcn*hwpevkqp.e+=
rnqv*]c.c_.]{ac.2_.)eqnqt).)dncem)+="vgzv*c."/
203,ykfa{.)z?c)+
rnqv*]e.e_.]{ae.2_.)eqnqt).)dncem)+="vgzv*e."/
203,ykfa{.)z?e)+
kh"*{ac,{ae"@"2+
hrtkpvh*)h*c+h*e+"@"2"^p)+=
gnug
yjkng"3
kv?kv-3=
d?*c-e+14="{ad?hgxcn*hwpevkqp."d+=
rnqv*]d.d_.]{ad.2_.)<).)eqnqt).)dncem)+="
rnqv*d.2.)q).)eqnqt).)dncem)+
kh"*cdu*e/c+>?vqngtcpeg+
dtgcm
gpf
kh"*kv@kvankokv+
dtgcm
gpf
kh"kv>6
vgzv*d."ykfa{142."]pwo4uvt*kv+_+
gpf
kh"*{ac,{ad>?2+
e?d="{ae?{ad=
gnug
c?d="{ac?{ad=
gpf
z?d=
7.4 Newton’s Iteration 137

rnqv*]z"z_.]2027,ykfa{"
204,ykfa{_.)eqnqt).)dncem)+
vgzv*z."2047,ykfa{.")Hkpcn"uqnwvkqp)+
rnqv*]z"*z/ykfaz,20226+_.]2027,ykfa{"
202;,ykfa{_.)eqnqt).)dncem)+
rnqv*]z"*z-ykfaz,20226+_.]2027,ykfa{"
202;,ykfa{_.)eqnqt).)dncem)+

To instantiate the above code, it is necessary to call the function and assign the essential
parameters, which would be stipulated as follows bis graph(‘fun des’, 2, 6, 0, 7, 100),
where the parameter fun des is a MATLAB file where Eq. (7.5) is described. For further
understanding, this small file is stipulated in Code 7.3

Code 7.3
hwpevkqp"{?hwpafgu*z+
{"?"*3/z0,equ*z+0,z+=

7.4 Newton’s Iteration

When the main objective is to identify the root of an equation that does not present
linearity, it is feasible to choose to implement Newton’s iteration methodology. This
methodology can be applied in the complex domain to identify a complex root. In addi-
tion, this methodology can also be used for nonlinear simultaneous equations, which will
be discussed in more detail in its corresponding section.
Newton’s iteration methodology is derived from the Taylor series method [2]. Let us
assume that we have a problem where we must find a root of f (x) = 0. Focusing only on
the first-order Taylor series term of f (x) close to an initial estimate, x1 can be represented
as follows:

f (x) ≈ f (x1 ) + f '(x1 )(x − x1 ) (7.6)

The above equation can be considered an approximation of f (x). Equating Eq. (7.6)
to zero gives the following approximation:
138 7 Roots of Nonlinear Equations

( )
f (x1 )
x2 = x1 − (7.7)
f '(x1 )

By repeating the same process n times, we obtain

( )
f (xn−1 )
xn = xn−1 − (7.8)
f '(xn−1 )

Figure 7.6 graphically shows how the Newton iteration methodology resulting from
Code 7.3 works. Taking as reference an initial value x0 , a line is drawn that passes tan-
gentially through the point identified with the ordinate pair (x 0 , f 0 ). x1 can be identified
as the intersection of the tangent line with the x-axis. With the x-axis, we now take the
point identified with the ordinate pair (x 1 , f 1 ), draw a tangent line to the new point, and
then iteratively implement the same process by taking the most updated point at each
iteration.

Fig. 7.6 Newton’s iteration methodology implementation

7.4 Newton’s Iteration 139

Code 7. 4
hwpevkqp"z"?"pgyaitcrj*hwpevkqp."uvoakpk."
okp.ocz."rqkpvu+'*hwpevkqp."uvoakpk."okp.ocz."
rqkpvu+
enh."jqnf"qhh
fgnaz?20223=
ykfaz"?"oczkowo"/ okpkowo="fz?*oczkowo/
okpkowo+1rqkpvu=
zr?okpkowo<fz<oczkowo="{r?hgxcn*hwpevkqp."zr+=
rnqv*zr."{r.)eqnqt).)dncem)+="zncdgn*)z)+="
{ncdgn*){)+=
jqnf"qp
{okp?okp*{r+="{ocz?ocz*{r+="ykfa{?{ocz/{okp=
{r?20,zr="rnqv*zr."{r.)eqnqt).)dncem)+
z?guvoakpk="zd?z-;;;="p?2=
yjkng"cdu*z/zd+@20222223
kh"p@522
dtgcm
gpf
{?hgxcn*hwpevkqp.z+=rnqv*]z.z_.]{.2_.)<).)eqnqt).)dncem
)+="rnqv*z.2.)q).)eqnqt).)dncem)+
kh"p>5~p@6."vgzv*z."ykfa{142.]pwo4uvt*p+_+."gpf
{aftkx?*hgxcn*hwpevkqp.z-fgnaz+/{+1fgnaz=
zd?z=
z?zd/{1{aftkx="p?p-3=

'rnqv*]zd.z_.]{.2_.)eqnqt).)dncem)+
gpf
rnqv*]z"z_."]2027,ykfa{"204,ykfa{_.)eqnqt).)dncem)+
vgzv*z."204,ykfa{.")Hkpcn"uqnwvkqp)+
rnqv*]z"*z/ykfaz,20226+_.]2027,ykfa{"
202;,ykfa{_.)eqnqt).)dncem)+
rnqv*]z"*z-ykfaz,20226+_.]2027,ykfa{"
202;,ykfa{_.)eqnqt).)dncem)+

To instantiate the above code, it is necessary to call the function and assign the
necessary parameters, which would be stipulated as follows new graph(‘fun des’, 2,
0, 5, 50), where the parameter fun des is a MATLAB file describing the equation
y = (0.01x + 1)sen(x) − (x − 0.01) − 0.0096 for a better understanding. This small
file is stipulated in Code 7.5.
140 7 Roots of Nonlinear Equations

Code 7.5
hwpevkqp"{?hwpafgu*z+
{?*2023,z-3+0,ukp*z+/*z/2023+0",*z0`4-3+0`*/3+/2022;8=

Exercise C Considering a partition wall, we do not know the temperature of its outer part,
but we know that the wall has a thickness of 0.05 m and that the temperature of its inner
part (T0 ) is 625,000. The temperature loss outside the wall due to radiation and convection
is a factor to consider. The equation representing the temperature (T1 ) is as follows:
k ( ) ( )
f (T1 ) = (T1 − T0 ) + εσ T14 − T∞
4
+ h T1 − T f = 0 (7.9)
∆x
In the above equation, ‘k’ is the thermal conductivity factor of the wall, which has a value
of 1.2 W/mK, ‘ε’ refers to the emissivity value, which has a value of 0.8, and T0 and T1
are identified as the inside and outside temperatures of the wall, respectively, which has a
temperature of 625,000 in its inner part and an unknown temperature in its outer part. In
addition, T∞ is associated with the temperature of the environment, which is implemented
with a value of 298,000; likewise, the air temperature has the same value. T f has the same
value, the factor h is the heat transfer coefficient with a value of 20 W/m2 K, the Stefan-
Boltzmann constant ‘σ ’ is defined as 5.10 × 10−8 W/m2 K and finally ∆x, is the wall
thickness with a value of 0.05 m.

We will proceed to the resolution of the problem by applying Code 7.4. Nevertheless,
before we must establish our function in a file.m to pass it as the argument to the function
new graph, this function will have to contain the explicit form each variable approached
previously and conclude with the definition of Eq. (7.9). The result of the file can be seen
in the following code:
7.5 Secant Method 141

Fun_ejer
m?304=
g?20:=
Vkph?4;:=
Vh?4;:=
j?42=
V2?847=
uki?708g/:=
vjkempguu?2027=
{?m1vjkempguu,*z/V2+-g,uki,*z0`6/Vkph`6+-j,*z/
Vh+=

Once the file where the function to be implemented is stipulated, it is necessary to

use the function new graph as we have done previously, only that the parameters will
be different given the focus of the problem. The call of the function will be as follows:
new graph (‘fun_ejer3’, 550, 400, 600, 50). These parameters, as we remember, are the
function to implement, the initial estimation, minimum coordinate of the graph, maximum
coordinate of the graph, minimum coordinate of the graph, the maximum coordinate of
the graph and the points implemented to plot the curve. The graphical result of Exercise
C can be seen in Fig. 7.7.

7.5 Secant Method

Within the methodologies discussed, some methods are considered variations of others;
in particular, the secant method is considered a variation of Newton’s iteration method-
ology. A difference approximation has been implemented to evaluate the derivative of
an established function. Nevertheless, it is necessary to establish that it is possible to
obtain the derivative of a function by implementing two consecutive previous values of
the same function. The iteratively established scheme [3], which is established from the
above concept, is as follows:
xn−1 − xn−2
xn = xn−1 − f n−1 , n = 2, 3, . . . (7.10)
f n−1 − f n−2
To start the algorithm, it is necessary to specify the value of x0 , and it is also indispens-
able to take into account the value associated with x1 , which can be set as x1 = x1 + ∆x,
where ∆x can take the value of a small number, such as 0.01, so that we can continue
with the next iteration until the established tolerance is met.
142 7 Roots of Nonlinear Equations

Fig. 7.7 Implementation of Newton’s iteration method in Exercise B

Exercise D Consider the mass ‘m’ of a bullet with a value of 0.002 kg; now imagine that
the bullet is fired vertically into the air and is reaching its terminal velocity. The equation of
the terminal velocity is defined by g M = F, where ‘M’ is the mass, ‘g’ is the acceleration
obtained due to gravity and ‘F’ is the drag force. Therefore, by encompassing the values
mentioned above, the equation can be written in this form:

(0.002)(9.81) = 1.4 ∗ 10−5 v 1.5 + 1.15 ∗ 10− v 2 (7.11)

In Eq. (7.11), the first term of the right-hand section alludes to the frictional drag force,
while the second term indicates the pressure drag force, in addition to the fact that it v
defines the value of the terminal velocity; therefore, the problem in question would be to
identify the terminal velocity employing the secant method. For this, we will implement
an initial estimation of v = 30m/s. It is important to clarify that the function to be
implemented will be the one previously defined as Eq. (7.11). By assigning to v0 a value
of 30 and v1 a value of 30.1 based on the initial estimation, we can start to iterate and
obtain the results with the equation defined previously. The terminal velocity is 377 m/s.
7.6 Method of Successive Substitutions 143

7.6 Method of Successive Substitutions

The methodology of successive substitutions refers to a wide set of methods that

lead problems to an iterative solution for nonlinear equations [4]. As discussed in the
previous sections, the secant methodology and Newton’s iteration methodology are imple-
mentations of the method of successive substitutions. The methodology of subsequent
substitutions is implemented in several numerical algorithms to find the solution to nonlin-
ear equations, including nonlinear simultaneous and differential equations. If it is possible
to consider the equation to be solved as f (x) = 0, then we can rewrite it as follows:

x = g(x) (7.12)

It is possible to be defined in an iterative method as follows:

xn = g(xn−1 ) (7.13)

Equation (7.13) identifies n as the number of iteration steps and the initial estimate we
have considered. x0 the initial estimate that we have considered thus far. The literature also
considers the successive substitution methodology a fixed-point iteration methodology.
g(x) However, it also has a notorious disadvantage since it does not converge in all cases
to an established shape. g(x) established. The following condition must be satisfied to
address the above disadvantage:
I I
I '(x) I
Ig I < 1 (7.14)

A systematic alternative to identify a form of g(x) is to establish the following

relationship:

xn = x − α f (x) (7.15)
144 7 Roots of Nonlinear Equations

Converting the above expression to an iterative environment, we have:

xn = xn−1 − α f (xn−1 ) (7.16)

where alpha is a constant, and this constant can be determined by substituting Eqs. (7.5)
into (7.3). We can see how the iteration converges when:

−1 < 1 − α f '(x) < 1 (7.17)

Or similarly

0 < α f '(x) < 2 (7.18)

Equation (7.18) makes two interesting points: the first is that it must have the same
sign as the f prime, and the second is that an optimal convergence speed is reached when
alpha is close to 1/ f '

Exercise E The determination of the critical size that a nuclear reactor must have is based on
the criticality equation. Consider the following equation as a simplification of the criticality
equation.

tan(0.1x) = 9.2e−x (7.19)

Applying the iterative methodology previously studied in Eq. (7.16) a

f(x) = tan(0.1x) = 9.2e−x (7.20)

To estimate a value of f’ considering x between 3 and 4, we have:

' f (4) − f (3)

f = = 0.40299 (7.21)
4−3
Finally, we assign the following value to alpha:
1
α= (7.22)
0.40299
With the above, we obtain the following table.
7.7 Nonlinear Simultaneous Equations 145

7.7 Nonlinear Simultaneous Equations

The need to solve nonlinear simultaneous equations is very frequent within the math-
ematical field. Therefore, this section will address two techniques to solve nonlinear
simultaneous equations. The first methodology to be addressed is the one known as iter-
ation with successive substitutions; in this method, the system of nonlinear equations is
conceived as an engineering system or natural phenomenon. It is very common that the
equations can be represented as linear if the solution obtained is very small; it should
also be clear that a nonlinear system of this type can be rewritten similarly to the linear
system except for the coefficients, which will be defined according to the solution.
The solution in its iterative form of a nonlinear system based on successive substitu-
tions can be defined as follows:

y = An−1 xn (7.23)

In Eq. (7.23), An−1 represents a matrix that has the coefficients and is defined accord-
ing to the previous solution, and it is a nonhomogeneous term that, in the same way, An−1
is defined according to the final solution and is considered the nth xn It is considered the
nth iterative solution.
As we have implemented in previous methodologies, we will apply an initial estima-
tion factor of the solution with which the coefficient matrix will be defined. Once the
coefficient matrix is defined, we can solve the equation as a linear system. Once the
solution is identified, the coefficient matrix is updated, and the equation is solved again.
If an instability occurs during the iterative process, it is allowed to apply the following
subrelaxation, where omega (ω) has a value between 0 and 1.
146 7 Roots of Nonlinear Equations

xn = ω A−1
n−1 y + (1 − ω)x n−1 (7.24)

The second methodology addressed in this section will be Newton’s iteration method.
By implementing this methodology, nonlinear equations can be represented linearly using
the Taylor series. Let us assume that a system of equations is defined as follows:

f i (x1 , x2 , . . . , xn ) = 0, i = 1, 2, . . . (7.25)

f i represents a linear function of the set of x j . By possessing an initial estimate of the

solution in question, it can be rewritten as follows:
∆

x j = x j + ∆x j (7.26)
∆

In Eq. (7.26), x j represents the initial estimate plus an unknown correction identified
with the symbol ∆x j . By expanding Eq. (7.5) to identify a bounded Taylor polynomial
∆

of the first order close to x j , the following equation can be obtained:

∑ ∂ fi ∆ ∆ ∆

∆x j = − f i (x 1 , x 2 , . . . , x n ) (7.27)
j ∂x j

The partial derivatives are identified based on the initial estimates. Equation (7.27) can
be rewritten in matrix form as follows:

J ∆x = − f (7.28)

J represents the Jacobian matrix assigned by the expression:

[ ]
∂ fi
J= (7.29)
∂x j

In addition,
⎡ ⎤ ⎡ ∆ ∆ ∆ ⎤
∆x1 f 1 (x 1 , x 2 , . . . , x n )
⎢ ⎥ ⎢ ∆ ∆ ∆
⎥
⎢ ∆x2 ⎥ ⎢ f 2 (x 1 , x 2 , . . . , x n ) ⎥
∆x = ⎢ ⎥f =⎢ ⎥ (7.30)
⎣ ... ⎦ ⎣ ... ⎦
∆ ∆ ∆

∆xn f n (x 1 , x 2 , . . . , x n )

A difference approximation can be implemented to evaluate the partial derivatives as

shown below, where σ x j has a very small value that is set arbitrarily
( ∆ ∆ ∆ ) ( ∆ ∆ )
∆

∂ fi fi x 1 , . . . , x j + σ x j , . . . , x n − fi x 1 , . . . , x j , . . . , x n
≈ . (7.31)
∂x j σxj

Exercise F The main objective of this exercise is to implement the methodology discussed
in this section to identify the solution to the following equation.
7.7 Nonlinear Simultaneous Equations 147

f 1 (x, y) = f 2 (x, y) = 0 (7.32)

The above is identified as follows, and the condition x > 0 is satisfied

( )
f 1 (x, y) = xex p(x y + 0.8) + exp y 2 − 3 (7.33)

f 2 (x, y) = x 2 − y 2 − 0.5exp(x y) (7.34)

To solve the above, we proceed to graph f 1 (x, y) = 0 and f 2 (x, y) = 0;, which can
be seen in Fig. 7.8, based on the methodology discussed in the section on contouring
two-dimensional functions in Chap. 2.

Fig. 7.8 Solution to Exercise E

148 7 Roots of Nonlinear Equations

Code 7.6
engct
enh
jqnf"qhh
z3?2<203<4=
{3?/4<203<4=
]z.{_?ogujitkf*z3.{3+=
h3?hwpaqpg*z.{+=
h4?hwpavyq*z.{+=
|?|gtqu*4.4+
eqpvqwt*z3."{3.h3.]2022."2022_.)eqnqt).)dncem)+
jqnf"qp
eqpvqwt*z3."{3.h4.]2022."2022_.)eqnqt).)dncem)+
zncdgn*)z/czku)+=
{ncdgn*){/czku)+=

fun_one.m
hwpevkqp"h"?"hwpaqpg*z.{+
h"?"z0,gzr*z0,{-20:+-gzr*{0`4+/5=
Hwpavyq0o
hwpevkqp"h"?"hwpavyq*z.{+
h"?"z0`4/{0`4/207,gzr*z0,{+=

By interpreting Fig. 7.8, it can be identified that two roots are in the positive domain
of x. The first root is located approximately when x has a value of 0.8 and y has a value
of 0.2, while the second root is found when x has a value of 1 and y has a value of −0.8.
7.7 Nonlinear Simultaneous Equations 149

Code 7.7
engct
enh
hrtkpvh*)^p)+
fz?2023=
f{?2023=
z?kprwv*)V{rg"vjg"kpkvkcn"guvkocvg"qh"z")+=
{?kprwv*)V{rg"vjg"kpkvkcn"guvkocvg"qh"{")+=
hqt"p?3<72
u?]z.{_=
zr?z-fz=
{r?{-f{=
L*3.3+?*hwpaqpg*zr.{+/hwpaqpg*z.{++1fz=
L*3.4+?*hwpaqpg*z."{r+/hwpaqpg*z.{++1f{=
L*4.3+?*hwpavyq*zr.{+/hwpavyq*z.{++1fz=
L*4.4+?*hwpavyq*z."{r+/hwpavyq*z.{++1f{=
h*3+?hwpaqpg*z.{+=
h*4+?hwpavyq*z.{+=
fu?/Lh)=
z?z-fu*3+=
{?{-fu*4+=
hrtkpvh*)p?'402h."z?'3407g."{?'3407g)."
p.z.{+
hrtkpvh*)h*3+?'3204g."h*4+?'3204g)."h*3+."
h*4+++
kh"*cdu*h*3++>302g/;"("cdu*h*4++>302;g/;+
dtgcm=
gpf
gpf

Implementing Code 7.7 twice with the initial estimate values of x being equal to 1 and
y with the value of 1 for its first implementation and x with the value of 1 and y with a
value of −1, we find that the solution is x = 07749, y = 0.1716 and x = 0.9687, y =
−0.8477
150 7 Roots of Nonlinear Equations

References

1. Eiger, A., Sikorski, K., & Stenger, F. (1984). A bisection method for systems of nonlinear equa-
tions. ACM Transactions on Mathematical Software (TOMS), 10(4), 367–377.
2. Roots of Nonlinear Equations. In: Hauser, J.R. (eds) Numerical Methods for Nonlinear Engineer-
ing Models (2009). Springer, Dordrecht. https://doi.org/10.1007/978-1-4020-9920-5_3
3. F.A. Potra, V. Ptak. Nondiscrete Induction and Iterative Processes (1984). Pitman
4. Bogoljubov, N. N., Mitropol' skij, J. A., Samojlenko, A. M., & Sneddon, I. N. (1976). Methods
of accelerated convergence in nonlinear mechanics. Delhi, India: Hindustan Publishing Corpora-
tion.
Spline Interpolation
8

Abstract

In this chapter, we will discuss the topic of spline interpolation. Spline interpolation
is a powerful technique used to approximate a smooth curve or surface that passes
through a given set of data points. This method proves particularly useful when we
have discrete data points and need to estimate the values between those points. The core
idea of spline interpolation involves constructing a spline function by piecing together
multiple polynomial segments. The goal is to ensure that the resulting curve is not
only smooth but also continuous. Typically, these polynomial segments are chosen to
be of low degree, such as cubic polynomials (polynomials of degree 3), although other
degrees can be utilized depending on the specific requirements. Within this chapter, we
will cover the most important spline interpolation methods. By studying these methods,
you will gain insight into how to select the appropriate spline and construct the desired
curve for your data. Additionally, to enhance your understanding, the chapter includes
several MATLAB programs. These programs serve as practical illustrations of how to
employ these techniques in various simple cases.

8.1 Introduction

Many observations in science and engineering are made by experiments in which physical
quantities are measured and recorded. With data, experts can use them in a variety of
ways. Often, the data are used to develop or evaluate mathematical formulas (equations)
that represent the data.

© The Author(s), under exclusive license to Springer Nature Switzerland AG 2024 151
E. Cuevas et al., Computational Methods with MATLAB® , Synthesis Lectures
on Engineering, Science, and Technology, https://doi.org/10.1007/978-3-031-40478-8_8
152 8 Spline Interpolation

This is done by curve fitting [1]. In this technique, the use of a specific model or equa-
tion, extracted by some theory, is assumed. Then, the model parameters are determined in
such a way that the model response best fits the data. Sometimes data points are used to
estimate the expected values between known points. This process is called interpolation.
Models can also be used to predict how the data may extend beyond the range in which
they were measured. This process is known as extrapolation.
A simple way to perform the adjustment is, as already seen in previous chapters of
this book, by means of polynomial adjustment. In this approach, by using a polynomial
model defined as:

f (x) = an x n + an−1 x n−1 + · · · + a1 x + a0 (8.1)

The elements an , an−1 , …, a1 and a0 define the degree or order of a polynomial, with
n being a nonnegative integer. The polynomial’s evaluation represents a curve that can be
linear or nonlinear, depending on the polynomial’s order. For instance, a linear function
is a first-order polynomial, and its graph is a straight line. Conversely, a curve represents
higher-order polynomials that are nonlinear functions. A quadratic polynomial (second
order) creates a concave-up or down curve (parabola), while a third-order polynomial
generates an inflection point, causing the curve to be concave up or down in different
regions. Generally, higher-order polynomials produce curves with more inflections.
A given set of n data points can be approximated to a polynomial curve of different
orders, up to an order of n − 1. As an example of polynomial fitting capabilities, an
attempt is made to operate on the following data: (1, 1), (2, 2), (3, 2.2), (4, 4), (5, 6), (6,
8), (7, 8.4), (8, 8.6), (9, 9) and (10, 10). The objective is to use MATLAB to obtain the
approximation that is achieved with polynomials of different orders. To do this, the data
are entered in the command line in this way:

x=1:10;

y=[1 2 2.2 4 6 8 8.4 8.6 9 10];

Then, we use the MATLAB ® function polyfit(x,y,Or), which allows us to

find the polynomial of order Or that fits the data x and y. With this function, we will
obtain the polynomials of degree 1, 2, 3, 4, 6 and 9 that approximate the data through the
following instructions:
8.1 Introduction 153

P1= polyfit(x,y,1);

P2= polyfit(x,y,2);

P3= polyfit(x,y,3);

P4= polyfit(x,y,4);

P6= polyfit(x,y,6);

P9= polyfit(x,y,9);

Once the polynomials are obtained, their level of approximation is analyzed. For this,
the response of the polynomial to a sequence of points in x from 1 to 10 in small steps of
0.1 is evaluated. To obtain the response y of each polynomial Pi to new points d, the MAT-
LAB function y = polyval(Pi,d) is used. Therefore, the answers are obtained by
the following instructions:

d=1:0.1:10;

R1=polyval(P1,d);

R2=polyval(P2,d);

R3=polyval(P3,d);;

R4=polyval(P4,d);

R6=polyval(P6,d);

R9=polyval(P9,d);

In Fig. 8.1, various polynomial orders are used for curve fitting using the same 10
data points. The graphs in the figure demonstrate that increasing the polynomial order
results in curve fitting closer to the points. It is feasible to obtain a polynomial that passes
precisely through all points, where the value of the polynomial equals that of the point
number. For n data points, the polynomial order that passes through all points is n − 1.
In the case of Fig. 8.1, it is the ninth-degree polynomial since there are 10 points.
154 8 Spline Interpolation

Fig. 8.1 The polynomial approximations obtained for grades a 1, b 2, c 3, d 4, e 6 and f 9

8.1 Introduction 155

Figure 8.1 illustrates that the same data set can be approximated using polynomials of
different orders, as shown in the (a) first, (b) second, (c) third, (d) fourth, (e) sixth, and
(f) ninth order fits. The choice of the best fitting polynomial is not straightforward and
depends on various factors, such as the data type and source, the engineering or scientific
application, and the purpose of the curve fitting. For instance, if the data contain inaccu-
racies, using a higher-order polynomial to fit the points closely may not be reasonable.
In such cases, having an idea of the curve that relates the points may suffice. Conversely,
if the data points are precise, using a higher-order polynomial for curve fitting might be
more suitable to represent the data accurately.
Although it is possible to obtain a polynomial of order n − 1 that passes through all n
data points, using higher-order polynomials for curve fitting is not always advisable. When
there are many data points, this polynomial can have a very high degree, and although it
produces exact values at all data points, it may significantly deviate between some points.
In Fig. 8.1, the ninth-order polynomial exhibits such behavior, with the curve diverging
between the first two and last two points, failing to follow the general trend of the data
points. Therefore, while a higher-order polynomial generates the exact values at all data
points, it is not reliable for interpolation or extrapolation.
When a set of n data points is provided, a single polynomial can be used for interpola-
tion between the points, producing exact values at the points and estimated (interpolated)
values between them. For small n, where the order of the polynomial is low, the interpo-
lated values are generally reasonably accurate. However, using a high-order polynomial
for interpolation with many points can lead to significant errors, as demonstrated in
Fig. 8.1, where a ninth-order polynomial is used for interpolation with ten data points.
In the figure, it is evident that the polynomial deviates significantly from the trend of the
data near the extremes, and using a ninth-order polynomial is not a reliable method for
interpolation.
To obtain a better data fit when dealing with a large number of points, it is recom-
mended to use many low-order polynomials instead of a single high-order polynomial.
Each low-order polynomial is valid only in one interval between two or more points, and
they are of the same order, but the coefficients are different in each interval. This approach
is called spline interpolation, where the points are connected by straight lines when first-
order polynomials are used and by curves when second- or third-order polynomials are
used. The points where two adjacent intervals’ polynomials meet are called nodes, and
spline interpolation gets its name from a drafter’s spline—a flexible rod used to interpo-
late over discrete points marked by pins. There are three types of spline interpolation:
linear, quadratic, and cubic.
156 8 Spline Interpolation

8.2 Linear Splines

Linear spline interpolation involves connecting data points with straight lines. The equa-
tion of the straight line between the first two points using the Lagrange format can be
seen in Fig. 8.2. Using the Lagrange format, the equation of the straight line connecting
the first two points is given by:

(x − x2 ) (x − x1 )
f 1 (x) = y1 + y2 (8.2)
(x1 − x2 ) (x2 − x1 )
For n given points, there are n − 1 intervals. Interpolation in interval i, which is between
points xi and xi+1 (xi ≤ x ≤ xi+1 ), is performed by the equation of the line joining point
(xi , yi ) with point (xi+1 , yi+1 ):

(x − xi+1 ) (x − xi )
f i (x) = yi + yi+1 (8.3)
(xi − xi+1 ) (xi+1 − xi )
In all intervals between pointsi = 1, 2, . . . , n − 1, Eq. 8.3 can be used for linear spline
interpolation. This approach provides a continuous interpolation because the two adjacent
polynomials have the same value at a common node, but there is a discontinuity in the
slope between the nodes.
To illustrate this method, the approximation to points (8, 11), (11, 9), (15, 10) and (18,
8) will be obtained by means of linear splines.
Since the fit involves four points, we will have three intervals that will be modeled
by three linear splines f 1 , f 2 , and f 1 , f 2 . These functions will be modeled according to
Eq. 8.3 by the following formulations expressed in Eqs. 8.4, 8.5, and 8.6:

Fig. 8.2 Linear splines

8.2 Linear Splines 157

(x − x2 ) (x − x1 )
f 1 (x) = y1 + y2
(x1 − x 2) (x2 − x1 )
(x − 11) (x − 8)
= 5+ 9
(8 − 11) (11 − 8)
5(x − 11) 9(x − 8)
= + (8.4)
−3 2
(x − x3 ) (x − x2 )
f 2 (x) = y2 + y3
(x2 − x3 ) (x3 − x2 )
(x − 15) (x − 11)
= 9+ 10
(11 − 15) (15 − 11)
9(x − 15) 10(x − 11)
= + (8.5)
−4 4
(x − x4 ) (x − x3 )
f 3 (x) = y3 + y4
(x3 − x4 ) (x4 − x3 )
(x − 18) (x − 15)
= 10 + 8
(15 − 18) (18 − 15)
10(x − 18) 8(x − 15)
= + (8.6)
−3 3
The models described by Eqs. 8.4, 8.5 and 8.6 characterize the interval fit. That is,
function f 1 models the data in the interval 8 ≤ x ≤ 11, function f 2 in the interval
11 ≤ x ≤ 15 and function f 3 at 15 ≤ x ≤ 18. In this way, if one wants to obtain the
estimate at point x = 12.7, one should use the function f 2 to obtain it; in this way, the
calculation is executed as follows:
9(x − 15) 10(x − 11)
f 3 (x) = +
−4 4
9(12.7 − 15) 10(12.7 − 11)
= +
−4 4
= 9.425 (8.7)

The function described in Program 8.1 generically implements the interpolation method
under the linear spline approach.
158 8 Spline Interpolation

%Program 8.1

function Yest = Spline_lineal(x,y,Xest)

% It estimates the interpolation with the linear spline method.

% Input variables:

% y Vector containing the y-coordinates

% Vector containing the x-coordinates.

% Xest point or vector in x to be estimated.

% Output value:

% Yest is the estimated value as a result of Xest.

n = length(x);

nd=length(Xest);

for j=1:nd

for i = 2:n

if Xest(j)< x(i)

break

end

end

Yest(j)=(Xest(j)-x(i))*y(i-1)/(x(i-1)-x(i))+...

...(Xest(j)-x(i-1))*y(i)/(x(i)-x(i-1));

end

To evaluate program code 8.1, we consider the points (1, 1), (2, 1.5), (3, 4.5), (4, 5)
and (5, 6) by means of the linear spline method as an example. To do this, we first enter
the data vectors in the command line:

x=1:5;

y=[1 1.5 4.5 5 6];

8.2 Linear Splines 159

Then, the vector where the interpolated values will be estimated is defined, which in
this case will be represented as a vector of data ranging from 1 to 5 in increments of 0.1.
This is achieved by defining the following vector by command line:

d=1:0.1:5;

With the data entered, the function is called, and its result is plotted by using the
following commands:

Ye = Spline_lineal(x, y, d);

plot(x,y,’o’);

hold on
plot(d,Ye);

The result is shown in Fig. 8.3.

Fig. 8.3 Interpolation of quadratic splines

160 8 Spline Interpolation

8.3 Quadratic Splines

Quadratic splines [3] use second-order polynomials for interpolation in the intervals
between points, as shown in Fig. 8.4. For n given points, there are n − 1 intervals, and the
equation of the polynomial in interval i, between points xi and xi+1 , is given by using
the standard form:

f i (x) = ai x 2 + bi x + ci (8.8)

This model is valid for the n − 1 intervals where the data curve is fitted. For curve
fitting or interpolation, in general, there are n − 1 equations, and since each equation
has three coefficients (ai , bi , ci ), a total of 3n − 3 parameters must be determined. The
coefficients are determined by applying the following conditions:

1. Each polynomial f i (x) must pass through the extremes of the interval, (xi , yi ) and
(xi+1 , yi+1 ), which means that f i (x) = yi _i and f i (xi+1 ) = yi+1 . Thus, we have:

ai xi2 + bi xi + ci = yi (8.9)

2
ai xi+1 + bi xi+1 + ci = yi+1 (8.10)

Since there are n−1 intervals, these conditions produce 2(n−1) = 2n−2 equations.
2. Within each interval of quadratic splines, the slopes of the adjacent polynomials
are equal, resulting in continuous slope for the curve passing through the interval.
Therefore, the first derivative of the ith polynomial can be expressed as:

Fig. 8.4 Quadratic splines

8.3 Quadratic Splines 161

df
f (x) = = 2ai x + bi (8.11)
dx
For n points, the first point in the interior is i = 2 and the last one is i = n − 1.
Considering the first successive derivatives at all points in the interior, it is generated
the following expression:

2ai−1 xi + bi−1 = 2ai xi + bi (8.12)

Considering that there are n − 2 points in the interior, this condition produces n − 2
equations.
Together, the two conditions described above result in 3n − 4 equations. However,
since there are n−1 polynomials, there are 3n−3 coefficients that need to be estimated.
Therefore, one more condition is required to solve for all the coefficients. A common
additional condition is to set the second derivative at either the first or last point to
zero.
3. Assuming that the second derivative at the first point (x1 , y1 ) is zero, the polynomial
in the first interval, between the first and the second point can be expressed as:

f 1 (x) = a1 x 2 + b1 x + c1 (8.13)

The second derivative of the polynomial f 1 (x) is f 1 (x) = 2a1 . When the second
derivative equals zero, it implies that a1 is zero, which means that a straight line with
constant slope connects the first two points.

The use of quadratic splines for interpolation ensures that the first derivative at the interior
points (nodes) is continuous. This requires solving a linear system of 3n − 4 equations for
the coefficients of the polynomials. In contrast, cubic splines have continuous first- and
second-order derivatives at the interior points, which enables the system to be expressed
in a way that requires only the solution of a linear system of n − 2 equations to estimate
its parameters, as will be explained in the next section.
To illustrate this method, the approximation to points (8, 5), (11, 9), (15, 10), (18, 8)
and (22, 7) will be obtained by means of linear splines. The example considers 5 points
(n = 5), so the data will be fitted through four intervals (i = 1, . . . , n − 1). Each of these
intervals is modeled by a quadratic function of the form f i (x) = ai x 2 + bi x + ci .
Thus, 4 polynomials are needed to fit all 5 points. Since each polynomial has three
coefficients, a total of 12 coefficients must be determined. The parameters to estimate
are a1 , b1 , c1 , a2 , b2 , c2 , a3 , b3 , c3 , a4 , b4 and c4 . It is important to consider that the
coefficient a_1 equals zero (as explained in condition 3). The other 11 coefficients are
determined from a linear system of 11 equations. In this way, eight equations are formu-
lated considering the condition that the polynomial in each interval passes through the
endpoints, as defined in Eqs. (8.9) and (8.10):

i =1 f 1 (x) = a1 x12 + b1 x1 + c1 = b1 8 + c1 = 5
162 8 Spline Interpolation

f 1 (x) = a1 x22 + b1 x2 + c1 = b1 11 + c1 = 9
i =2 f 2 (x) = a2 x22 + b2 x2 + c2 = a2 (11)2 + b2 11 + c2 = 9
f 2 (x) = a2 x32 + b2 x3 + c2 = a2 (15)2 + b2 15 + c2 = 10
i =3 f 3 (x) = a3 x32 + b3 x3 + c3 = a3 (15)2 + b3 15 + c3 = 10
f 3 (x) = a3 x42 + b3 x4 + c3 = a3 (18)2 + b3 18 + c3 = 8
i =3 f 4 (x) = a4 x32 + b4 x3 + c4 = a4 (18)2 + b4 18 + c4 = 8
f 4 (x) = a4 x52 + b4 x5 + c4 = a4 (22)2 + b4 22 + c4 = 7

Three additional equations are also obtained from the condition (Eq. 8.12) that the
slopes (first derivative) at the interior nodes of the polynomials in the adjacent intervals
are equal. Thus, it follows that:

i =2 2a1 x2 + b1 = 2a2 x2 + b2
b1 = 2a2 11 + b2
b1 − 2a2 11 − b2 = 0
i =3 2a2 x3 + b2 = 2a3 x3 + b3
2a2 15 + b2 = 2a3 15 + b3
2a2 15 + b2 − 2a3 15 + b3 = 0
i =4 2a3 x4 + b3 = 2a4 x4 + b4
2a3 18 + b3 = 2a4 18 + b4
2a3 18 + b3 − 2a4 18 + b4 = 0

These equations formulating all the conditions produce a system of 11 linear equations,
which in matrix form will be defined as:
⎡ ⎤⎡ ⎤ ⎡ ⎤
8 1 0 0 0 0 0 0 0 0 0 b1 5
⎢ ⎥⎢ ⎥ ⎢ ⎥
⎢ 11 1 0 0 0 0 0 0 0 0 0 ⎥⎢ c1 ⎥ ⎢ 9 ⎥
⎢ ⎥⎢ ⎥ ⎢ ⎥
⎢ 0 0 112 11 1 0 0 0 0 0 0 ⎥⎢ a2 ⎥ ⎢ 9 ⎥
⎢ ⎥⎢ ⎥ ⎢ ⎥
⎢ 0 0 152 15 1 0 0 0 0 0 0 ⎥⎢ b ⎥ ⎢ 10 ⎥
⎢ ⎥⎢ 2 ⎥ ⎢ ⎥
⎢ ⎥⎢ ⎥ ⎢ ⎥
⎢ 0 0 0 0 0 152 15 1 0 0 0 ⎥⎢ c2 ⎥ ⎢ 10 ⎥
⎢ ⎥⎢ ⎥ ⎢ ⎥
⎢ 0 0 0 0 0 182 18 1 0 0 0 ⎥⎢ a3 ⎥ = ⎢ 8 ⎥ (8.14)
⎢ ⎥⎢ ⎥ ⎢ ⎥
⎢ 2 ⎥⎢ ⎥ ⎢ ⎥
⎢ 0 0 0 0 0 0 0 0 18 18 1 ⎥⎢ b3 ⎥ ⎢ 8 ⎥
⎢ ⎥⎢ ⎥ ⎢ ⎥
⎢ 0 0 0 0 0 0 0 0 222 22 1 ⎥⎢ c3 ⎥ ⎢ 7 ⎥
⎢ ⎥⎢ ⎥ ⎢ ⎥
⎢ 1 0 −22 −1 0 0 0 0 0 0 0 ⎥⎢ a ⎥ ⎢ 0 ⎥
⎢ ⎥⎢ 4 ⎥ ⎢ ⎥
⎢ ⎥⎢ ⎥ ⎢ ⎥
⎣ 0 0 30 1 0 −30 −1 0 0 0 0 ⎦⎣ b4 ⎦ ⎣ 0 ⎦
0 0 0 0 0 36 1 0 −36 −1 0 c4 0
8.4 Cubic Spline 163

The system described by the matrix model of Eq. 8.14 can be formulated as a relation
of type Ax = b. Considering this, the system can be solved by MATLAB® using the
command.

Coefficients=A\B;

The coefficients obtained correspond to the following values:

b1 c1 a2 b2 c2 a3 b3 c3 a4 b4 c4
1.33 −5.66 −0.27 7.29 −38.4 0.05 −2.50 35.00 0.062 −2.75 37.25

With these values, the spline polynomials found for each interval will be defined as
follows:

f 1 (x) = 1.33x − 5.66 8 ≤ x ≤ 11

f 2 (x) = −0.27x 2 + 7.29x − 38.4 11 ≤ x ≤ 15
f 3 (x) = 0.05x 2 − 2.5x + 35 15 ≤ x ≤ 18
f 3 (x) = 0.062x 2 − 2.75x + 37.25 18 ≤ x ≤ 22

8.4 Cubic Spline

The cubic spline method [4] uses third-order polynomials for interpolation between data
points. For n data points, there are n − 1 intervals, and each third-order polynomial has
four coefficients. The computation of all coefficients can be time-consuming, as there are
different forms of the polynomial that can be used (standard, Lagrange, Newton). Two
different developments of cubic splines are presented in this section. The first one uses
the standard form of the polynomial, which is easier to understand and use but requires
solving a system of 4n − 4 linear equations. The second one uses a variation of the
Lagrange form, which is more sophisticated but requires the solution of a system of only
n − 2 linear equations.
164 8 Spline Interpolation

8.4.1 Cubic Splines with Standard Form Polynomials

When interpolating n given points with cubic splines, there are n − 1 intervals, and each
interval i between points xi and xi+1 is represented by a third-order polynomial in the
standard form, as shown in Fig. 8.5:

f i (x) = ai x 3 + bi x 2 + ci x + di (8.15)

In total, there are n − 1 equations, and each equation has four coefficients, so the total
number of coefficients that need to be determined is 4(n − 1) = 4n − 4. The coefficients
are determined by applying the following conditions:

1. Each polynomial f i (x) must pass through the extremes of the interval, (xi , yi ) and
(xi+1 , yi+1 ), which means that f i (x) = yi , and f i (x) = yi . Thus, it follows that:

ai xi3 + bi xi2 + ci xi + di = yi (8.16)

3
ai xi+1 + bi xi+1
2
+ ci xi+1 + di+1 = yi+1 (8.17)

These conditions will be valid for all intervals of the points 1, . . . , n − 1. Since
there are n -1 intervals, this condition produces 2(n − 1) = 2n − 2 equations.
2. The first derivative of the ith polynomial is continuous at the interior nodes, which
means that the slopes of the adjacent intervals are equal. This condition ensures that
the slope is continuous when changing from one polynomial to another. Thus, the first
derivative of the ith polynomial is:

Fig. 8.5 Cubic splines

8.4 Cubic Spline 165

d fi
f 1 (x) = = 3ai x 2 + 2bi x + ci (8.18)
dx
For n points, the first interior point is i = 2, and the last one is i = n − 1. Equating
the first derivatives at each interior point we obtain:

3ai−1 xi2 + 2bi−1 xi + ci = 3ai xi2 + 2bi xi + ci (8.19)

Since there are n − 2 interior points, this condition produces n − 2 additional

equations.
3. In the interior knots, the second derivatives of the adjacent polynomial intervals should
be equal. This implies that as the curve moves through an interior knot and transitions
from one polynomial to another, the rate of change of the slope (curvature) should be
smooth. The second derivative of the polynomial in the ith interval can be expressed
as:
d 2 fi
f i (x) = = 6ai x + 2bi (8.20)
dx2
For n points, the first interior point is i = 2, and the last is i = n − 1. Equating the
second derivatives at each interior point produces:

6ai−1 xi + 2bi−1 = 6ai xi + 2bi (8.21)

Since there are n − 2 interior points, this condition produces n − 2 additional

equations.

The combination of the three conditions generates a set of 4n − 6 equations for n − 1

polynomials, which have 4n − 4 coefficients. Therefore, two more equations or conditions
are necessary to solve for the coefficients. Typically, these additional conditions set the
second derivative to zero at the first and last points, resulting in two more equations. This
can be expressed as follows:

6a1 x1 + 2b1 = 0
(8.22)
6an−1 xn + 2bn−1 = 0

A type of cubic spline where the second derivatives at the endpoints are set to zero is
referred to as natural cubic splines. By applying all the necessary conditions, a system of
4n − 4 linear equations is formed to determine the 4n − 4 coefficients.
166 8 Spline Interpolation

8.4.2 Cubic Splines Based on Lagrange Polynomials

The cubic spline process using the Lagrangian form [5] starts with the second derivative
of the polynomial. Figure 8.6 shows the spline interpolation with cubic polynomials in
(a), the first derivatives of the polynomials and in (b) their second derivatives in (c).
The figure shows the interval i with the adjacent intervals i − 1 and i + 1. The second
derivative of a third-order polynomial is a linear function. This means that within each
spline, the second derivative is a linear function of x (see Fig. 8.6b). For interval i, this
linear function can be written in Lagrangian form as the following formulation:
 x − xi+1  x − xi 
f 1 (x) = f i (xi ) + f i (xi+1 ), (8.23)
xi − xi+1 xi+1 − xi

where the second derivative values of the third-order polynomial at the nodes of the ith
interval are f i (xi ) and f i (xi+1 ). The third-order polynomial in interval i can be deter-
mined by integrating Eq. (8.23) twice. The resulting expression contains two integration
constants. These two constants can be determined from the condition that the values of
the polynomial at their respective nodes are known: f i (xi ) = yi and f i (xi+1 ) = yi+1 .
Once the integration constants have been determined, the equation of the third-order
polynomial in the interval i is defined by:
 
f i (xi ) f i (xi+1 )
f i (x) = (xi+1 − x)3 + (x − xi )3
6(xi+1 − xi ) 6(xi+1 − xi )
 
yi f (xi )(xi+1 − xi )
+ − i (xi+1 − x)
xi+1 − xi 6
 
yi+1 f (xi+1 )(xi+1 − xi )
+ − i (x − xi ) (8.24)
xi+1 − xi 6

In each interval (xi ≤ x ≤ xi+1 ; i = 1, . . . , n − 1), Eq. (6.82) involves two unknown
elements, which are f i (xi ) and f i (xi+1 ). These are the values of the second derivative
at the endpoints of the interval. Equations can be obtained to relate the values of the
second derivatives with the n − 2 interior points by considering the continuity of the first
derivatives of polynomials between adjacent intervals at the interior points.

f i (xi+1 ) = f i+1 (xi+1 ), (8.25)

This condition is applied in Eq. 8.24 to write the expressions involving f i (x) and
f i+1 (x). To do this, we derive the expressions and substitute the derivatives in Eq. (8.25).
These operations produce (after algebraic simplification) the following equations:
8.4 Cubic Spline 167

Fig. 8.6 Third order

polynomials: a polynomial,
b first derivative and c second
derivative

  
(xi+1 − xi ) f i (xi ) + 2(xi+2 − xi ) f i (xi+1 ) + (xi+2 − xi+1 ) f i (xi+2 )
yi+2 − yi+1 yi+1 − yi (8.26)
=6 −
xi+2 − xi+1 xi+1 − xi

This is a system of n − 2 linear equations containing n unknown variables.

For n data points, there are n − 1 intervals. The cubic polynomial in each interval is
given by Eq. 8.24, which implies a total of n − 1 polynomials.
To determine the n values of the second derivative at each data point for n given data
points, the n−1 polynomials contain n coefficients, ranging from f i (x1 ) to f i (xn ). These
coefficients represent the values of the second derivatives of the polynomials at the data
points. The continuity of the second derivatives of the polynomials of adjacent intervals
168 8 Spline Interpolation

at interior nodes means that there are n values of the second derivative to be determined
for n data points.
To find the n unknown coefficients from f i (x1 ) to f i (xn ), the equations defined in
Eq. 8.26 form a system of n − 2 linear equations. To solve for the coefficients, two more
equations are needed. Typically, the second derivative at the first and last points of the
data is set to zero, which gives a total of n unknown coefficients and n − 2 + 2 = n
equations. This type of cubic spline is called a natural cubic spline.
 
f i (x1 ) = 0 and f i (xn ) (8.27)

By imposing the conditions given in Eq. 8.26, a system of linear equations with n
unknown coefficients (from f i (x1 ) to f i (xn )) is obtained. To solve this system, two
additional conditions are needed, and natural cubic splines are obtained by setting the
second derivatives at the first and last points to zero. Once these conditions are applied and
the linear system is solved, the obtained coefficients can be substituted into the equations
defined by Eq. 8.24 to obtain the cubic splines.
The form of Eqs. (8.24) and (8.26) can be simplified by defining a variable h i that
symbolizes the length of the ith interval (normally, the intervals are not of the same
length):

h i = xi+1 − xi (8.28)

and defining ai as the second derivative of the polynomial at point xi :

ai = f  (xi ) (8.29)

with these definitions, the equation of the polynomial in the ith interval is defined as:
ai ai+1
f i (x) = (xi+1 − x)3 + (x − xi )3
6h i 6h i
(8.30)
yi ai h i yi+1 ai+1 h i
+ − (xi+1 − x) + − (x − xi )
hi 6 hi 6

and the system of linear equations to be solved for the parameters of ai is defined by:

h i ai + 2(h i + h i+1 )ai+1 + h i+1 ai+2

yi+2 − yi+1 yi+1 − yi (8.31)
=6 −
h i+1 hi

To interpolate using cubic splines, first, we use the model in Eq. 8.31 to generate a
system of n − 2 equations with n − 2 unknown elements, ranging from a2 to an−1 (note
that for natural cubic splines, a1 and an are equal to zero). This system of equations
is derived from the equations in (8.31) and results in a tridiagonal system of equations
that can be efficiently solved using an appropriate numerical method. After solving the
8.4 Cubic Spline 169

system of equations, we can write the cubic polynomials for each interval using Eq. 8.30.
By using these steps, we can solve the cubic spline interpolation problem efficiently. The
following examples demonstrate the process of solving the fitting problem using cubic
splines.
To illustrate this method, the fit to points (8, 5), (11, 9), (15, 10), (18, 8) and (22, 7)
will be obtained by means of linear splines. The example considers 5 points (n = 5), so
the data will be fitted through four intervals (i = 1, . . . , n − 1). Each of these intervals is
modeled by a function of the form:

ai ai yi ai h i
f i (x) = (xi+1 − x)3 + (x − xi )3 + − (xi+1 − x)
6h i 6h i hi 6
yi+1 ai h i
+ − (x − xi ) (8.32)
hi 6

where h i = xi+1 − xi . There are four functions that involve five unknown coefficients: a1 ,
a2 , a3 , a4 and a5 . However, for natural cubic splines, a1 and a5 are zero, and the remain-
ing three coefficients are determined from a linear system of three equations defined by
Eq. 8.31. To apply this method, we need to calculate the values of h i for the given points
in the problem. They have the following values:

h 1 = x2 − x1 = 11 − 8 = 3
h 2 = x3 − x2 = 15 − 11 = 4
(8.33)
h 3 = x4 − x3 = 18 − 15 = 3
h 4 = x5 − x4 = 22 − 18 = 4

Thus, the linear system of three equations is defined as:

y3 − y2 y2 − y1
i =1 a1 h 1 + 2(h 1 + h 2 )a2 + a3 h 3 = 6 −
h2 h1
10 − 9 9 − 5
3 · 0 + 2(3 + 4)a2 + 4a3 = 6 −
4 3
= 14a2 + 4a3 = −6.5
y4 − y3 y3 − y2
i =2 a2 h 2 + 2(h 2 + h 3 )a3 + a4 h 3 = 6 −
h3 h2
8 − 10 10 − 9
4a2 + 2(4 + 3)a3 + 3a4 = 6 −
3 4
= 4a2 + 14a3 + 3a4 = −5.5
y5 − y4 y4 − y3
i =3 a3 h 3 + 2(h 3 + h 4 )a4 + a5 h 4 = 6 −
h4 h3
3a3 + 2(3 + 4)a4 + 4 · 0
= 3a3 + 14a4 = 2.5
170 8 Spline Interpolation

From this development, we have that the complete system of linear equations is defined
by the following matrix relation:
⎡ ⎤⎡ ⎤ ⎡ ⎤
14 4 0 a2 −6.5
⎢ ⎥⎢ ⎥ ⎢ ⎥
⎣ 4 14 3 ⎦⎣ a3 ⎦ = ⎣ −5.5 ⎦ (8.34)
0 3 14 a4 2.5

This system Ax = b can be solved in MATLAB through the following command

sequence:

A=[14 4 0; 4 14 3 0; 0 3 14];

b=[6.5; -5.5; 2.5];

Coefficients=A\B;

The coefficients obtained correspond to the following values:

a2 a3 a4
−0.3665 −0.3421 0.2519

Once the coefficients are determined, the spline interpolation polynomials are defined
as follows:
a1 a2 y1 a1 h 1
i =1 f 1 (x) = (x2 − x)3 + (x − x1 )3 + − (x2 − x)
6h 1 6h 1 h1 6
y2 a2 h 1
+ − (x − x1 )
h1 6
0 −0.3665 5 0·3
= (11 − x)3 + (x − 8)3 + − (11 − x)
6·3 6·3 3 6
9 −0.3665 · 3
+ − (x − 8)
3 6
= (−0.02036)(x − 8)3 + 1.667(11 − x) + 3.183(x − 8)
a2 a3 y2 a2 h 2
i =2 f 2 (x) = (x3 − x)3 + (x − x2 )3 + − (x3 − x)
6h 2 6h 2 h2 6
y3 a3 h 2 y3 a3 h 2
+ − (x − x2 ) + − (x − x2 )
h1 6 h1 6
−0.3665 −0.3421 9 −0.366
= (15 − x)3 + (x − 11)3 + − (15 − x)
6·4 6·4 4 6
8.5 Matlab Functions 171

10 −0.3421 · 4
+ − (x − 11)
4 6
= (−0.01527)(15 − x)3 + (−0.01427)(x − 11)3
+ 2.494(15 − x) + 2.728(x − 11)
a3 a4 y3 a3 h 3
i =3 f 3 (x) = (x4 − x)3 + (x − x3 )3 + − (x4 − x)
6h 3 6h 3 h3 6
y4 a4 h 3
+ − (x − x3 )
h3 6
−0.3421 0.2519 10 −0.342
= (18 − x)3 + (x − 15)3 + − (18 − x)
6·3 6·3 3 6
8 0.2519 · 3
+ − (x − 15)
3 6
= (−0.019)(18 − x)3 + (0.014)(x − 15)3 + 3.504(18 − x) + 2.504(x − 15)
a4 a5 y4 a4 h 4
i =4 f 4 (x) = (x5 − x)3 + (x − x4 )3 + − (x5 − x)
6h 4 6h 4 h4 6
y5 a5 h 4
+ − (x − x4 )
h4 6
0.2519 0 8 0.25 · 4
= (22 − x)3 + (x − 18)3 + − (22 − x)
6·4 6·4 4 6
7 0·4
+ − (x − 18)
4 6
= (0.0105)(22 − x)3 + (1.832)(22 − x)3 + 1.75(x − 18)

where the functions are valid in the following intervals:

f 1 (x) 8 ≤ x ≤ 11
f 2 (x) 11 ≤ x ≤ 15
f 3 (x) 15 ≤ x ≤ 18
f 4 (x) 18 ≤ x ≤ 22

8.5 Matlab Functions

The objective of this section is to illustrate how interpolation by means of the splines
method can be carried out using functions built in MATLAB. The idea is to first analyze
the function in charge of performing this type of interpolation and then to use it together
with other auxiliary functions in the construction of illustrative examples that allow us to
see the power of these approaches.
The main function that performs spline interpolation in MATLAB is spline. This
function has the following format:
172 8 Spline Interpolation

yy = spline(x,y,xx)

This command considers cubic spline interpolation to estimate the yy values corre-
sponding to the vector xx values of the underlying function connecting the y and x
points. If y is a matrix, then the data are taken as vector values, and the interpolation is
performed for each row of y. For this case, length (x) must be equal to size (y, 2), and
yy is of size (y, 1) and length (xx). This structure allows us to obtain the interpolation
yy of the points xx.
Another way to use the spline function is through the following formulation:

pp = spline(x,y)

Using this command, the interpolation data are stored in the structure defined by pp.
The interpolation data refer to the number of intervals (defined by the points defined by
the x and y points) and the values of the coefficients of the cubic polynomials defined
for each polynomial of the interval. Once the structure is stored, it can be used to make
interpolations by using the MATLAB function ppval defined by:

v = ppval(pp,xq)

As mentioned above, this function allows the interpolation of the points defined by xq,
considering the information of the polynomials defined in the intervals associated with the
pp structure.
Another important function for defining interpolation structures is mkpp. This function
allows the definition of different polynomials, each polynomial being only valid in a
defined interval. Its structure is defined by the following format:

pp = mkpp(breaks,coefs)

This instruction constructs a piecewise polynomial function in the PP structure from

the intervals defined in breaks and with coefficients specified in coefs. The vector
breaks are vectors of length L + 1 with strictly increasing elements representing the
beginning and end of each of the L intervals. coefs is an L × k matrix with each row
8.5 Matlab Functions 173

coefs(i,:) containing the coefficients of the terms, from highest to lowest exponent,
of the polynomial of order k in the interval [breaks(i), breaks(i + 1))].
To illustrate the use of the mkpp function, a simple example will be made. Consider
that it is desired to create a structure that allows interpolation by parts of a space from 0 to
15 divided into 3 intervals: from 0 to 4, from 4 to 10 and from 10 to 15. The polynomials
defined in each interval are expressed as follows:

0≤x ≤4 y1 = x 3 − x 2 + x + 1
4 ≤ x ≤ 10 y2 = x 2 − 2x + 53
10 ≤ x ≤ 15 y3 = −x 4 + 6x 3 + x 2 + 4x + 77

Since the highest order polynomial is 4, the entire structure will be defined as fourth
order. With this in mind, the MATLAB code that would generate the structure would be
defined by the following segment:

breaks = [0 4 10 15];

coefs = [0 1 -1 1 1; 0 0 1 -2 53; -1 6 1 4 77];

pp = mkpp(breaks,coefs)

The structure generated in pp is defined as follows:

pp = struct with fields:

form: 'pp'

breaks: [0 4 10 15]

coefs: [3x5 double]

pieces: 3

order: 5

dim: 1
174 8 Spline Interpolation

The structure defines the intervals, coefficients and general characteristics of the
interpolator by section. Once the structure containing the number of polynomials, their
coefficients and the interval in which they are valid has been defined, this structure can
be used to interpolate information.
As an example of its use, the piecewise structure will be evaluated at many points in
the interval [0,15], and the results will be plotted. At the same time, some vertical dashed
lines will be drawn at the break points (the boundaries of the intervals of validity) where
the polynomials meet. This is accomplished by executing the following code segment in
MATLAB:

xq = 0:0.01:15;

plot(xq,ppval(pp,xq))

line([4 4],ylim,'LineStyle','--','Color','k')

line([10 10],ylim,'LineStyle','--','Color','k')

After executing this code, the graph shown in Fig. 8.7 is generated.

Fig. 8.7 Use of the mkpp function to create interpolation structures

8.5 Matlab Functions 175

Once the operation and use of the ppval and mkpp functions have been explained,
some examples of the function for data interpolation and curve fitting will be made. The
information contained in this part is similar to that contained in the information provided
by Mathworks. However, this information is commented on and explained from a different
perspective.
The first example is intended to fit a sinusoidal curve. For this purpose, a set of initial
points are defined that define the interpolation intervals. For this example, we have 8
original points defined by the following code segment:

x = [0 1 2.5 3.6 5 7 8.1 10];

y = sin(x);

Considering these elements as the original ones to construct the spline, we try to deter-
mine the values of the function for a distribution of points from 0 to 10 in an equispacing
of 0.25. This operation can be performed by executing the following code segment:

xx = 0:.25:10;

yy = spline(x,y,xx);

plot(x,y,'o',xx,yy)

Once this code is executed, the graph shown in Fig. 8.8 is obtained.
As a second example, curve fitting is performed when the slopes of the end points
are known. To do this, the vector of function values can be specified with two additional
elements, one at the beginning and one at the end, to define the slopes of the start and
end points. As an example, a sequence of points (−4, 0), (−3, 0.15), (−2, 1.12), (−1,
2.36), (0, 2.36), (1, 1.46), (2, 0.49), (3, 0.06) and (4, 0) is considered. This sequence of
points can be defined by the following code segment:

x = -4:4;

y = [0 .15 1.12 2.36 2.36 1.46 .49 .06 0];

176 8 Spline Interpolation

Fig. 8.8 Use of the spline function for interpolation

Once the points are defined, the data are interpolated using the spline function,
and the results are plotted. For this, two additional values [0 and 0] must be specified as
additional information to indicate that the slopes of the initial and final points are zero. To
perform the exercise, the ppval function is used to estimate the spline fit over 101 points
within the interpolation interval. This operation is performed by executing the following
code segment:

cs = spline(x,[0 y 0]);

xx = linspace(-4,4,101);

plot(x,y,'o',xx,ppval(cs,xx),'-');

After executing this code, the graph shown in Fig. 8.9 is produced.
References 177

Fig. 8.9 Use of the spline function for interpolation with slope conditions

References

1. Ferguson, James C, Multivariable curve interpolation, J. ACM, vol. 11, no. 2, pp. 221–228, Apr.
1964
2. Bartels, Beatty, and Barsky, An Introduction to Splines for Use in Computer Graphics and Geo-
metric Modeling, 1987.
3. Davis, B-splines and Geometric design, SIAM News, vol. 29, no. 5, 1997.
4. Epperson, History of Splines, NA Digest, vol. 98, no. 26, 1998.
5. Elaine Cohen, Richard F Riesenfeld, Gershon Elber, Geometric Modeling with Splines: An
Introduction, CRC Press, 2019.
Differential Equations
9

Abstract

In this chapter, the focus is on differential equations, which are equations involving
unknown functions and their derivatives. These equations allow us to understand the
relationship between these unknown functions and how they change over time or in
response to different variables. Differential equations have widespread applications in
fields such as physics, engineering, economics, and biology, as they provide a math-
ematical framework for modeling various phenomena. The chapter explores the main
methods for solving differential equations using numerical techniques. While some
differential equations can be solved analytically with explicit formulas, many real-
world problems require numerical methods for obtaining solutions. Numerical methods
involve approximating the solution by dividing the problem into smaller steps and com-
puting the values of the unknown function(s) at each step. This approach allows us to
tackle complex problems that may not have a straightforward analytical solution. To
illustrate these numerical techniques, the chapter includes several MATLAB codes.
The codes demonstrate how to apply these techniques to solve differential equations
in different scenarios. By working through the examples, readers can gain a practical
understanding of how to implement and use numerical methods to solve differential
equations using MATLAB.

9.1 Introduction

Differential equations are mathematical models that describe the way in which variables
and their relationships of change are associated. In their association, two types of variables
are differentiated: dependent and independent. In their descriptions, differential equa-
tions show how the dependent variable or variables change according to variations of

© The Author(s), under exclusive license to Springer Nature Switzerland AG 2024 179
E. Cuevas et al., Computational Methods with MATLAB® , Synthesis Lectures
on Engineering, Science, and Technology, https://doi.org/10.1007/978-3-031-40478-8_9
180 9 Differential Equations

the independent variable. Many phenomena in the natural sciences and engineering can
be described through differential equations.
Differential equations can be classified according to three different criteria: type, order
and linearity.
There are two different kinds of differential equations: ordinary and partial differential
equations. In an ordinary differential equation, all ratios of change are given as a function
of a single variable. Examples of this type of equation include the following:
du
− 2y = x (9.1)
dx
du dv
+ −x =0
dx dx

d2 y dy
+2 + 4y = sinx
dx2 dx
On the other hand, in partial differential equations, the ratios of change are expressed
as a function of different variables; for this reason, the formulation ∂ x/∂ y, which is
understood as the variation of x in terms of y, is frequent in this type of equation. Some
examples of differential equations include models such as:
∂u ∂u
+ =0 (9.2)
∂x ∂y
∂u ∂u
y +x = 4u
∂x ∂y

∂ 2u ∂u ∂ 2u
+ 5 = 2
∂t 2 ∂t ∂x2
The differential equations according to their order are grouped considering the order
of the maximum derivative found in the expression. In this way, the differential equation
is defined as
( 2 )6
d4 y d2 y d y
+3 2 + = sinx (9.3)
dx4 dx dx2

It is of fourth order. The last way to classify differential equations is according to their
linearity. Under these conditions, there are two types of differential equations: linear and
nonlinear. A differential equation is said to be linear if it can be described generically as
follows:
dn y d n−1 y
bn (x) + b n−1 (x) + · · · + b0 (x)y = f (x) (9.4)
dxn d x n−1
9.2 Symbolic Solution of Equations 181

This formulation has two important properties: (A) The dependent variable y and all its
derivatives are raised to the power 1. Furthermore, all expressions relating to the depen-
dent variable represent a linear operation (not exponential or trigonometric relations).
(B) Each coefficient depends only on the variable x. If these properties are not satisfied,
the differential equation is said to be nonlinear. Some examples of nonlinear differential
equations are presented below:
dy
− x y = y2 (9.5)
dx
( )5
d2 y dy
− = sinx
dx2 dx
( 3 )( )
d y dy
3
+ 2y = 0
dx dx

In this chapter, we will analyze several methods to solve numerically ordinary dif-
ferential equations where all the dependent variables are related to a single independent
variable. In the first part of this chapter, we will study problems with initial conditions.
For their solution, the Runge–Kutta and predictor-correction methods will be explained.
The last part addresses the shooting method and the finite difference method for solv-
ing differential equations with boundary conditions. A differential equation with initial
conditions refers to the case in which its solution depends on its starting point. On the
other hand, the differential equation is a boundary condition differential equation if to
determine its solution, it is necessary to define its starting and end points.

9.2 Symbolic Solution of Equations

Although the objective of this chapter is to cover the most important methods for solv-
ing differential equations numerically, MATLAB® has support for solving differential
equations in symbolic or closed form. In this section, we will use this functionality. It
is important to consider that these methods are useful for simple differential equations
because if they are complex, they cannot be solved using these approaches.
Suppose we want to solve the differential equation defined as:
dy
+ 4y = e−t (9.6)
dx
The steps to follow are as follows. First, the dependent variable is defined through the
following code:
182 9 Differential Equations

syms y(t)

This instruction says that the symbolic math toolbox will be used and defines the
dependent variable. Then, the differential equation is defined in code. For this, a variable
is declared to represent it. Then, the differential equation is described. In its description,
it is important to consider that the equal = symbol of the differential equation will be a
double equal ==. This is because there will already exist an = symbol that is used to
associate the differential equation with the variable that represents it. Another aspect to
consider is that the differential operator dy/dt is defined in code as diff(y). Thus, the
definition of the differential equation would be:

ode = diff(y)+4*y == exp(-t);

In this definition, the variable ode will represent the differential equation in subse-
quent operations. To solve the differential equation, it is necessary to define the initial
condition(s). This is done with the following line of code:

cond = y(0) == 1;

Once the initial condition is established, the differential equation is solved with the
following instruction:

ySol(t) = dsolve(ode,cond)

With this, the solution y found and assigned to ySol(t) would be:

1 −t 2 −4t
y(t) = e + e (9.7)
3 3
9.3 Euler’s Method 183

9.3 Euler’s Method

The initial numerical method for the solution of ordinary differential equations is Euler’s
method, which is simple to understand and easy to program. Despite its low accuracy, this
method is still popular in the solution of differential equations. Because of these charac-
teristics, its analysis and understanding provide ideal conditions for the further treatment
of more elaborate methods. Consider the following first-order differential equation.

y'(t) + ay(t) = r (9.8)

Assuming the initial condition y(0) = y0 , the analytical solution to 9.1 is:
( r ) −at r
y(t) = y0 − e + (9.9)
a a
The solution described in 9.9 can be obtained by using conventional analytical methods
such as the Laplace transform. Finding the solution to differential equations by analyt-
ical methods is only possible in simple cases, since in most cases (at least in those of
theoretical or practical importance), there is no deterministic solution despite the use of
powerful symbolic computation packages. This is the reason why it is necessary to solve
it numerically through numerical methods.
To transform a differential equation into a formulation more suitable for its numerical
solution, the derivative operator y'(t) = dy/dt is substituted by an approximation. This
substitution is performed using the model shown in 9.10.

y(t + h) − y(t)
y'(t) = (9.10)
h
In this approximation, the variable is called the step size or increment that determines
the precision of the solutions to the differential equation. According to this approximation,
the differential equation shown in 9.8 is transformed as shown in 9.11.

y(t + h) − y(t)
+ ay(t) = r
h
y(t + h) = (1 − ah)y(t) + hr (9.11)

Thus, the solution to the differential equation is found by iterating the equation in
difference 9.11 step by step, increasing t in terms of h, starting at h = 0, where its initial
condition y(0) = y0 is involved.
184 9 Differential Equations

y(h) = (1 − ah)y(0) + hr = (1 − ah)y0 + hr (9.12)

y(2h) = (1 − ah)y(h) + hr = (1 − ah)2 y0 + (1 − ah)hr + hr

∑
2
y(3h) = (1 − ah)y(2h) + hr = (1 − ah)3 y0 + (1 − ah)k hr
k=0

These iterations form a numerical sequence y(kh), which represents the numerical
solution to differential Eq. 9.8. To exemplify the solution to a specific case, the parameters
of the differential Eq. 9.8 are considered to be determined by the values a = 1, r = 1,
and y0 = 0. Under such conditions, the analytical solution of the differential equation is
given by the model shown in 9.13.

y(t) = 1 − e−t (9.13)

Program 9.1 shows the MATLAB® coding to numerically solve the differential Eq. 9.1
with a h = 0.25.

%Program 9.1

clear all

a=1;

r=1;

y0=0;

h=0.25;

tf=2;

seq=0.0:0.25:2;

y(1)=y0;

for k=2:length(seq)

y(k)=(1-a*h)*y(k-1)+h*r;

end

plot(seq,y)
9.3 Euler’s Method 185

Fig. 9.1 Solution obtained by Euler’s method with h = 1, 0.5, 0.25 and analytical solution

The sequence found by Program 9.1 is shown in Fig. 9.1. The graph also shows the
analytical solution and the solutions for h = 1 and 0.5 for comparative purposes.
The accuracy of the solution is apparently related to the value of h. The smaller it
is, the closer it is to the real solution. The smaller the value is, the closer it is to the
real solution. However, this is not entirely true since in many cases, as will be seen later,
its small value will only make the method take longer to find the numerical solution to
the differential equation, since the rounding errors will result in a considerable difference
between the numerical Euler solution and the real solution.
When trying to solve differential equations numerically, it is often useful to manipulate
them algebraically to transform them in such a way that they acquire the following form:

y'(t) = f(t, y) (9.14)

Under this formulation, the ratios of change of the dependent variable are a function
of the independent variable and ratios of the dependent variable. Considering that the
differential equation to be solved numerically is in the format defined by 9.14, its solution
under Euler’s method is calculated by the following update equation:

yk+1 = yk + f(t, yk )withy(t0 ) = y0 (9.15)

186 9 Differential Equations

9.4 Trapezoidal Method

Another interesting method to numerically solve first-degree differential equations is the

trapezoidal approach. This approach starts from the integration of the differential equation
such that:

y'(t) = f(t, y) (9.16)

∫ tk+1
t
y(t)|tk+1
k
= y(tk+1 ) − y(tk ) = f(t, y)dt
tk
∫ tk+1
y(tk+1 ) = y(tk ) + f(t, y)dt
tk

Assuming that the value of the variables in f(t, y) is constant in the interval [tk+1 , tk ],
the integral of the final expression in 9.16 can be solved by means of the trapezoidal rule.
Applying this rule, one has the following as a solution:
h
y(tk+1 ) = y(tk ) + {f(tk , yk ) + f(tk+1 , yk+1 )} (9.17)
2
However, in the evaluation of the function f(tk+1 , yk+1 ), it is necessary to know the
value of the variable yk+1 , which at instant k is not known. To solve this problem, the
value of the variable yk+1 is replaced by the numerical approximation obtained by the
Euler method (discussed in the previous section).

yk+1 ∼
= yk + f(t, yk ) (9.18)

Considering Eq. 9.18 as an approximation, the final solution determined by the

trapezoidal method is defined as follows:
h
y(tk+1 ) = y(tk ) + {f(tk , yk ) + f(tk+1 , yk + f(tk , yk ))} (9.19)
2
This numerical solution is also known as the Heun method. This algorithm can be
considered a predict-and-correct approach. This is because the method first predicts the
value of yk+1 through Euler’s method so that yk+1 ∼ = yk + f(t, yk ) at the k-step. Then,
this value is corrected. Program 9.2 shows the MATLAB coding to numerically solve the
differential Eq. 9.1 by the trapezoidal method with a h = 0.25.
9.5 Runge–Kutta Method 187

%Program 9.1

clear all

a=1;

r=1;

y0=0;

h=0.25;

tf=2;

seq=0.0:0.25:2;

y(1)=y0;

for k=2:length(seq)

y(k)=(1-a*h)*y(k-1)+h*r;

end

plot(seq,y)

The sequence found by Program 9.2 is shown in Fig. 9.2. The graph also shows the
analytical solution and the solutions for h = 0.25 for comparative purposes.
From Fig. 9.2, it can be clearly seen how the trapezoidal method is more accurate than
the Euler method.

9.5 Runge–Kutta Method

Although the trapezoidal method has a better accuracy than the Euler method, its accuracy
is not sufficient for most real engineering problems. One of the most widely used methods
for solving differential equations numerically is the Runge–Kutta fourth-degree algorithm.
This method, although more complex, maintains excellent accuracy with a truncation error
of less than h 4 .
188 9 Differential Equations

Fig. 9.2 Solution obtained by the trapezoidal method with h = 0.25 and the analytical solution (real
solution)

Similar to the trapezoidal method, the Runge–Kutta method is intended to solve the
equation numerically:
∫ tk+1
y(tk+1 ) = y(tk ) + f(t, y)dt (9.20)
tk

In the trapezoidal method, the trapezoidal rule is used to solve the integral
∫ tk+1
tk f(t, y)dt. However, the Runge–Kutta method uses Simpson’s integration rule to solve
it. Under Simpson’s rule, an integral is calculated under the following formulation:
∫ tk+1 ( )
∼ ∆h ( )
f (x) = f k + 4 f k+1 + f k+1 (9.21)
tk 3 2

h
∆h =
2
9.5 Runge–Kutta Method 189

In Simpson’s rule 9.14, the value of f ( k+1 ) is replaced by the average value between
2
two successive values of the function ( f k2 + f k3 )/2. Under these conditions, the numeri-
cal solution of a differential equation is calculated by means of the Runge–Kutta algorithm
through the following procedure:
h
yk+1 = yk + (fk1 + 2fk2 + 2f k3 + fk4 ) (9.22)
6
where the values of fk1 , fk2 , fk3 y fk4 are calculated as follows:

fk1 = f(tk , yk ) (9.23)

h h
fk2 = f(tk + , yk + fk1 )
2 2
h h
fk3 = f(tk + , yk + fk2 )
2 2

fk3 = f(tk + h, yk + hfk3 )

To exemplify the use of the Runge–Kutta method, we will consider the differential
equation that models the phenomenon of series charging of a capacitor, as described in
Fig. 9.3. It is assumed that initially, the capacitor is unloaded (q = 0). By closing the
switch at t = 0, current I begins to flow through the circuit. In this process, the capacitor
C accumulates charge q. Once the capacitor acquires its maximum charge, the current
disappears from the circuit (I = 0) as if the switch were reopened. This phenomenon is
modeled by the following differential equation:
dq q
R = VE − (9.24)
dt C
Assuming that the values of the parameters of the differential equation described in
9.17 are R = 2, C = 0.8 and VE = 9, simulate by the Runge–Kutta method the loading
phenomenon for tf = 5. Program 9.3 shows the MATLAB coding to numerically solve the
differential Eq. 9.19 by the trapezoidal method with a h = 0.25. The numerical solution
found by Program 9.3 is shown in Fig. 9.4.
190 9 Differential Equations

Fig. 9.3 Electric circuit exhibiting the phenomenon of charging a capacitor

Fig. 9.4 Solution obtained by the Runge–Kutta method with h = 0.25

9.6 Predictor–Corrector Method 191

%Program 9.3

clear all

h=0.25;

VE=10;

R=2;

C=0.8;

t0=0.0;

tf=10;

y0=0.0;

t=t0:h:tf;

y(1)=y0;

for i=2:length(t)

k1=h*((VE/R)-(y(i-1)/(R*C)));

k2=h*(((VE/R)-(y(i-1)/(R*C)))+(k1/2));

k3=h*(((VE/R)-(y(i-1)/(R*C)))+(k2/2));

k4=h*(((VE/R)-(y(i-1)/(R*C)))+k3);

y(i)=y(i-1)+(k1+2*k2+2*k3+k4)/6;

end

plot(t,y)

9.6 Predictor–Corrector Method

The most popular method based on the prediction and correction approach is the Adams–
Bashforth method. The method involves three processes: prediction, correction and
modification. The first process aims to approximate the function f(t, y) by a third-order
Lagrange polynomial, such that the relations are four prior points.

(tk−3 , fk−3 ), (tk−2 , fk−2 ), (tk−1 , fk−1 ), (tk , fk ) (9.25)

192 9 Differential Equations

Once the polynomial approximation is obtained, it is substituted into the integral of

the formulation shown in
∫ tk+1
y(tk+1 ) = y(tk ) + f(t, y)dt (9.26)
tk

Thus, the prediction p(tk+1 ) of y(tk+1 ) is determined by the following expression:

∫ h
p(tk+1 ) = y(tk ) + f(t, y)dt (9.27)
0
∫ h h
f(t, y)dt = (−9fk−3 + 37fk−2 − 59fk−1 + 55fk )
0 24
In the correction process, the obtained prediction p(tk+1 ) is added as one of the new
points to obtain a new polynomial approximation of four points. Thus, the new set of
points will be composed of:

(tk−2 , fk−2 ), (tk−1 , fk−1 ), (tk , fk ), (tk+1 , fk+1 ) (9.28)

fk+1 = f(tk+1 , p(tk+1 ))

With these new points, we obtain a corrected approximation c(tk+1 ) (best approxima-
∫t
tion) of the integral tk+1
k
f(t, y)dt, which will be calculated as follows:
∫ h
c(tk+1 ) = y(tk ) + f(t, y)dt (9.29)
0
∫ h h
f(t, y)dt = (fk−2 − 5fk−1 + 19fk + 9fk+1 )
0 24
The third process corresponds to the modification. In this part, the aim is to cor-
rect between the prediction process and the correction calculation. This operation allows
the accuracy of the numerical solution to increase. Thus, the modification m(tk+1 ) is
calculated by the following equation:

251
m(tk+1 ) = p(tk+1 ) + (c(tk+1 ) − p(tk+1 )) (9.30)
270
To date, the processes of prediction, correction and modification have been described.
Under these conditions, it is now possible to define the complete Adams–Bashforth
algorithm for the numerical solution of differential equations. The steps refer to the
computation of the three processes in a sequential way prediction → modification →
correction. In summary, the algorithm presents the following calculations:
9.6 Predictor–Corrector Method 193

h
p(tk+1 ) = y(tk ) + (−9fk−3 + 37fk−2 − 59fk−1 + 55fk ) (9.31)
24
251
m(tk+1 ) = p(tk+1 ) + (c(tk+1 ) − p(tk+1 ))
270
h
c(tk+1 ) = y(tk ) + (fk−2 − 5fk−1 + 19fk + 9fk+1 )
24

9.6.1 Functions for Generic Implementation

In this section, the Adams–Bashforth solution method will be implemented generically.

This means that contrary to past implementations, this time, the differential equation to be
solved can be defined. To introduce a differential equation into the code, it is necessary
to explain how to define it.
To define a differential equation in the MATLAB code, it is first necessary to work it
algebraically in such a way that it is defined as follows:

y'(t) = f(t, y) (9.32)

There are two ways to define function 9.32, implicitly or as a function. In this first
exercise, we will do it implicitly, which is the simpler of the two methods.
To define a function implicitly, first define any variable that will represent the function
in the following code. The name of the variable is indistinct, governed only by the MAT-
LAB rules for describing a variable. Then, the @ symbol is added. After the symbol in
parentheses, the variables that participate in the function are written. The order is impor-
tant since the order will be the way the data will be sent when the function is evaluated.
After specifying the variables, the function itself is defined. Assume that the differential
equation to be defined is:

y'(t) = 1 − y(t) (9.33)

194 9 Differential Equations

Its implicit definition in MATLAB would be:

fun=@(t,y) 1-y

This line of code is associated with the variable fun the function 1 − y(t) that repre-
sents the differential equation. Once the function associated with the differential equation
is defined, it is necessary to evaluate it. There are two methods for its evaluation. The
simplest is to simply call it in code by passing it the parameters in the sequence in which
they were defined. Therefore, to obtain the Val evaluation of the function fun at t = 0
and y = 3, it will be written in code:

Val=fun(0,3)

Another method for its evaluation, MatLAB, has the feval function. Its generic
syntax is as follows:

Eval = feval(fun,x1,...,xM)

The function feval evaluates the function fun using as arguments the values ×
1,...,xM. These values must be in the order in which the function was defined
implicitly.
Having explained the use of the feval function and the way to implicitly define a
differential equation in the code, we are now ready to implement the Adams–Bashforth
method in a generic way. The method is coded as a function in Program 9.4.
9.6 Predictor–Corrector Method 195

%Program 9.4

function [t,y]=ode_AB(f,tspan,y0,N)

%función para resolver numéricamente y'(t)=f(t,y(t))

%mediante en método Adams-Bashforth

% donde tspan=[t0,tf] es el intervalo de simulación

% y0 el valor inicial y N el número de pasos

y0=y0(:)';

h=(tspan(2)-tspan(1))/N; tspan0=tspan(1)+[0 3]*h;

for k=1:4, F(k,:)= feval(f,t(k),y(k,:)); end

p= y(4,:); c= y(4,:); KC= 251/270;

h24=h/24; h241=h24*[1 -5 19 9]; h249=h24*[-9 37 -59 55];

for k=4:N

p1 = y(k,:) + h249*F; % Eq.(6.4.8a)

m1 = pk1 + KC*(c-p); % Eq.(6.4.8b)

fk1 = feval(f,t(k+1),m1)];

c1 = y(k,:) + h241*[F(2:4,:); fk1(:)']; % Eq.(6.4.8c)

y(k+1,:) = c1 - KC*(c1-p1); % Eq.(6.4.8d)

p=p1; c=c1;

fk11 = feval(f,t(k+1),y(k+1,:)});

F = [F(2:4,:); fk11(:)'];

end
196 9 Differential Equations

9.7 System of Differential Equations

To date, the numerical solution to a differential equation has been found. However, in
many situations, the ratios of change of variables are related in two different systems
or functions. Such models are called systems of differential equations. A system of
differential equations is defined generically with a set of equations such as:
dy dx
= f (t, x, y), = g(t, x, y) (9.34)
dt dt
with initial values x(0) and y(0). To numerically solve this system of differential equa-
tions follows the same principle as with a simple differential equation. Only in the case
of a system is a numerical method performed for each of the derivatives of the variables.
Perhaps the main difference with the case of a single differential equation results in the
fact that the updating of the solution values is interleaved with the values of the variables
previously calculated. This is important because it must be taken into account that the
values of the variables are related by the system as a whole.
Thus, its solution through the Runge–Kutta method will be determined by the
calculations defined in the following table:
dy dx
= f (t, x, y) = g(t, x, y)
dt dt
k1 = h · f (t,
( x, y) ) q1 = h · g(t,
( x, y) )
1 1 1 1 1 1
k2 = h · f t + h, x + q1 , y + k1 q2 = h · g t + h, x + q1 , y + k1
( 2 2 2 ) ( 2 2 2 )
1 1 1 1 1 1
k3 = h · f t + h, x + q2 , y + k2 q3 = h · g t + h, x + q2 , y + k2
2 2 2 2 2 2
k4 = h · f (t + h, x + q3 , y + k3 ) q4 = h · g(t + h, x + q3 , y + k3 )

1 1
y(k + 1) = y(k) + (k1 + 2k2 + 2k3 + k4 ) x(k + 1) = x(k) + (q1 + 2q2 + 2q3 + q4 )
6 6
(9.35)

From this calculation definition, it is clear that the updating of the values ki and qi are
handled interleaved for the two variables x and y.
As an example of implementation, the following system of differential equations will
be assumed.
dy dx
= ax − by, = −ax (9.36)
dt dt
With initial conditions x(0) = 100, y(0) = 0,, we would have as analytical solutions
the following functions:

x(t) = x(0) · e−at , (9.37)

9.7 System of Differential Equations 197

a ( )
y(t) = x(0) e−at − ebt
b−a
Assuming that a = 0.1, b = 0.2, x(0) = 100, y(0) = 0 and t f = 10, Program 9.5 finds
the numerical solution according to the Runge–Kutta method for the system of differential
equations defined in 9.36.

'Rtqitco";07

clear all

tf=10;

a=0.1;

b=0.2;

h=0.25;

t=0:h:tf;

x0=100;

y0=0;

f=@(t,x,y) -a*x;

g=@(t,x,y) a*x-b*y;

x=zeros(length(t),1);

y=zeros(length(t),1);

x(1)=x0; y(1)=y0;

for i=1:length(t)-1

k1=h*f(t(i),x(i),y(i));

q1=h*g(t(i),x(i),y(i));
198 9 Differential Equations

k2=h*f(t(i)+h/2,x(i)+k1/2,y(i)+q1/2);

q2=h*g(t(i)+h/2,x(i)+k1/2,y(i)+q1/2);

k3=h*f(t(i)+h/2,x(i)+k2/2,y(i)+q2/2);

q3=h*g(t(i)+h/2,x(i)+k2/2,y(i)+q2/2);

k4=h*f(t(i)+h,x(i)+k3,y(i)+q3);

q4=h*g(t(i)+h,x(i)+k3,y(i)+q3);

x(i+1)=x(i)+(k1+2*k2+2*k3+k4)/6;

y(i+1)=y(i)+(q1+2*q2+2*q3+q4)/6;

end

plot(t,x)

hold on

plot(t,y)

Figure 9.5 shows the evolution of the two variables when solved by the Runge–
Kutta method. The method is effective and allows us to numerically solve the system
of differential equations.

Fig. 9.5 Solution obtained by the Runge–Kutta method for a system of differential equations with
h = 0.25.
9.8 Second Degree Differential Equations 199

9.8 Second Degree Differential Equations

When a differential equation presents a double derivative, it is said that the differen-
tial equation is of second order. A second-degree equation is defined generically by the
following expression:

y '' (t) + ay'(t) + by(t) = s(t) (9.38)

In differential Eq. 9.38, the elements a and b are constants, while s(t) is a function
in terms of t, y(t) or y'(t). Finding the solution of a second-order differential equation
requires the definition of two initial conditions involving the dependent variable y(0) = y0
and its first derivative y'(0) = y0 '. If s(t) is independent of y(t), the differential Eq. 9.38
is linear. On the other hand, if s(t) is defined in terms of either y(t) or y'(t), Eq. 9.38
is nonlinear. In this subsection, we will explain how to solve this type of differential
equation numerically using the Runge–Kutta method.
The first step in solving a second-degree differential equation is to decompose the
equation into a system of first-degree equations, as discussed in Sect. 9.7. For this, an
auxiliary variable x(t) is considered that decomposes the equation into two. With this
decomposition, the system of differential equations has the following form:

x(t) = y'(t) (9.39)

x'(t) + ax(t) + by(t) = s(t)

With this transformation (the inclusion of an auxiliary variable), the initial condition
comes from the second condition. In this way, x(0) = y0 '. Under these circumstances, the
set of differential equations defined in 9.39 can be expressed under the following system:

y'(t) = x(t) (9.40)

x'(t) = s(t) − ax(t) − by(t) = f (t, x, y)

Thus, its solution through the Runge–Kutta method will be determined by the
calculations defined in the following table:
200 9 Differential Equations

dy dx
= f (t, x, y) = g(t, x, y)
dt dt
k1 = h · f (t,
( x, y) ) q1 = h · g(t,
( x, y) )
1 1 1 1 1 1
k2 = h · f t + h, x + q1 , y + k1 q2 = h · g t + h, x + q1 , y + k1
( 2 2 2 ) ( 2 2 2 )
1 1 1 1 1 1
k3 = h · f t + h, x + q2 , y + k2 q3 = h · g t + h, x + q2 , y + k2
2 2 2 2 2 2
k4 = h · f (t + h, x + q3 , y + k3 ) q4 = h · g(t + h, x + q3 , y + k3 )

1 1
y(k + 1) = y(k) + (k1 + 2k2 + 2k3 + k4 ) x(k + 1) = x(k) + (q1 + 2q2 + 2q3 + q4 )
6 6
(9.41)

As an example of implementation, the following second degree equation will be

assumed:
d2 y dy
2
+ 2g + ω02 y = 0 (9.42)
dt dt
This differential equation models a damped harmonic oscillator where its initial con-
ditions establish its initial position y(0) and its initial velocity y'(0). As a first step,
this second-degree differential equation is decomposed into a set of two first-degree
differential equations. From this decomposition, the set is defined as:

y'(t) = x(t)

x'(t) = f (t, x, y) = −2gx(t) − ω02 y (9.43)

Assuming ω0 = 2, g = 0.5, x(0) = 100, y(0) = 1.5, y'(0) = 0 and t f = 8, Pro-

gram 9.6 finds the numerical solution according to the Runge–Kutta method for the
second-degree differential equation defined in 9.42 with a step h of 0.08. Program 9.6
finds the numerical solution according to the Runge–Kutta method for the second-degree
differential equation defined in 9.42.
9.8 Second Degree Differential Equations 201

%Program 9.6

clear all

h=0.08;

tf=8;

t=0:h:tf;

w0=2;

g=0.5;

f=@(t,y,x) -2*g*x-w0*w0*y;

y=zeros(length(t),1);

x=zeros(length(t),1);

y(1)=1.5; x(1)=0;

for i=1:length(t)-1
k1=h*x(i);

q1=h*f(t(i),y(i),x(i));

k2=h*(x(i)+q1/2);

q2=h*f(t(i)+h/2,y(i)+k1/2,x(i)+q1/2);

k3=h*(x(i)+q2/2);

q3=h*f(t(i)+h/2,y(i)+k2/2,x(i)+q2/2);

k4=h*(x(i)+q3);

q4=h*f(t(i)+h,y(i)+k3,x(i)+q3);

y(i+1)=y(i)+(k1+2*k2+2*k3+k4)/6;

x(i+1)=x(i)+(q1+2*q2+2*q3+q4)/6;

end

plot(t,y)

Figure 9.6 shows the evolution of the variable y(t) when solved by the Runge–
Kutta method. The method is effective and allows us to numerically solve the system
of differential equations.
202 9 Differential Equations

Fig. 9.6 Solution obtained by the Runge–Kutta method for the second-degree differential equation
with h = 0.08

MatLAB has already built in functions for the numerical solution of differential
equations. Among them, the most popular one is the ode45 function. This function
implements the Runge–Kutta method, which was already explained in Sect. 9.4. Its syntax
of use is defined by:

[t,y] = ode45(odefun,tspan,y0)

where tspan = [t0 tf] is a vector that indicates the start time t0 and end time
t f in which the equation or system of differential equations y' = f (t, y) is integrated.
The element y0 represents the initial conditions that will be considered in the numerical
solution.
odefun represents the function where the equation or system of differential equations
is defined. This function can be defined in implicit form, as already seen in this chapter,
or it can be defined as an external function. The definition of a function in external form
usually occurs in two scenarios:
9.8 Second Degree Differential Equations 203

When the differential equation is very complex from the point of view of the variables
it relates.
When it is a system of differential equations.

To demonstrate the use of MATLAB’s ode45 function, the following second-degree

differential equation will be solved:
( )
y '' − μ 1 − y 2 y' + y = 0 (9.44)

With initial conditions:

' '
y(0) = 2 ; y (0) = 0 (9.45)

To solve this differential equation, as already explained, the first step is to decom-
pose the equation as a system of two first-degree differential equations. Performing this
decomposition, the system would be defined by the following formulation:

y1 ' = y2 (9.46)

( )
y2 ' = μ 1 − y12 y1 ' + y1

This is the set of differential equations to be defined as an external function. Consid-

ering the constant μ = 1, the function is called sist_seg.m and is shown in Program
9.7.

'Rtqitco";09"

function dydt = sis_seg(t,y)

dydt = [y(2); (1-y(1)^2)*y(2)-y(1)];

This function encodes the System of Differential Equations. The elements y1 and y2
are the inputs of the system, while the vector dydt represents the change in the associated
variables in each iteration.
Once the function is defined, the MATLAB function ode45 is used to solve the system.
For this purpose, a simulation period from 0 to 20 is considered. This is done through the
following command:

[t,y] = ode45(‘sist_seg’,[0 20],[2 0]);

204 9 Differential Equations

Fig. 9.7 Solution obtained by the Runge–Kutta method for the second-degree differential equation
with the Ode45 function

Plotting the variable y yields a graph of the numerical solution for this differential
equation (Fig. 9.7).

9.9 Differential Equations with Boundary Value

Up to this point, we have considered the numerical solution of differential equations for
which their initial values are known. However, many engineering applications require
that their boundary values be specified. That is, not only the initial value but also the final
value of a function that is related to its rate of change is needed. A simple example of
these problems can be visualized through the following case defined by a second-degree
differential equation:
( )
d2 y dy
= f t, y, (9.47)
dt 2 dt

This differential equation is solved on an interval t\in[a,b] with the following general
boundary conditions:
9.9 Differential Equations with Boundary Value 205

dy(a)
α1 y(a) + β1 = γ1 (9.48)
dt
dy(b)
α2 y(b) + β1 = γ2
dt
Under these conditions, the solution of Eq. 9.47 over a specific interval satisfies the
initial and final values described in the conditions defined in 9.48. Figure 9.8 graphically
exhibits a generic representation of a boundary condition problem. In the graph, the points
y(a) and y(b) are specified as boundary conditions, while the dynamic behavior between
those points is described by a second-order differential equation.
For the definition of the boundary problems, information is required at the initial time
t = a and at the final time t = b. This is the fundamental difference from the initial
condition methods (t = a) discussed in the previous sections. The methods for solving
boundary problems are modifications of the single condition methods.
To solve a problem of boundary conditions of a generic differential equation such as
9.37, it is necessary to define an operating interval t ∈ [a, b] and the values that the
function will have at the extremes of the interval, such that y(a) = α and y(b) = β.
Numerical methods for solving second-degree differential equations require knowledge
of the initial conditions y(a) and y'(a). To solve a boundary condition problem, the value
of the condition y'(a) is not fixed. This value is modified with the objective that the
numerical solution of the differential equation obtains its final value y(b) = β. This

Fig. 9.8 Visualization of the generic solution of a boundary condition problem

206 9 Differential Equations

'
Fig. 9.9 Solutions to the boundary problem with y(a) = α and y (a) = Vi . In the figure, two values
V1 and V2 are used to illustrate the behavior of the solutions and their inaccuracy in trying to reach
the required final value y(b) = β

method is known as the shooting method. Under this method, iteratively, the value of
y'(a) = Vi is modified until the final solution of the differential equation y(b) acquires
the value of β. This process is shown in Fig. 9.6. The behavior of the solutions in the
figure for the two values V1 and V2 suggests that it is possible to use a bisection method
to find the optimal value Vi that allows the final solution to be y(b) = β. The case of
the solution found with y'(a) = V1 makes it clear that its final value y(b) is too low to
satisfy the value of β. On the other hand, the solution that corresponds to y'(a) = V1
yields a final solution with a value much higher than β. Under these conditions, it can
be assumed that if the value of Vi is modified with values between [V1 , V2 ], the optimal
value of Vi can be found so that the solution y(b) reaches the value of β (Fig. 9.9).
'
The firing method allows a sweep of the condition y (a) = Vi in such a way that the
final boundary value of the solution to the differential equation is found. The algorithm
has the following steps:
9.9 Differential Equations with Boundary Value 207

1 The differential equation is solved using any method of numerical solution of differential
equations with the initial conditions y(a) = α and y'(a) = Vi
2 The value of the solution at the end point y(b) is evaluated and compared with the required
value β
3 The value of Vi is adjusted up or down depending on the strategy. This process is carried out
until a good approximation of the final value is reached and y(b) = β

To exemplify the implementation and use of the shooting approach, we consider

solving the following differential equation:
( )
y '' + x 2 − sint y' − (cos2 t)y = 5 (9.49)

where t ∈ [0, 1] considering the following boundary conditions:

y(0) = 3 (9.50)

y(1) = 5

According to the firing method, the first step is to write the second-degree differential
equation described in 9.49 as a system of first-degree differential equations, such that the
system is defined as follows:

y1 ' = y2 (9.51)

( )
y2 ' = −(t 2 − sint)y2 + cos2 t y1 + 5

y1 (0) = 3, y1 (1) = 5

The objective is to use the Vi value of the initial condition y1 '(0) to find a dynamic
behavior that allows that at the end, the value of the function y(1) is 5. The algorithm for
finding the optimal Vi value is the firing method.
As indicated by the triggering algorithm, the value of Vi must be modified at each inter-
action so that at each new value, it is analyzed whether the dynamic behavior complies
with the boundary conditions.
Program 9.8 presents the implementation of the firing method to solve the boundary
conditions problem. Program 9.9 shows the explicit definition of the system function
describing the set of differential equations formulated in 9.51.
208 9 Differential Equations

In Program 9.8, it is clear how starting from an initial value of Vi = −3, this value is
checked to see if the final value of the solution is greater or less than the value specified
by the boundary condition β = 5. If the value is less than 5, the value of Vi is increased
to favor the final value of the function to increase. On the other hand, if the value is
greater than 5, the value of Vi is decreased so that the final values of the function tend to
decrease.
Figure 9.10 shows the values obtained by the shooting method. As shown in Fig. 9.10,
the method is effective in finding the final value of the differential equation.

Fig. 9.10 Solution obtained by the firing method to obtain the solution to the boundary problem
described in 9.49 and 9.50
9.9 Differential Equations with Boundary Value 209

"

'Rtqitco";0:"

xspan=[0 1];

V=-3; dV=0.5;

for j=1:100

y=[3 V];

[t,ysol]=ode45('sistema', xspan,y);

if abs(ysol(end,2)-5)<10^(-6)

break

end

if ysol(end,2)<5

V=V+dV;

else

V=V-dV;

dV=dV/2;

end

end

plot(ysol)

'Rtqitco";0;"

function sal=sistema(t,y)

sal=[y(2); -(t^2-sin(t))*y(2) + (cos(t)^2)*y(1) + 5];

210 9 Differential Equations

9.10 Graphical Methods

As seen in this chapter, many differential equations cannot be solved in a closed or analyt-
ical way, which is why a large number of algorithms have been developed to solve them
numerically. In this section, we will highlight alternative methods to the numerical ones
for the analysis of the solutions of differential equations. These are graphical methods that
allow us to understand the nature of the solutions by analyzing various points of action.
The graphical method follows the fundamental concept that a derivative is considered
the slope of a straight line that is tangent to a curve given a specific point. If a differential
equation is considered as:

y' = f (x, y) (9.52)

which has an initial condition y(x0 ) = y0 . For any point (x, y), it is possible to plot a
small line whose slope is equal to f (x, y). This graphical representation is known as the
potential field of Eq. 9.39. Under this method, starting from the initial point (x0 , y0 ), the
solution curve can be constructed by extending the initial segment in such a way that the
tangent of the solution is parallel in the direction of the potential field through which the
curve passes.
In the graphical method, it is common to first draw the lines with a constant slope
(isolines) where f (x, y) = C. To illustrate this technique, we consider obtaining the
solutions of Eq. 9.53 by the graphical method:
dx
= x − t2 (9.53)
dt
where its exact solution is determined by:

x(t) = Cet + t 2 + 2t + 2 (9.54)

where C is an arbitrary constant. Program 9.10 shows the code for finding the potential
field by the graphical method.
9.10 Graphical Methods 211

'Rtqitco";032"

clear

[t,x] = meshgrid(-2:0.2:3,-1:0.2:2);

slope = x - t.*t;

length = sqrt(1 + slope .* slope);

quiver(t,x,1./length,slope./length,0.5)

quiver(t,x,1./length,slope./length,0.5)

axis equal tight

hold on

tt =-2:0.2:3;

for ind =-10:1:10

xreal = ind*exp(tt) + tt.*tt + 2*tt + 2;

plot(tt,xreal)

hold on

end

xlim([-2 3])

ylim([-1 2])

Figure 9.11 shows the solution of the graphical method for each possible point (x, y)
within the limits x [−2,3] and y [−1,2]. In the graph, one can also visualize some
solutions for fixed values ranging from [−10, 10].
212 9 Differential Equations

Fig. 9.11 Graphical method for the analysis of differential equations

References

1. Won Y. Yang, Wenwu Cao, Jaekwon Kim, Kyung W. Park, Ho-Hyun Park, Jingon Joung, Jong-
Suk Ro, Han L. Lee, Cheol-Ho Hong, Applied Numerical Methods Using MATLAB, 2020 John
Wiley & Sons, Inc.
2. Burden, Richard L. and J. Douglas Faires, Numerical Analysis, 7th ed., Brooks/Cole
3. Fausett, Laurene V., Applied Numerical Analysis using MATLAB, Prentice Hall, Upper Saddle
River, New Jersey, 1999.
4. Lindfield, George R. and John E.T. Penny, Numerical Methods using MATLAB, 8th ed., Prentice
Hall, Inc., Upper Saddle River, New Jersey, 2000.
5. Mathews, John H. and Kurtis D. Fink, Numerical Methods Using MATLAB, Prentice Hall, Inc.,
Upper Saddle River, New Jersey, 1999.
6. Schilling, Robert J. and Sandra L. Harris, Applied Numerical Methods for Engineers using
MATLAB and C, Brooks/Cole Publishing Company, Pacific Grove, CA, 2000.
7. Recktenwald, Gerald W., Numerical Methods with MATLAB, Prentice Hall, Inc., Upper Saddle
River, New Jersey, 2000.
8. Stoer, J. and R. Bulirsch, Introduction to Numerical Analysis, Springer-Verlag, New York, 1980