7.2 Integration Using Straight Lines
227
We shall now take a while to derive the method we are going to use to integrate
the function on the grid of points. We shall use N points and as such we use
the code
#
step = (b-a)/(N-1);
x = a:step:b;
f = integran

226
7. Numerical Integration
In this case the area under the curve is approximately
7 4 the area of the boxes
(the fourth row up has a few part squares as does the fth row). We obviously
need a scheme which is slightly more robust.
7.2 Integration Using S

6.9 Tasks
223
Task 6.21 Determine the eigenvalues
1 0
0 1
0 0
1 0
of the matrix
0 1
0 0
.
1 0
0
1
Task 6.22 Prove that An = PDn P1 using induction where n N and P
and D are comprised of the eigenvectors (as columns) and eigenvalues of A
respectively.

7
Numerical Integration
7.1 Introduction
In this chapter we shall discuss techniques whereby functions can be integrated;
these are quite classical. We will give the derivation of the approximations and
will mention the likely failings of the techniques.

222
6. Matrices
Task 6.17 Determine whether the following systems have solutions (and if
these are unique):
3x + 2y = 7
3x 2y = 7;
x+y+z
x+yz
x+y
=1
=0
= 0;
x+y+z+a+b+c=1
xy+z+a+b+c=1
x+yz+a+b+c=1
x+y+za+b+c=1
x+y+z+ab+c=1
x + y + z + a + b c = 1.
Task 6.

