Skip to the content of the web site.
## Laboratory 9

# 9.1 4th-order Runge-Kutta and the Dormand-Prince method

## 9.1.1 Additional Matlab Functionality

## 9.1.2 4th-order Runge Kutta

# 9.2 Linear Regression and Least Squares

## 9.2.1 Derivation

## 9.2.2 Matlab and Least Squares

## 9.2.3 Best-Fitting Least-Squares Curves

This laboratory will investigate two topics:

- MATH 211: 4th-order Runge-Kutta and the Dormand-Prince method, and
- MATH 215: Least Squares.

In this topic, we will discuss higher-order methods for solving initial-value problems; however, to start, we will begin with
additional features in Matlab. Rather than teaching 1001 different Matlab commands in the first class, it is more useful
to introduce commands when they their use is apparent. It's like learning bridge and the various
contract bridge bidding conventions: until you have played
for a while, it's very difficult to understand the purpose of, for example,
a *negative double*.

We will begin with defining functions in Matlab directly (as opposed to creating explicit `.m` files. For example, consider
defining the function :

y = @(t)(3.2*t.*cos(1.52*t).*exp(-0.74*t));

You can now treat this as any other user-defined function:

y(1) ans = y = @(t)(3.2*t.*cos(1.52*t).*exp(-0.74*t)); t = linspace( 0, 10, 1000 ); plot( t, y(t), 'r' )

It is also possible to define multi-variate functions but, in this case, the variables should be specified:

f9a = @(t, y)(-t*y + t - 1) f9a = Inline function: f9a(t,y) = -t*y + t - 1 f9a(0,1) % the slope at (0,1) ans = -1

Other useful commands:

`clc`- Clear the command window: begin working with a blank screen.
`clear M`- Unassign the variable
`M`. Often useful if you accidentally assign a useful function`plot = [1 2 3]'`.plot = [1 2 3]'; plot( x, y ) % fails ??? Subscript indices must either be real positive integers or logicals. clear plot plot( x, y ) % okay

Without any arguments,`clear`unassigns all variables. `diary`- Records all the commands in Matlab and their corresponding output.

At all times, you can always use `help cmd` to get more information about the Matlab function

Recall that Heun's method is as follows:

- Find the slope at by calculating ,
- Find the slope at by calculating , and
- Use the average of the slopes to find the next approximation .

We have now sampled the slope at two points: and . One may reasonably expect that further samples may be used to reduce the error even further. The 4th-order Runge-Kutta algorithm proceeds as follows:

- ,
- ,
- , and
- .

The approximation of is found by taking a weighted average of the four slopes:

.

You will notice that the coefficients in the numerator add up to .

We could compare Euler's, Huen's and 4th-order Runge-Kutta, for example:

Method | Value |
---|---|

Euler | 0 |

Huen | -0.5 |

4th-order Runge-Kutta | -0.2557977040 |

Exact | -0.2484078630 |

however, this would be unreasonable: Runge-Kutta uses four evaluations of the slope while Euler uses only one and Heun uses two. Thus, it would be more reasonable to allow for four steps with with Euler's method and two steps with with Huen's method:

Method | Steps | Value |
---|---|---|

Euler | 0.75, 0.51953125, 0.2407703400, -0.1389852180 | -0.1389852180 |

Huen | 0.4687500000, -0.3137178270 | -0.3137178270 |

Runge-Kutta | -0.2557977040 | -0.2557977040 |

Exact | -0.2484078630 |

Similarly, we can approximate , and this yields

Method | Steps | Value |
---|---|---|

Euler | 0.8750000000, 0.7600097656, 0.6452477295, 0.5221857591 | 0.5221857591 |

Huen | 0.7597656250, 0.5063127134 | 0.5063127134 |

Runge-Kutta | 0.5187200043 | 0.5187200043 |

Exact | 0.5187154179 |

The error of 4th-order Runge Kutta is actually ; however, as with Euler and Heun, if you use it in a multi-step method, the error is reduced by one order to .

To test your code for `rk4` which implements the 4th-order
Runge Kutta algorithm, you could try the following.
The calling sequence should be the same as for `heun.m` from the previous
laboratory.

In order to ensure that your function is working, you must run the following tests with the differential equation :

2*exp(-0.5) ans = 1.213061319425267 g = @(t, y)(-y); [ts, ys] = rk4( g, [0, 0.5], 2, 2 ) % y(0) = 2 and approximate y(0.5) with two points (one step) ts = 0 0.500000000000000 ys = 2.000000000000000 1.213541666666667 [ts, ys] = rk4( g, [0, 0.5], 2, 3 ) % y(0) = 2 and approximate y(0.5) with three points (two steps) ts = 0 0.250000000000000 0.500000000000000 ys = 2.000000000000000 1.557617187500000 1.213085651397705 [ts, ys] = rk4( g, [0, 0.5], 2, 5 ) % y(0) = 2 and approximate y(0.5) with five points (four steps) ts = 0 0.125000000000000 0.250000000000000 0.375000000000000 0.500000000000000 ys = 2.000000000000000 1.764994303385417 1.557602445491486 1.374579721615834 1.213062689100529 [ts, ys] = rk4( g, [0, 0.5], 2, 1025 ); % y(0) = 2 and approximate y(0.5) with 1025 points (1024 steps) ys(end) ans = 1.213061319425268

Consider the following attempt to approximate for the same initial-value problem described above:

[ts, y1] = euler( g, [0, 0.1], 2, 2 ) ts = 0 0.100000000000000 y1 = 2.000000000000000 1.800000000000000 [ts, y2] = heun( g, [0, 0.1], 2, 2 ) ts = 0 0.100000000000000 z1 = 2.000000000000000 1.810000000000000

The correct solution is . Consider the value

dy1 = abs( y1(end) - 2*exp(-0.1) ) dy2 = abs( y1(end) - y2(end) )

The first is the absolute error of the approximation of using Euler's method.
Recall that Heun's method is significantly better than Euler's method and therefore the second value, `dy2`
will represent an approximation of `dy1`. What is the relative error of
`dy2` as an approximation to `dy1`?

Would it be reasonable to use a higher-order approximation to approximate the error of a lower-order approximation in the manner which has been described here?

In this topic, we will investigate the problem of finding the best-fitting line or curve to data with noise.

Suppose that we have a collection of data points which were read form an input source with noise, for example

.

A plot of this data is in Figure 1.

Figure 1. Twenty noisy data points.

The following is the Matlab code for the data.

xs = [0.3279750853 2.464595439 2.717370800 2.736334037 2.953520867... 3.427599402 3.982231181 3.985682362 4.269594141 4.413798254... 4.544029287 5.781998567 7.164764166 7.331776715 7.467984805... 7.777045621 7.916769076 9.051432635 9.492962523 9.904912835]'; ys = [2.407577618 3.203292359 3.390732556 3.559771990 3.749224638... 4.192705730 4.556825783 3.923328666 4.511085438 4.962406558... 4.925002895 5.507543465 6.464099937 5.957215200 6.616048328... 6.907485826 6.750967379 7.633998055 8.116679158 8.334500002]';

Clearly the data appears to be approximately linear, and if we wanted to approximate this data, we could find the interpolating polynomial. Unfortunately, this turns out to be the interpolating polynomial shown in Figure 2.

Figure 2. The polynomial of degree 19 interpolating the 20 points shown in Figure 1.

As you may guess, this is completely inappropriate—just the constant coefficient of the interpolating polynomial is on the order of . Consequently, we would like to find an approximating line. The question is, which is the best possible line? For example, Figure 3 shows two possible lines which appear to be a reasonable approximation of the data.

Figure 3. Two possible lines approximating the data in Figure 1.

The question is, what criteria will we use to find the best possible line? First, the line we are finding is of the form (we will use the Greek letters alpha and beta as these are standard in the literature). Given any line we have an error associated with the line and how well it approximates our points. Figure 4 shows a zoom of the points through and let us suppose our approximating line is shown in blue. In this case, the point is the approximation of the point . Similarly, the point on the line approximates the point .

Figure 4. A zoom on some of the points on the right-hand plot of Figure 4.

The black lines in Figure 4 represent the errors. We could, for example, attempt to choose a line which minimizes the sum of all the absolute errors, that is, we could minimize

.

Unfortunately, for theoretical and practical reasons, this will not work. Instead, we choose to minimize the sum of the squares of the errors:

.

Here we have a function of only two variables: and —remember that all the and are already given: they are constants.

From first-year calculus, you have learned to to find a minimum: differentiate, equate to zero, and solve for the variable. In this case, we have two variables and therefore we must have two partial derivatives, equate each partial derivative to zero, and simultaneously solve for the two variables. Using the chain rule, we get:

and

.

Expanding these two equations and dividing by 2, we get

.

Because and are not related to the index, we may expand this to get

.

Now we come back to Math 215 *Linear Algebra*: this is the system of linear equations

.

Without proof, the determinant of the matrix is zero if and only if .

In order to quickly calculate the matrix in Matlab, notice that if we define

M = [xs xs.^0]; M' * M ans = 719.2291 107.7124 107.7124 20.0000 M'*ys ans = 657.5218 105.6705 p = (M' * M) \ (M' * ys) p = 0.6355 1.8609

The last is our solution to our minimization problem; that is, the best-fitting line (according to our definition of minimizing the sum of the squares of the errors) is

.

This line is called the *best-fitting least-squares line* and the process
to find it is called *linear regression*. A plot of the line together with
the points is shown in Figure 5.

Figure 5. The best-fitting least-squares line passing through our data.

The line can be plotted creating, for example, 100 points using `xs = linspace( 0, 10, 100 )`
and plotting those points against `polyval( p, xs )`.

It is also possible to calculate the sum of the squares of the error and demonstrate that even a slight change in the coefficients will produce a line which has a larger sum of squares of the errors:

yp = polyval( p, xs ); norm(yp - ys)^2 ans = 1.033686445547913 % A slightly different polynomial ypp = polyval( [0.6354 1.8608]', xs ); norm(yp - ys)^2 ans = 1.033698337287500

You will notice that the run time of calculating is and thus, it is only an calculation; however, Matlab even provides a short cut:

With the above matrix `M`, suppose you asked

M \ ys

What would you expect? From linear algebra, we have a system of 20 equations but only two unknowns. Except in the most degenerate cases, no answer exists to such a system. Rather than returning an error, Matlab assumes you are looking for the best-fitting least-squares line:

(M' * M) \ (M' * ys) ans = 0.635520678531010 1.860852407821792 M \ ys ans = 0.635520678531010 1.860852407821796

There is a slight difference in the answers, but it is only in the last decimal place.

What happens, however, if we don't want to find the best straight line? Consider the data in Figure 6.

Figure 6. Another set of ten data points with noise.

The data points are

and the exact data points are

xs = [0.007814731122 0.008848343076 0.04502115842 0.07825884434 0.3100781521... 0.3216720831 0.3359594404 0.6089038628 0.8011831826 0.8078659323]'; ys = [1.543701192 1.557800839 1.420389441 1.402311968 1.338448518... 1.296099571 1.306500894 1.805512249 2.561774915 2.512208861]';

An important observation is that **you can fit a straight line to any data**, but
what is important is: does that line have any significance what-so-ever?

For example, the best-fitting line is given by

M = [xs xs.^0]; p = M \ ys p = 1.2307 1.2652

and this is plotted in Figure 7; however, is this really the best curve we fit to the data? Probably not.

Figure 7. The best-fitting least-squares line passing through the data in Figure 6.

Note that the data initially dips and then appears to rise again. We could attempt to fit a quadratic function to this data:

M = [xs.^2 xs xs.^0]; p = M \ ys p = 4.0835 -2.0560 1.5472

If we plot this quadratic function , as is shown in Figure 8, we see that it appears to be a good fit.

Figure 8. The best-fitting least-squares quadratic function passing through the data in Figure 6.

Is this the best-fitting curve? Should we try a cubic function, or a quartic? With each extra term, the fit will be better and we know that if we try a polynomial of degree nine, the error will be zero—it will be a perfect fit.

This is a question about having a model for the data: does the model suggest a quadratic function? Are we taking the data from a system which suggests that the sampled data would behave quadratically as the variable is changed?

As another example, assume we are testing a resistor by placing various currents through the resistor and measuring the voltages. For example,

These points are given in terms of and :

current = [0.03465445112 0.02764957747 0.08053249405 0.09868266998 0.08720882767 ... 0.06688377664 0.08705955749 0.04056978469 0.04319976668 0.05855719319]'; voltage = [4.379276551 3.301032219 10.01909723 11.62485446 10.32856468 ... 7.831413695 10.55234404 4.939755206 5.418582132 7.110929634]';

Figure 9. Current-voltage readings taken across a resistor.

You should note that the data, in this example, is not ordered, nor need it be.

If we find the best-fitting least-squares line, shown in Figure 10, we have

M = [current current.^0]; p = M \ voltage p = 117.3677 0.2151

Figure 10. The best-fitting least-squares line .

This says that the voltage across the resistor is defined by

which should instantly raise the alarm in the mind of any electrical or computer engineering student (if it didn't, go back and review ECE 140). If the current is zero, the voltage must be zero.

Thus, rather than using the model that , we should use the model , in other words, we want to find the best line which passes through the data but also passes through the origin.

M = [current]; r = M \ voltage r = 120.3697

Thus, the estimate on the resistance of the resistor is and this best-fitting line passing through the origin is shown in Figure 11.

Figure 11. The best-fitting least-squares line .

To practice linear regression, you could try the following: find the best-fitting line (using linear regression) which passes through the points

x1 = [1.3862 1.4929 2.5751 5.4722 6.9908 8.4072 8.9090 9.5929]; y1 = [4.6006 4.8822 5.4627 6.9211 7.9681 8.7056 9.3139 9.4706];

Plot both the eight points and a plot of the best-fitting straight line you found using the
`polyval` function. Use an appropriate interval slightly larger than the x-values of the points which were found.

Find the best-fitting line (using linear regression) which passes through the points

x2 = [1.3862 1.4929 2.5751 5.4722 6.9908 8.4072 8.9090 9.5929]; y2 = [0.4191 1.1449 4.5667 0.7619 4.1291 2.6917 4.9807 0.3909];

Plot both the eight points and a plot of the best-fitting straight line you found using the
`polyval` function. Use an appropriate interval slightly larger than the x-values of the points which were found.

Two questions:

- Is this the best-fitting line passing through these points? (yes or no)
- Just because you can fit a line to data, does that necessarily mean the line has any significance? (yes or no)

Find the best-fitting line which passes through the points

x3 = [1.3862 1.4929 2.5751 5.4722 6.9908 8.4072 8.9090 9.5929]; y3 = [3.5101 6.6099 2.1850 -16.1321 -33.2007 -47.7779 -55.9318 -67.5378];

Plot both the eight points and a plot of the best-fitting straight line you found using the
`polyval` function. Use an appropriate interval slightly larger than the x-values of the points which were found.

Now, find the best-fitting quadratic polynomial which passes through the eight points.
Plot both the eight points and a plot of the best-fitting quadratic polynomial you found using the
`polyval` function. Use an appropriate interval slightly larger than the x-values of the points which were found.