程序代写代做代考 matlab algorithm my MAT420 Lecturenotes 5 Programming

my MAT420 Lecturenotes 5 Programming
In this lecture we consider different implementations of various fondamental problems using Matlab, and compare their ease of implementation and efficiency.
5.1 Powers
Compute an for a real (possibly a matrix) and n integer
􏰐
an = a · an−1 recursive algorithm
The strategy corresponds to a (1, n − 1)-split of n (one-sided). The corresponding number of operations is
Wn =1+Wn−1 =…=n−1+W1 =n−1. For example, for n = 9:
a9 =a·􏰉a8􏰈
a8 =a·􏰉a7􏰈
a7 =a·􏰉a6􏰈
a6 =a·􏰉a5􏰈 a5 =a·􏰉a4􏰈
a4 =a·􏰉a3􏰈 a3 =a·􏰉a2􏰈 a2 = a · a
requires 8 multiplications.
power recursive.m recursive power function
1 function a = power_recursive(a,n)
2 if n==1, return; end % nothing to do: output a same as input
3 a = a*power_recursive(a,n-1); % for n>1 use the recursive definition
5.1

page 5.2 􏰝5.1. POWERS
􏰐
divide-and-conquer algorithm
This strategy splits n (approximately) evenly (more balanced). The number Wn of multiplications depends on n.
1. If n is even the cost Wn of computing an is the cost of computing an/2 plus one multiplication.
2. In particular, for n = 2p,
Wn =1+Wn/2 =1+1+Wn/4 =…=1+…+1+W1 =p=log2n≪n−1
􏰆 􏰅􏰄 􏰇
p times
3. If n is odd, the number of multiplications needed can be shown to be less than 2 log2 n if intermediate computations are saved and reused (dynamic program- ming). For example, for n = 9:
a9 =􏰉a4􏰈·􏰉a5􏰈 a5 =􏰉a2􏰈·􏰉a3􏰈 a4 =􏰉a2􏰈·􏰉a2􏰈 a3 =a·􏰉a2􏰈
a2 = a · a
requires 5 multiplications.
power divideconquer.m recursive power function with balanced partitioning
􏰐
􏰑BDW
Draft August 22, 2013
1 2 3 4 5 6 7
function a = power_divideconquer(a,n)
if n==1, return; end % nothing to do: output a same as input
if mod(n,2) % remainder of n/2; = 1 (True) if n odd, = 0 (False) if n even
a = power_divideconquer(a,(n-1)/2)*power_divideconquer(a,(n+1)/2); % n odd else
a = power_divideconquer(a,n/2)ˆ2; % n even end
mixed algorithm
an = 􏰘an/2 · an/2 if n even, a(n−1)/2 · a(n+1)/2 if n odd
for large n. For example, for n = 8: a8 =􏰉a4􏰈·􏰉a4􏰈
a4 =􏰉a2􏰈·􏰉a2􏰈 a2 = a · a
requires 3 multiplications.
an = 􏰘an/2 · an/2 if n even, a·an−1 if n odd
This strategy combines the recursive and divide-and-conquer approaches and can be more efficient than either one. For example, for n = 9:
a9 =a·􏰉a8􏰈
a8 =􏰉a4􏰈·􏰉a4􏰈 a4 =􏰉a2􏰈·􏰉a2􏰈
a2 = a · a