6.9 Tasks
221
Task 6.12 Calculate the quantities
3
1 1 2
2 and
3
0 1
1
5 2
1
2
4
2
0
1
1 1
2 1
.
Task 6.13 Calculate, by hand, the quantities
a
c
b
d
1
0
0
1
and
1
0
0
1
a b
c d
.
As we saw earlier a matrix with ones on the leading diagonal (running from

220
6. Matrices
We have also seen that it is possible to use dot arithmetic, which acts
element by element on matrices and vectors.
Task 6.7 Show that for all square matrices A, the matrices B = A + AT and
C = A AT are symmetric and anti-symmetric respect

6.9 Tasks
219
9
A = ones(3);
B = 2*ones(3,2);
C = 3*ones(2,3);
6
[A B]
[A B]
[A C B]
[A C]
[A; C]
[A; B]
8
7
At this stage you are now in a position to perform relatively simple operations on matrices; we shall start with a couple of examples of addition

218
6. Matrices
0
0
0
3
0
0
3
4
0
3
4
5
0
4
5
You should recall the nth diagonal is shorter by |n| than the leading diagonal.
In order to understand this code it may help you to type r(1:4) (which gives
the rst four elements of r) and r(2:5) (which gives

6.9 Tasks
217
of the eigenvalues are real and positive the solution will expand as t increases.
Similarly if both are real and negative the solution will contract. This analysis
is very powerful and allows us to understand the structure of these problems

216
6. Matrices
the previous expression for An we nd that An+1 can also be expressed in terms
of a similar polynomial. In fact An can be written as a polynomial of degree
n 1 in A. Since the exponential form of the matrix involves only powers of
the matri

6.8 Exponentials of Matrices
215
We now have to dene and construct exp(A). We do this using the summation
1 i
exp(A) =
A.
i!
i=0
Again we can use the diagonalised form of the matrix to yield
exp(A) = P
i=0
1 i
D
i!
P1 ,
but this gives
exp(A) = P exp(D)P1

214
6. Matrices
We have elected to use the points xi = (i 1)/n since these will not coincide
with eigenvalues unless the elements of the matrix involve . We can now
consider the roots of this polynomial to determine the eigenvalues of the matrix
(using ro

6.7 Characteristic Polynomials
213
have coecients a0 , , an . Consequently we have
n
ai Ai = 0,
i=0
where we use the natural convention that A0 = I. We can extract one term
from the summation so that
n
ai Ai = 0.
a0 I +
i=1
Notice that p(0) = a0 so that a

212
6. Matrices
8.4127
-4.8097
stopcrit =
3.2411e-15
=
V =
0.0769
0.6050
0.0143
0.7674
0.1974
0.1400
0.5780
0.0270
-0.8034
0.0100
8.4127
0
0
-4.8097
D =
6.7 Characteristic Polynomials
We shall pause here to discuss the characteristic polynomial of a matri

210
6. Matrices
Example 6.17 Calculate the matrix
1
A= 0
1
0
1
0
1
1
0
raised to the hundredth power.
This is accomplished using the code
5
a = [1 0 -1; 0 1 1; 1 0 0];
[v,d] = eig(a);
d = sparse(d);
d100 = d.100;
a100 = v*d100*inv(v);
4
2
3
We also menti

208
6. Matrices
which are solved to give = 2/5 and = 1/5. Hence
1
1
=
1
2
e=4 + e=1 .
5
5
Now we multiply both sides by A to give
A
1
1
=
=
=
=
=
2
1
Ae=4 + Ae=1 ,
5
5
2
1
(4)e=4 + (1)e=1 ,
5
5
8
1
e=4 e=1 ,
5
5
8
1
2
1
,
3
1
5
5
3
5
.
We could have used

6.6 Specic MATLAB Commands
209
where P is a matrix formed with the eigenvectors of the matrix A as columns
(that is V) and D is a diagonal matrix with the eigenvalues on its diagonal.
Note: the order of these must be the same.
This is veried for the above

6.5 Eigenvalues and Eigenvectors
207
which gives
x + 2y
= x
y = x,
3x + 2y
= y
x = y.
Hence the eigenvector is
e=1 =
1
1
.
We could have determined these using the command
9
[V,D] = eig(A)
6
V =
-0.7071
0.7071
-0.5547
-0.8321
D =
-1
0
0
4
8
7
This retur

6.5 Eigenvalues and Eigenvectors
205
In this simple case it is easy to see how these quantities are determined but
consider the more involved example:
Example 6.16 Determine the eigenvalues and eigenvectors of the matrix
1
3
2
2
,
and comment on the eect

206
6. Matrices
ans =
4
-1
>
The rst command gives the coecients of the characteristic polynomial (which
is the name of the polynomial in , see page 214 for more details). The second
one determines the roots of polynomial whose coecients are stored in p.

204
6. Matrices
full one if the result is full). A full matrix can be retrieved by using the full
command.
6.5 Eigenvalues and Eigenvectors
Although this is slightly beyond the scope of this text we shall briey discuss
the ideas of eigenvalues and eigenve

6.4 Matrix Decomposition
203
matrix should be exploited.
MATLAB has a few commands which can use the properties of these matrices: perhaps the most simple one is sparse. We will meet examples in due
course but we shall simply show this example.
> A = spar

202
6. Matrices
9
%
% mtestmat.m
%
nn = 100:50:800;
n = length(nn);
6
for i = 1:n
nn(i)
[t1(i),t2(i)] = testmat(nn(i);
end
plot(nn,[t1; t2])
xlabel(Matrix size,FontSize,15)
ylabel(CPU time (secs),FontSize,15)
text(600,10,Decomposition,FontSize,14)
text(50

6.4 Matrix Decomposition
201
where we have introduced E = L1 P. This actually occurs earlier in the series
of operations since EA = U. The algorithm can be written as:
(i) Calculate L, U
(ii) Calculate Ly = b (y = Eb)
(iii) Solve Ux = y.
This method is ec

200
6. Matrices
1.0000
0.3333
0.3333
0
1.0000
0.2500
0
0
1.0000
3.0000
0
0
2.0000
1.3333
0
1.0000
0.6667
0.5000
U =
P =
0
1
0
1
0
0
0
0
1
Now we have the three matrices L, U and P, such that
PA = LU.
Notice it is possible to call this routine with only tw

6.4 Matrix Decomposition
199
and
2x + 3y
xy
= 2,
=
8.
These equations can be written as
2
1
3
1
x
y
=
7
1
2
8
This system can be solved using the MATLAB code:
#
A = [ 2 3; 1 -1];
b = [7 -2; 1 8];
x = A\b;
"
This gives
.
!
x =
2.0000
1.0000
4.4000
-3.6000

198
6. Matrices
The above technique involves other methods than those we used in the
simple two-by-two case (for instance pivoting). It is not the intention of this
text to provide a comprehensive discussion of these matrix operations, merely
to make you