MATLAB Applications Part 1: Numerical Methods

Introduction

• modelling: creating representations of a scientific phenomena
• scientific computing: using the capabilities of a computer to analyze these models and solve difficult problems

1. Intro to Numerical Analysis

Around 2600 years ago, Pythagoras and the Greeks believed that the universe could be expressed entirely through natural numbers and ratios between them. Some time later, their world view was disturbed when they discovered that is irrational using the Greek ladder technique. And yet, a few millenia before even this, the Babylonians and Indians had methods to approximate .

The Babylonians were sophisticated mathemeticians, known for developing their own numerical system. A Babylonian tablet, YBC (Yale Babylonian Collection) 7289, is one example of just how impressive their work was. In base-60 (can you think of anything in our lives today that uses the number 60?), one diagonal contains the numbers . After converting to decimal, this refers to the following fraction, which is incredibly close to :

Why did they record ? Well, notice how the square also has a side labeled 30. Multiplying these two values together gets the value along the other diagonal, . Given a side length of 30, they used the approximation of to find the diagonal of the square.

Note: In Babylonian mathematics, division by 2 was equivalent to multiplication by 30, since 2 and 30 are reciprocal values in their number system . This allowed them to scale values depending on what they were solving for.

a is the approximation to

B is a "correction" term such that

improves the estimate

Let's try evaluating this for :

Let's try evaluating this for :

is your current guess for

is then your improved guess for

If we plug in , we get:

In other words, the Babylonian method is Newton's method for a function in the form of which has the root .

1.1 Numbers and Rounding: Unsigned and Signed Integers

Types of Numbers	Definition	Description
Natural numbers		The set of all positive whole numbers
Unsigned integers		An extension of the natural numbers that includes zero
Integers		An extension of the unsigned integers that includes negative whole numbers
Rational numbers		An extension of the integers to include fractions
Real numbers		An extension of the rational numbers to include irrational numbers
Complex numbers	ℂ =	An extension of the real numbers to include imaginary numbers

The 8-bit binary number is equivalent to the 4-bit binary number 1001. There are different ways to represent binary numbers depending on the number of bits.

Try representing the decimal number 16 through binary. That's right, you cannot. 4-bit binary numbers can only represent a total of 16 () possible values, including 0. These 4-bit binary numbers are also known as nibbles.

8-bit binary numbers can represent a total of 256 () values, including 0.

in binary

in binary

in binary

in binary

in binary

in binary

Number of Digits d	Form	Min Value	Max Value
1		0
2
3
n

Letwhere n is an integer and .

Then there is some such that

This gives us n, the number of decimal digits needed to represent .

Let's test this in MATLAB. How many digits are needed to represent ?

Number of Digits b (bits)	Form	Min Value	Max Value
1		0
2
3
n

Letand .

Then there is some such that

This gives us n, the number of bits needed to represent .

Let with . We have the expression:

This is an increasing function with where z is the domain of the function.

Let . We can simplify:

Keep in mind that . And the domain of . Thus,

We have proved the propisition that if , then .

where

x = -10:0.01:10;

y = floor(x);

figure

scatter(x,y)

xlabel('Actual Value')

ylabel('Floored Value')

title('Floored Value Versus Actual Value')

Let's look at the floor and rem functions. If take the floor of , we get the following:

The actual value of is , but floor returns 12. This is because it rounds down.

V = -10:1:10;

M = [mod(V(:),3),rem(V(:),3)];

figure

plot(V(:),M)

xlabel('Input')

ylabel('Output')

title('Difference Between Rem and Mod')

legend('mod','rem','location','best')

Aspect	mod(a,m)	rem(a,m)
Purpose	returns the remainder after division	returns the remainder after division
Rounding	floor(a/m) (rounds toward )	fix(a/m) (rounds toward 0)
Output	Always 0 or same sign as a	Always 0 or same sign as m
a > 0	Returns positive number	Returns negative number
m = 0	Returns a	Returns 0

where

Given the condition and the fact we want an integer quotient , should we be using floor division or fix? We should be using floor. Thus, because . So we now have the following:

5 fits into 32 exactly 6 times, meaning that is accounted for.

We can see from this that there is a remainder of .

We can also multiply the fractional leftover by the denominator (divisor) 5 to get 2.

Thus, .

27 can be represented by bits b:

	= decimal
	1
	2
	4
	8
	16
	32
	64
	128
	256
	512
	1024
	2048
	4096
	8192
	16384
	32768
	65536

1	1	0	1	1

Thus, 27 in the binary form of is 11011.

149 can be represented by bits b:

1	0	0	1	0	1	0	1

Thus, 149 in the binary form of is 10010101.

149 can be represented by bits b:

1	0	0	1	0	1	1	0	0	0	1	1	1	0	1

