Wichita State University Logo

Math 451: Computational Mathematics

2.2 Structured Programming

2.2.1 Local Functions

Local functions are defined in an M-file just like scripts. There are two kinds of local functions

  • The first local function in an M-file has the same function name as the filename of the file in which it resides.

Like scripts these local functions are named procedures . The function is invoked using the filename of the M-file. The scope of a named procedure is the MATLAB search path. Any script, function, or a statement executed from the command window can invoke this local function.

  • Subsequent local functions in an M-file have function names that may not match the filename or the name of any other local function defined in the M-file.

These local functions are invoked using their function name. These local functions have file scope . That is they may be invoked by any MATLAB statement that resides in the same file as the local function .

2.2.2 Structured Programming

Each of these functions types are structures for encapsulating a single algorithm. The use of the function structure can be very powerful. It allows the programmer to implement a complicated algorithm, test it and verify its accuracy. Then any MATLAB statement within the scope of the function may utilize the algorithm simply by invoking the function and passing appropriate inputs.

  1. Functions create abstractions that communicate an algorithm through careful choice of the name of the function.
  2. Functions group a function block of complex computations into an action that can be used to communicate the flow of control while hiding the details of the implementation.
  3. Functions create a separate workspace to implement algorithms that can be reused by many seemingly different applications.

The structured programming paradigm utilizes language keywords to create logical blocks of code within a function as well.

  • A programming structure may implement a decision element of an algorithm to select one of several code paths for execution.

  • A programming structure may encapsulate a block of several MATLAB statements that perform complex computations to achieve a single result for an algorithm.

  • A programming structure may execute a code block repeatedly to implement an algorithm.

2.2.3 Conditional Structures

Conditional structures allow the state of the workspace to determine the flow of control in a function. The basic arithmetic conditional statement executes an associated code block only if the expression or condition evaluates to true . If the the condition evaluates to false , then the conditional structure passes control to the next code block defined by the structure or after the end of the structure.

Basic Conditional Flowchart

Thus, the code block in the body of the conditional structure is executed only when the value of the expression in the condition statement resolves to true when the statement executes. This means that one will not know whether the code block of the condition structure will execute when writing the function. The state of the variables in the workspace will determine the value of the expression at runtime . When the condition resolves to true during execution we say that the condition is satisfied and the code block is executed. When the condition resolves to false during execution we say that the condition was not satisfied and the code block is skipped.

2.2.4 A Simple Conditional Structure

Our example function quadraticformula contains several conditional structures.

The first basic conditional structure uses the MATLAB keyword if followed by a MATLAB logical expression . 

if nargin == 1
    b = a(2);
    c = a(3);
    a = a(1);
end

The conditional statement 

if <MATLAB expression>

declares the beginning of the conditional structure and it is paired with an 

end

statement. The code block between these two statements make up the body , conditional block or if block of the conditional structure. The MATLAB expression in the if statement must result in a logical value. Logical expressions are usually the result from a logical operators or relational operators . In our example the relational operator == determines the equality of two values. The MATLAB expression 

>> nargin == 1

determines the equality of the value returned by the function nargin and the literal value 1. If we use function notation instead of the operator, this is equivalent to the MATLAB statement 

>> ans = eq(nargin(),1)

The MATLAB if statement requires a MATLAB logical expression like this one and determines the next MATLAB statement to execute based on the value of the expression. The conditional block of code will be executed when the value of the expression is true . Most MATLAB statements may be excuted in the conditional block. However functions may not be defined in another programming structure such as a conditional structure.

2.2.5 A More Complicated Example

The statements in the conditional block may return control of execution to the invoking script or function. Consider the conditional structure in our example 

%...
    if ~all(isnumeric([a b c]))
        error('%s: %s. %s', ...
            'Error', 'Invalid expression', ...
            'The input values must be real scalars.')
    end

The error function raises and exception . This will immediately release flow of control to the MATLAB Interpreter. Normally, after the last MATLAB statement in the conditional block executes, the end statement is executed. The end statement terminates execution of the conditional structure and the flow of execution proceeds to the MATLAB statement following the end statement. That is not the case in this example.

In this example, the MATLAB interpreter executes the most recent exeption control block registered with the interpreter. In this case control will be returned to the MATLAB command window. We will study how to capture exceptions with our own exception control block in subsequent chapters.

For now raising a error exception by calling error will terminate the execution of quadraticformula and execution will resume with the invoking process.

2.2.6 Conditional Structures with Multiple Code Blocks

If the logical expression in an if statement is empty or false , then the if statement jumps to either the next conditional statement or the end statement. A conditional structure may contain more than one conditional statement ! A MATLAB conditional statement may utilize an elseif statement or an else statement. The elseif statement also requires a logical expression to determine the flow of execution of the conditional structure. An else statement has no expression and always executes its code block.

In quadraticformula, the second conditional structure defines multiple conditional statements. 

    % Inputs should be numbers only
    if ~all(isnumeric([a b c]))
        error('%s: %s. %s', ...
              'Error', 'Invalid expression', ...
              'The input values must be real scalars.')

    % Inputs should be real numbers (no complex polynomials!)
    elseif ~all(isreal([a b c]))
        error('%s: %s. %s', ...
              'Error', 'Invalid expression', ...
              'The input values must be real valued.')

    % Compute the two roots
    else

        if abs(a) > eps
            % Compute the discriminant once because it is required to
            % determine both roots.
            discriminant = sqrt(b^2 - 4*a*c);

            % Calculate the two roots and assign the values to the vector x
            x = (-b + [ discriminant; -discriminant ])./(2*a);
        end
    end