􏰝5.1.
POWERS page 5.3
requires 4 multiplications.
power mixed.m recursive power function with mixed partitioning
function a = power_mixed(a,n)
if n==1, return; end % nothing to do: output a same as input
if mod(n,2) % remainder of n/2; = 1 (True) if n odd, = 0 (False) if n even
a = a*power_mixed(a,n-1); % n odd else
a = power_mixed(a,n/2)ˆ2; % n even end
2
symmetric) and powers of a do not grow in size. Note that the default maximum number of recursive levels is 500 in Matlab, which needs to be increased here for the first algorithm (this is risky, since it increases the memory stack needed for the computation, which can lead to Matlab (or the computer) crashing)
>> format long % more digits in output compared to default format short >> set(0,’RecursionLimit’,1021) % reset recursion limit
>> A = [1 1;1 0];
>> r = max(abs(eig(A))); % spectral radius of A
>> a = A/r; % spectral radius of a is 1, equal to 2-norm of A (symmetric) >> n = 1019;
>> tic, power_recursive(a,n), toc
ans =
0.723606797749976 0.447213595499956
0.447213595499956 0.276393202250020 Elapsed time is 0.029440 seconds.
>> tic, power_divideconquer(a,n), toc ans =
0.723606797749909 0.447213595499915
0.447213595499915 0.276393202249994 Elapsed time is 0.005978 seconds.
>> tic, power_mixed(a,n), toc
ans =
0.723606797749939 0.447213595499933
0.447213595499933 0.276393202250006
Elapsed time is 0.002562 seconds.
>> tic, aˆn, toc % compare with naive power using Matlab ˆ ans =
0.723606797749908 0.447213595499914
0.447213595499914 0.276393202249994 Elapsed time is 0.000034 seconds.
1 2 3 4 5 6 7
The quantity a can be a scalar or a (square) matrix. We compare the 3 strategies for a=1A, A=􏰖1 1􏰗, and n=1019,
(5.1) with r = 1+√5 (Golden ratio and spectral radius of A), so that ρ(a) = ∥a∥2 = 1 (a is
r10
􏰑BDW Draft August 22, 2013

page 5.4 􏰝5.1. POWERS
As expected, the mixed approach is faster than the divide & conquer approach, which is itself faster than the standard recursive strategy. The reason the command an is (much) faster is twofold: 1) it is executed as a (pre-compiled) C procedure by Matlab, so that no checking of syntax is performed during execution, and 2) it is performed using a loop inside a single call, rather that without a loop but using recursive programming. An implementation of the basic strategy using such a loop would read something like:
power direct.m power function with for loop
1 function an = power_direct(a,n)
2 an = a;
3 for k = 1:n-1
4 an = a*an;
5 end
Note that, while the recursive approach is “top-down”, i.e., starts with an in terms of an−1 in terms of an−2 . . . the direct approach is “bottom-up”, i.e, the powers of a are built from the ground up, starting with a1 = a. This is a problem when implementing a strategy like the mixed approach, because it requires advanced knowledge of the sequence of powers of a needed to be computed. This can be achieved by pre-determining odd/even steps from the top down without actually computing the powers:
power mixed direct.m power function, mixed recursive strategy, with for loop
1 2 3 4 5 6 7 8 9
10
11
12
13
14
15
16
17
function an = power_mixed_direct(a,n) steps = [];
while n>1
if mod(n,2)
n = n-1; steps = [steps, 1]; % n odd
else
n = 9 example steps: []
n = n/2; steps = [steps, 0]; % n even
[1]
[1,0]
[1,0,0]
[1,0,0,0]
9
8
4
2
1
a9
a8
a4
a2
a :an
end
an=a; for k =
end
length(steps):-1:1 if steps(k)
an = a*an; % odd step else
an = anˆ2; % even step end
end
>> tic, power_direct(A,n), toc ans =
0.723606797749976 0.447213595499956
0.447213595499956 0.276393202250020 Elapsed time is 0.001347 seconds.
>> tic, power_mixed_direct(A,n), toc
ans =
0.723606797749939 0.447213595499933
􏰑BDW
Draft August 22, 2013