Thus, 19229 in the binary form of is 100101100011101.

In , notice how we go from 1-bit integers to a 2-bit integer.

Let's try adding in binary:

Let's try adding in binary:

Notice how ex1 yields the expected . But ex2 gives . The number 300 in binary cannot be represented in 8-bits. This is called overflow.

And then add:

Addition and subtraction in two's complement is pretty simple. Here's what it looks like for a positive integer and a negative integer. We are adding .

Let's try adding two integers in two's complement that result in a sum of 0. We are adding .

Overflow also exists when representation integers via two's complement. Let's say we add up .

91 in binary is 01011011. Flip the bits and then add to get 10100101. We now have .

113 in binary is 01110001. Flip the bits and then add to get 10001111. We now have .

in 16-bit is 1111111110100101.

Now, let's combine it all to evaluate .

Notice how we can go all the way to but only go up to 127. We can represent one more negative number than positive because there there is not enough space to represent the positive counterpart of the most negative number.

Decimal	Binary
	100
	101
	110
	111
0	000
1	001
2	010
3	011

The allowed range is . 100 is deserved for . We cannot go beyond 3 mathematically. Instead, we'd wrap around back to . You can think of two's complement system for a certain number of bits like a circle.

Here, there are a certain number of combinations .

If you are performing a sequence of arithmethic, it's okay if the intermediate operations result in overflow as long as the final result is within the bounds. For instance, let's say we want to do in this system.

For the first operation , there is overflow because .

For the second operation , we "revert" the overflow and get an answer .

1.2 What Are Fixed-Point and Floating-Point Numbers?

It's a good time to talk about how fractions work in these number systems. For any , we have where and represent the numerator and denominator respectively. Thus .

For instance, the fraction can be represented by a 4-digit decimal number:

Let be nonnegative integers. Then, we can write in afinite decimal expansion

Each digit is for all . The total number of decimal digits is thus

Let's look at an example for .

The integer part of our number 8142 is expressed using decimal digits therefore . We can see that decimal digits after the decimal point. All in all,

We cannot write finite decimal expansions for certain fractions because of repeating decimals. for instance is written as . For this, we'd have a different form where

Let's try this out both ways. We have .

Since has 7 decimal digits in the fractional part and we've already moved over 4 digits, let's shift over another 3 digits:

Now we can do to get where the fractional parts cancel out.

More in line with our general notation, we can also express the number as a summation.

We can rewrite the sum as such, where we factor out :

The inside can be represented as a geometric series of general form . Let and :

We can see that . We can then say that

for a sequence of repeating digits

Let I be the integer part , A be the non-repeating part after the decimal point, and B be the repeating decimal part. We now have

Trying this out for our previous number , we can see that , with , and with .

The above is for infinite repeating decimals, but what happens if we have a finite decimal? We should be able to distinguish between the two. We can say that the decimal expansion of is finite if and only if its fractional representation can be expressed with a denominator in the form of where j and k are nonnegative integers.

If we expand , we we will eventually arrive at the following:

In other words, a finite decimal can be written as and since , we can conclude that the fraction representation of a finite decimal must have a denominator that is a product of twos and fives.

Binary representations are what the computer deals in, so let's now work with binary fractions. First, we will define the general form of a finite binary expansion. Let be nonnegative integers. Then can be written in the form

As this is a binary representation we have each bit where . Similarly to before, we know that the total number of bits .

Here's an example dealing with the fraction . It can be written as . In binary form, we get . In binary, recall that decimal 3 is equivalent to binary 11.

For another example, we have . This can be rewritten as In binary form, we get . Once again, recall that decimal 19 is equivalent to binary 10011.

Let's try this out for an example where .

Try plugging in the example we did above with .

where .

Fixed-point arithmetic is essentially the same as integer arithmetic, but with a scale factor. For instance, in Q4.4 format (4 bits for integer, 4 bits for fraction), our scale factor is . We would multiply the stored integers by to get our actual number. We would do the reverse to go from the actual number to the stored integer.

Actual Number	Scaling	Integer Stored By Fixed-Point
		52
		88
		16

• Value v refers to the number assigned to the data type of the constructor.
• Signedness s refers to whether or not the number can be negative. It stores only logical 1 or 0.
• Word length w refers to the total number of bits used to store the number where w . This means you are limited to a maximum of 16 bits to index individual bits.
• Fraction length f refers to the number of fractional bits which is then used as the scale factor. 15 is the default.
• slope refers to the scalar you multiply the stored integer by to get the actual value (your real-world number). It is represented by the equation v = slope ∗ integer + bias. It differs from f because you can use non-binary scalars.
• bias lets you shift the scaling and slopeadjustmentfactor lets you use a potentially non-binary slope which must be . fixedexponent is also a scalar that lets you further adjust slope where slope = slopeadjustmentfactor ∗ 2^fixedexponent. These three parameters are not used much, so we will not discuss them further.