The flowchart for our algorithm

Nested Conditional Flowchart

emphasizes that the algorithm consists of three mutually exclusive code blocks.

  1. A code block that executes when the argument list contains non-numeric values

  2. A code block that executes when the argument list contains complex values

  3. A code block that computes the roots of a valid real-valued quadratic polynomial

The first two code blocks invoke the MATLAB function error and terminate execution of the function. The third code block is and else block. The else block executes only if the second conditions of the if block and elseif block are not satisfied.

2.2.7 Numerical Instability

Our quadraticformula function appears to be working fine, but there is an issue if we ask it to perform particular types of computations. For instance, suppose we ask it to solve the equation

$$ x^2 - \left(10^8 + 10^{-8}\right)x + 1 = 0. $$

In [17]: 
format long
quadraticPolynomial = [1 -(1e8 + 1e-8) 1];
r = quadraticformula(quadraticPolynomial)

r =

   1.0e+08 *

   1.000000000000000
   0.000000000000000

One can see that the numerical results must be incorrect because quadraticformula is claiming that one of the roots is 0. Plugging 0 into our equation should give 1, not 0, which is evidence of an issue. Let's figure out what the roots should be using exact (by hand) arithmetic :

$$ \begin{aligned} r &= \dfrac{\left(10^8 + 10^{-8}\right) \pm \sqrt{\left(10^8 + 10^{-8}\right)^2 - 4(1)(1)}}{2(1)} \\ &= \dfrac{\left(10^8 + 10^{-8}\right) \pm \sqrt{\left(10^{16} + 2 + 10^{-16}\right) - 4}}{2} \\ &= \dfrac{\left(10^8 + 10^{-8}\right) \pm \sqrt{\left(10^{16} - 2 + 10^{-16}\right)}}{2} \\ &= \dfrac{\left(10^8 + 10^{-8}\right) \pm \sqrt{\left(10^{8} - 10^{-8}\right)^2}}{2} \\ &= \dfrac{\left(10^8 + 10^{-8}\right) \pm \left(10^{8} - 10^{-8}\right)}{2} \end{aligned} $$
From here, we can see that the two roots of our quadratic equation are $10^8$ and $10^{-8}$. If we plug either value into the original quadratic equation, we get zero.

$$ (10^8)^2 - \left(10^8 + 10^{-8}\right)(10^8) + 1 = 0, \qquad\qquad (10^{-8})^2 - \left(10^8 + 10^{-8}\right)(10^{-8}) + 1 = 0. $$
These computations have an "infinite" amount of precision, but we can see that if we use these values in a numerical computation the result fails to be exactly zero, or even remotely close to it in one case.

In [21]: 
polyval(quadraticPolynomial,r)

ans =

  -0.490116119384766
  -0.490116119384766

It is caused by the fact that we are performing a addition operation involving terms of the form $10^8+10^{−8}$. This number has more than 16 significant digits, beyond the double precision floating point limit. When MATLAB adds these numbers, it forgets one of the significant figures because they are spread too many place values apart. The built-in function polyval is not immune to this.

$$ \begin{align*} \ \ 100000000&.0000000{\color{red}0} \\ +\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 0&.00000001 \\ \hline \\ \ \ 100000000&.0000000{\color{red}1} \end{align*} $$
This computation yields a result that requires 17 digits of accuracy to represent! When your computer processor performs this addition the result becomes

$$ \begin{align*} \ \ 100000000&.0000000 \\ +\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 0&.00000001 \\ \hline \\ \ \ 100000000&.0000000 \end{align*} $$
In the finite arithmetic of double precision floating point mathematics,

$$ 10^8 + 10^{-8} = 10^8 $$
This phenomenon is known as catastrophic cancellation. An addition of numbers that stretches beyond 16 decimal digits to compute are truncated to 16 decimal digits.

2.2.8 Numerically Stable Computations

While it is not possible to solve every problem, there is a way to make our root finding function quadraticformula more numerically stable , meaning that it loses less precision and introduces smaller errors into its calculations. For our example problem

$$ x^2 - \left(10^8 + 10^{-8}\right)x + 1 = 0, $$
our algorithm (the standard quadratic formula) can only find one of the roots. There is a modification that we can make which will allow it to find both. First, we determine the root with the largest absolute value using

$$ r_1 = \dfrac{-b - \text{sign}(b)\sqrt{b^2 - 4ac}}{2a}, $$
where

$$ \text{sign}(y) = \left\{\begin{matrix} \ 1 & y \ge 0\, \\ -1 & y \lt 0. \end{matrix}\right.$$
The second root is found with the identity

$$ r_1 r_2 = \frac{c}{a}. $$

2.2.9 Homework 2.2

Create a flowchart for this new computation of the roots of a quadratic polynomial. Include the parameter checking above in your algorithm. You algorithm should return both roots in a $1\times 2$ array.

Creative Commons Logo - White


Your use of this self-initiated mediated course material is subject to our Creative Commons License .


Creative Commons Attribution-NonCommercial-ShareAlike 4.0

Creative Commons Logo - Black
Attribution
You must give appropriate credit, provide a link to the license, and indicate if changes were made. You may do so in any reasonable manner, but not in any way that suggests the licensor endorses you or your use.

Creative Commons Logo - Black
Noncommercial
You may not use the material for commercial purposes.

Creative Commons Logo - Black
Share Alike
You are free to share, copy and redistribute the material in any medium or format. If you adapt, remix, transform, or build upon the material, you must distribute your contributions under the same license as the original.