􏰝5.2. RECURRENCES – FIBONACCI page 5.5 0.447213595499933 0.276393202250006
Elapsed time is 0.000340 seconds.
The corresponding timings are lower compared to the recursive strategies. On the other hand, the programing can become substantially more difficult and harder to read/maintain.
P5.1.1 Show that the matrix a from (5.1) satisfies
􏰙1+1 1􏰚
liman= 2 12√5 1 √51 n→∞ √5 2 − 2√5
(hint: determine the spectral/eigenvalue decomposition a = XDX−1 of a). Which of the above algorithms provides the best accuracy for a1019?
P5.1.2 The recursive program power recursive requires a memory stack of 1018 recursion levels (the function with n = 1019 calls the function with n = 1018, which call the function with n = 1017, and so on until the function with n = 2 calls the function with n = 1, which does not call anything). How tall is the memory stack for the power divideconquer and power mixed programs when n = 1019?
P5.1.3 Modify power mixed direct.m to obtain a direct implementation of the divide & con- quer strategy. Compare the timing for the matrix a and the value n = 1019 used in this section.
5.2 Recurrences – Fibonacci
As an example of 3-term recurrence with constant coefficients, consider the Fibonacci num- bers defined by F0 = 0,
F1 = 1,
Fn=Fn−1+Fn−2 forn>1
These numbers have many interpretations. For example, Fn+1 represents the number of ways to set-up a wall 2 feet high with n (2D) bricks of size 2 × 1 feet:
F1 =1 F2 =1 F3 =F2+F1 =2 F4 =F3+F2 =3 Proof. The initial values match and the recurrence follows from
is or
Fn+1 = Fn + Fn−1
􏰞
(n − 1) × 2
(n − 2) × 2
n×2
The following function determines Fn recursively, but is not practical because of too many redundant computations, see P5.2.1 :
􏰑BDW Draft August 22, 2013

page 5.6 􏰝5.2. RECURRENCES – FIBONACCI fibonacci recursive.m Generate Fibonacci number using slow recursion
function F = fibonacci_recursive(n) if n==0
F = 0;
elseif n==1
F = 1;
else
F = fibonacci_recursive(n-1)+fibonacci_recursive(n-2); end
The non-recursive implementations
fibonacci direct.m Generate Fibonacci number with for loop
Fn
Fk
Fk−1 Fk−2
F0
function F = fibonacci_direct(n)
F = zeros(n+1,1);
F(2) = 1; % F(2) is second F value, i.e. F_1 = 1 for k = 3:n+1
F(k) = F(k-1)+F(k-2);
end
F = F(n+1);
and
F array
fibonacci direct2.m Generate Fibonacci number with for loop and low storage
function F = fibonacci_direct2(n) if n==0
F = 0;
else
Fn1 = 0; F = 1; for k = 2:n
Fn2 = Fn1; Fn1 = F;
F = Fn1+Fn2;
end end
F Fk−1 Fk−1
Fn1 Fk−2 Fk−1 Fk−1
Fn2 Fk−3 Fk−2 Fk−2
The later implementation requires only 3 storage values, as previously computed numbers are overwritten with new numbers.
Improved recursive implementation
Proof.
􏰖Fn+1 Fn 􏰗=An with A=􏰖1 1􏰗 Fn Fn−1 10
􏰖Fn+1 Fn 􏰗=􏰖1 1􏰗􏰖 Fn Fn−1􏰗=…=􏰖1 1􏰗n􏰖F1 F0 􏰗=An􏰖F1 Fn Fn−1 1 0 Fn−1 Fn−2 1 0 F0 F−1 F0
if we introduce F−1 = F1 −F0 = 1 so that the recurrence also holds for n = 1. 􏰑BDW
􏰞
F0 􏰗=An􏰖1 0􏰗=An F−1 0 1
Draft August 22, 2013
Fk