This would create a fixed-point number that is unsigned, 4 bits total, and of fraction length 6. Thus, the fixed-point number a is found by . Then, the next stored integer is . But wait, how is it possible to store a fraction length that is larger than the total number of bits? To store thie value, MATLAB automatically expands the number of bits to prevent overflow; in this case the word length becomes 5 bits total.

For multiplication, the result is "scaled twice." You multiply the two integers and the scale factor is . Thus for two integers a and b, you would do . In binary, the scaling factor is often a power of two.

For example, let one stored integer be 125 and the scaling factor be .

Let another stored integer be 33 scaled by .

We can multiply these together to yield .

The scale factor is then .

With this, we can then do to get .

To check our work, let's convert the original stored integers back to real numbers by doing and . Finally, . The two results match.

For fixed-point division, let's just jump straight into an example. We are dividing the stored integer 3124 with a scale factor of by the stored integer 9753 with a scale factor of .

The numerator is given by .

The denominator is given by .

The stored quotient is thus .

The final scale factor is given by . Thus, our calculated real result is .

For = stored integer values and , we have the general formulas for fixed-point multiplication and division:

Operation	Formula	Accuracy (generally)
Multiplication		Exact
Division		Approximate

The number has 3 significant digits and has 9 significant digits. In the real world, it's generally not useful for large numbers to be represented with more significant digits than small ones. If these were experimental data, it would mean that the in is the digit we are not fully confident in. But if we were to change it to say , it wouldn't make too much of a difference. If we did the same for and it became , that would be orders of magnitude more noticeable. Thus, this way of representing numbers is not very effective for scientific computing.

Using the same number line, let's start over. Let's zoom in on the interval from 1 to 2. For the following interval, there are values because of binary representation. But the exact number isn't as important as the concept. Just imagine that there is a certain amount of values on this interval.

Let's say we were to multiply those values by 2. They would then end up on the interval from 2 to 4:

There are still values on the new interval, they are just spaced out twice as far apart. The same idea applies if we were to multiply our original by :

where e is the exponent range.

We can also say that each interval has evenly spaced values in and that the spaces are relative to the size of the interval.

Suppose we have a real number and is the closest value in floating point to x. It then follows that the error is bounded by:

.

Thus, we can remove . We are making the upper bound more flexible which is a safe approximation:

This special value machine epsilon can then be calculated as which is the relative precision of floating-point arithmetic in this case. In real computing, floating-point arithmetic has around 16 decimal digits of precision, so our generalization checks out.

The set 𝔽 of floating-point numbers can be expressed as where and exponent range . The mantissa is the part of a floating-point number that represents its significant digits. It then follows that

where d is a fixed integer called the binary precision.

The binary representation of uses d bits thus . We will elaborate on this definition later.

where is the unsigned, fixed-point number encoded using 8 digits in binary with l bits for the integer part of the binary expansion and f bits for the fractional part of the binary expansion such that

with nonnegative integers

Input is an 8-bit binary word. Output is a rational number with where and . The output can be expressed as a finite binary expansion.

For instance, let . We will set and which lets us see that

This is then equivalent to .

If we were to set and , we would get a different value:

There are different combinations of values for and f such that .

.

where

Binary Word	Binary Value	Decimal Value
0000
0001
0010
0011
0100
0101
0110
0111
1000
1001
1010
1011
1100
1101
1110
1111

The range for x is . In general form that is .

.

where

Binary Word	Binary Value	Decimal Value
0000
0001
0010
0011
0100
0101
0110
0111
1000
1001
1010
1011
1100
1101
1110
1111

The range for x is . In general form that is .

.

where

Binary Word	Binary Value	Decimal Value
0000
0001
0010
0011
0100
0101
0110
0111
1000
1001
1010
1011
1100
1101
1110
1111

The range for x is . In general form that is .

.

where

Binary Word	Binary Value	Decimal Value
0000
0001
0010
0011
0100
0101
0110
0111
1000
1001
1010
1011
1100
1101
1110
1111

The range for x is . In general form that is .

.

where

Binary Word	Binary Value	Decimal Value
0000
0001
0010
0011
0100
0101
0110
0111
1000
1001
1010
1011
1100
1101
1110
1111

The range for x is . In general form that is .

From all of these cases, we can conclude that for , we can write the interpreted value of the finite binary expansion as

where the total number of bits

Assuming

where

We now get an approximation in the range of that can be encoded in our data class using a n-bit binary string.

Let . Expanding out gets us