􏰝5.2. RECURRENCES – FIBONACCI
page 5.7
It is sufficient to consider n = 2m even. Then
􏰖F2m+1 F2m
F2m 􏰗 = A2m = (Am)2 = 􏰖Fm+1 F2m−1 Fm
Fm 􏰗2 Fm−1
Fm(Fm+1 + Fm−1)􏰗 Fm2 + Fm2 −1
(m − 1) ×2
(m − 1) ×2
m×2
m×2
2m × 2
is
= 􏰖 Fm2 +1 + Fm2 Fm(Fm+1 + Fm−1)
or
F2m+1 = Fm+1 × Fm+1 + Fm × Fm
(2m − 1) × 2
is
(m − 2) or ×2
(m − 1) ×2
F2m = Fm × Fm+1 + Fm−1 × Fm Note that F2m = Fm(Fm+1 + Fm−1) = (Fm+1 − Fm−1)(Fm+1 + Fm−1) = Fm2 +1 − Fm2 −1.
Thus 􏰘Fn2/2+1 − Fn2/2−1 if n even Fn= F2 +F2 ifnodd
(5.2)
(n+1)/2 (n−1)/2
Note that (5.2) cannot be used for n = 2 (set F2 = 1 as part of initialization). A recursive
implementation reads
fibonacci recursive2.m Generate Fibonacci number using faster recursion
function F = fibonacci_recursive2(n) if n==0
F = 0;
elseif n==1
F = 1;
elseif n==2
F = 1;
else
if mod(n,2)
F = fibonacci_recursive2((n-1)/2)ˆ2+fibonacci_recursive2((n+1)/2)ˆ2; else
F = fibonacci_recursive2(n/2+1)ˆ2-fibonacci_recursive2(n/2-1)ˆ2; end
end
Evaluation by rounding explicit expression
An explicit expression is F = 1 􏰛􏰒1+√5 􏰓n − 􏰒1−√5 􏰓n􏰜 n√52 2
m×2
(m − 1) ×2
Proof. Solutions of the recurrence of the form Fn = rn must satisfy r2 = r + 1 (characteristic equation), i.e., r = r1 = 1+√5 ≈ 2
1.618 (“golden ratio”) or r = r2 = 1−√5 ≈ −0.618. Then Fn = arn + brn with F0 = 0 = a + b and F1 = 1 = ar1 + br2, i.e., 212
a = 1 = −b, which leads to the given formula. 􏰞 r1 −r2
This formula does not evaluate to an integer for larger values of n due to round-off errors but can still be used to obtain Fn, see P5.2.4 . See the Wikipedia page for additional properties and occurences in nature of Fibonacci numbers.
􏰑BDW Draft August 22, 2013

page 5.8 􏰝5.3. POLYNOMIAL EVALUATION P5.2.1 Show that the number Wn of additions needed to determine Fn with fibonacci recursive
for0≤k≤n.
n n
exceeds Fn (i.e., growth exponentially with n).
P5.2.2 ThebinomialcoefficientCn isCn =􏰉n􏰈= n!
k k k k!(n−k)! a) Show that Cn = Cn−1 + Cn−1 for 0 < k < n. k k k−1 b) Write a Matlab function implementing this recurrence (initialize C0 = Cn = 1). c) Use it to evaluate 􏰉20􏰈. 10 P5.2.3 The sequence {an}n≥0 defined by a0 = 3, a1 = 2, and an = 2an−1 − 8an−2 for n > 1.
a) Determine an explicit expression for an in terms of n.
b) Write a Matlab function in the style of fibonacci direct2 to compute an. c) Compare the explicit and computed values of a40; of a60; what happens?
P5.2.4 Show that F can be obtained by rounding the quantity 1 􏰒1+√5􏰓n to the nearest n √5 2
integer for any n. Compare F70 obtained this way with the values obtained using the functions described in this section.
5.3 Polynomial evaluation
9
The polynomial
can be na ̈ıvely evaluated using peval.m evaluates p(x)
p(x)=a0 +a1x+…+anxn
function p = peval(a,x)
% naive evaluation of polynomial p = a(1)+a(2)*x+…+a(n+1)*xˆn n = length(a)-1;
p = 0*x+a(1);
for k = 1:n
p = p + a(k+1)*x.ˆk; % coefficientwise multiplication requires . if x is array end
The number of operations performed is
Wn =n+1+2+…+n=n+n(n+1)=n(n+3)=O(n2) 22
Horner’s rule
A more efficient way to evaluate p(x) is to write it as a nested product p(x) = (· · · (((an)x + an−1)x + an−2)x + · · · )x + a0
and evaluate it starting with an and working the way from inside-out. The resulting method is called Horner’s rule.
􏰑BDW Draft August 22, 2013

􏰝5.3. POLYNOMIAL EVALUATION horner.m evaluates p(x) with Horner’s rule
function p = horner(a,x)
% Horner’s rule to evaluate p = a(1)+a(2)*x+…+a(n+1)*xˆn
n = length(a)-1;
p = 0*x+a(n+1); % adding 0*x gives p the same size as input x for k = n:-1:1
p = p.*x+a(k); % coefficientwise multiplication requires . if x is array end
The work associated to Horner’s rule is
Wn = 2n = O(n),
i.e., much less than the na ̈ıve implementation.
page 5.9
MATLAB’s polyval function
Matlab includes its own function to evaluate polynomials. A slight difference is that the
coefficient array a must be passed in reverse order (so that a(1) = an is used first).
100000
􏰃 xk 1 1
Example 5.1 Evaluate p(x) = k2 at x = −1, − 2 , 2 , 1 with peval, horner, and polyval: k=1
>> format long
>> a = [0,1./(1:100000).ˆ2]; % a_0 = 0 >> x = [-1,-.5,.5,1];
>> tic, peval(a,x), toc
ans =
-0.822467033374095 -0.448414206923646 Elapsed time is 0.082717 seconds.
>> tic, horner(a,x), toc
ans =
-0.822467033374114 -0.448414206923646 Elapsed time is 0.030736 seconds.
>> tic, polyval(a(end:-1:1),x), toc
ans =
-0.822467033374114 -0.448414206923646 Elapsed time is 0.031449 seconds.
0.582240526465012
0.582240526465013
0.582240526465013
1.644924066898242
1.644924066898226
1.644924066898226
Horner’s rule is more efficient. Note that polyval yields similar timing and identical values! The values obtained with peval at x = ±1 are somewhat different. Since1
􏰃(−1)k π2
􏰃1 π2
k2 = −12 ≈ −0.822467033424113 and k≥1
k2 = 6 ≈ 1.644934066848226 k≥1
1The second series is the value ζ(2) of the Riemann zeta function, and was found by Euler (1737) by comparing powers of x3 in
x3 x5 􏰔 x2􏰕􏰔 x2 􏰕
(Taylor) x− 3! + 5! −…=sinx=x 1−π2 1−(2π)2 ··· (sincesinvanishesat0,±π,±2π,…)
The relation between the first and second series can be derived without knowing their value, see P5.3.2 . 􏰑BDW Draft August 22, 2013

page 5.10 􏰝5.4. BISECTION FOR NONLINEAR EQUATION Horner’s rule is more accurate here. This is not surprising since small numbers are added
together first, compared to larger values of ak added first in peval.
P5.3.1 Write a recursive version horner recursive.m of horner.m. Test it on the polynomial p and x values from Example 5.1, but with 1, 000 terms instead of 100, 000.
􏰃 (−1)k P5.3.2 Show that k2
k≥1
5.4 Bisection for nonlinear equation
Bisection is a numerical technique used to determine roots of nonlinear equations in 1D by repeatedly dividing an interval containing a root until the size of the interval is small enough.
bisect.m solves f(x) = 0 using bisection
1 function ab = bisect(f,ab,tol,alfa)
2 a = ab(1); b = ab(2); % assumes a> format long
>> f = @(x)xˆ2-2;
>> ab = [1,2]; % assumes fa*fb<0 >> tol = 1e-15;
>> tic, bisect(f,ab,tol,0.5), toc ans =
2
1
√2􏰟
1.414213562373095 1.414213562373096
Elapsed time is 0.002212 seconds.
>> tic, bisect(f,ab,tol,0.05), toc
ans = 0
1.414213562373095 1.414213562373095 Elapsed time is 0.000828 seconds.
>> tic, bisect(f,ab,tol,0.95), toc
ans =
1.414213562373095 1.414213562373095 Elapsed time is 0.000903 seconds.
>> tic, bisect(f,ab,tol,Inf), toc
ans =
1.414213562373095 1.414213562373095 Elapsed time is 0.000691 seconds.
-1
-2
012
􏰑BDW
Draft August 22, 2013