Let's recast it with and . We have a new object represented by . You may recognize the hat notation from unit vectors in physics, when you have to approximate inconsistent systems in linear algebra, or to denote the Fourier transform of a function. In this case, we are using it to refer to the programming object used to store our approximation of x.

Notice how this has created an extra bit. We could round down to or round up to . Let's say we round down to get our approximation:

With , it follows that .

For a different method of dealing with the extra bit, we can round up. We will refer to this instance as .

Let's look at another instance of storing x in :

Previously, the least significant bits (LSB) were extra. In this case, the most significant bit (MSB) is extra. We could round down to or round up to . Let's say we round down to get our approximation and then find absolute error:

This makes sense because we approximated our exact 4-bit integer value (1010) as a 3-bit integer (110) and .

Data Type	Binary Value	Absolute Error	Type of Error
		0	None
			LSB
			LSB
		8	MSB
		4	MSB

Let's move on to another example. Let . Then,

This can be further simplified as . If we were to store this using 8 bits, what would happen? Well, for one, it would have to be an approximation with a good amount of error. This is inevitable, but we can still try to find the encoding with low error. Applying a strategy from the previous example, let's try a data type where and . We will also round up:

Because , it then follows that

We now have an absolute error in our approximation of . To convert this back into binary, notice that . Then we have as our absolute error in binary.

In our previous work, when we were given a number of bits , one possible way to represent it was to say that such that half of the bits in the number represented the integer part and half represented the fractional part. This would mean , , , or . We can analyze the ranges of these representations. Let's look at .

The lower bound is .

The upper bound is .

Thus, . More generally, . If you do this for the other possibilities, you will notice that the upper bounds, or the max integer values, will increase very rapidly.

After all of these examples, we can better see some disadvantages of fixed-point. One of these drawbacks is wasted bits. Using to encode would result in the following string of bits . In high performance environments, this is a waste of storage that will also hurt performance.

Another drawback is that the range is very specific. For instance, has an upper bound of . The size of this would be okay for applications, but it's not nearly robust enough. If we needed to represent larger values, we could increase l, but this would result in us having less precision with the fractional part of a number. Remember that having more significant digits with lower numbers is generally more important than having the same with large numbers.

Let's say we wanted to represent . In fixed-point binary, this would become . All of these zeros do not really add anything important to the representation in terms of precision. They are only there to signify the location of the radix point. And in fixed-point decimal, this would become . We can represent the number much more efficiently by allowing the radix point to float:

We can then float the radix point four places to the right to get .

We have the number which is equal to in binary. To represent this in fixed-point, we would need 8 bits for the fractional part. In floating-point, we can represent it as follows:

,

,

where . We then have

,

where

Aspect	Fixed-Point	Floating-Point
General form
Stores	Integer part and fractional part	Sign bit, mantissa, exponent
Analogous to	Decimal numbers scaled by a fixed decimal point	Scientific notation
Range	Limited by fixed scaling	Controlled by dynamic exponent scaling
Resources	Very little, inexpensive	Generally fast, depends on hardware but is more expensive
Created by	Early computer architects and engineers	Numerical analysts
Used for	Digital signal processing, controls systems	Modern scientific computing

1.3 IEEE Standard for Floating-Point Arithmetic (IEEE 754)

Let's work with quarter-precision floats. Let the total number of bits be 8 and let the number be normalized. Then, we have the bitwise format . The 8-bit binary word looks as follows:

Let's try it out. Suppose . It then follows that

Since , . means that . And means that . We are thus going to store :

where + is the sign bit and 1 is the implied normalized leading bit. We can now say that , , and .

Let's try an example with a larger number :

We can now say that , , and . The value is stored as:

Recall that IEEE 754 doesn't store the actual exponent e but rather encodes something called the biased exponent u using the bias k. In this case, . Thus . Plugging our values in, we can find the actual value being stored .

We can now come up with a table that summarizes how this works. Recall that u is the value that is stored by the computer for the 3-bit exponent field . The exponent value is the 3-bit interpreted value of e in decimal found by . This is known as excess-3.

Exponent Bits	Biased Exponent u	Exponent Value e
000	0
001	1
010	2
011	3	0
100	4	1
101	5	2
110	6	3
111	7	4

Quarter-precision 8-bit Word	Sign Bit	Exponent Code	Numerical Value (floating-point binary)	Numerical Value as Sum (decimal)	Numerical Value (decimal)
11100000	1	110
11011111	1	101
11011110	1	101
11011101	1	101
11011100	1	101
11011011	1	101
11011010	1	101
11011001	1	101
11011000	1	101
...	...	...	...	...	...
01011000	0	101
01011001	0	101
01011010	0	101
01011011	0	101
01011100	0	101
01011101	0	101
01011110	0	101
01011111	0	101
01100000	0	110
...	...	...	...	...	...
01101111	0	110