Introduction to Scientific
Computing
-- A Matrix Vector Approach Using Matlab
Written by Charles F.Van Loan
陈 文 斌
Wbchen@fudan.edu.cn
复旦大学
网上的材料
? http://www.cs.cornell.edu/cv
? http://www.netlib.org/index.html
? http://www.scicomput.com
Chapter 1Power tools of the trade
1,Vectors and Plotting
2,More Vectors,More Plotting and New
Matrices
3,Building Exploratory Environments
4,Error
5,Designing Functions
Six
windows
intro
demo
Vectors and Plotting
Setting Up Vectors
Setting Up Vectors
? X=[10.1 20.2 30.3]
? X=[10.1; 20.2; 30.3]
? X=[10.1; 20.2; 30.3];
? X=[10.1 20.2 30.3]'
? X=[0,05,10,15,20,25,30,35,40,45,50 …
.55,60,65,70,75,80,85,90,95 1.0]
n=21;
h=1/(n-1);
for k=1:n
x(k)=(k-1)*h;
end
n=21;
h=1/(n-1);
x=zeros(1,n);
for k=1:n
x(k)=(k-1)*h;
end
1 "Dimension
mismatch" error
2,Efficient(memory)
In Matlab,variables are not declared by the
user but are created on a need-to-use basis
by a memory manager,Moreover,from
Matlab's point of view,every variable is a
complex matrix indexed from unity.
length
? u=[10 20 30];
? n=length(u);
? v=[10;20;30;40];
? m=length(v);
? u=[50 60];
? p=length(u);
help length
LENGTH Length of vector.
LENGTH(X) returns the length
of vector X,It is equivalent to
MAX(SIZE(X)) for non-empty
arrays and 0 for empty ones.
help
? help who
? help whos
? help clear
? help for
? help zeros
? help ;
? help []
help help
Regular vectors
? x=20:24
? x=20:2:29;
? x=10:-1:1
? x=3:2
colon notation
<Starting Index>:<stride>:<Bounding Index>
linspace
? x=linspace(a,b,n)
? xk=a+(k-1)*(b-a)/(n-1)
? x=linspace(0,1,10) 0 0.1111 0.2222 0.3333
0.4444 0.5556 0.6667
0.7778 0.8889 1.0000
Linspace(Left Endpoint,Right Endpoint,Number of Points)
logspace
? x=logspace(-2,2,6)
? x=logspace(a,b,n)
0.0100 0.0631 0.3981
2.5119 15.8489 100.0000
m=linspace(a,b,n)
for k=1:n
x(k)=10^m(k);
end
nkx nkabak,.,,,1,10 )1/()1)(( ?? ????
linspace(a,b) linspace(a,b,100)
logspace(a,b) logspace(a,b,50)
format
? x=.123456789012345*logspace(1,5,5)'
Format
short
Format
long
Format
short e
Format
long e
x?
Evaluating Functions
n=21;
x=linspace(0,1,n);
y=zeros(1,n);
for k=1:n
y(k)=sin(x(k));
end
n=21;
x=linspace(0,1,n);
y=sin(x);
vectorization
help elfun
Vectorization
? Speed
? Clarity
? Education
0 5 10 15 20 25-1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
Example
m=5;
n=4*m+1;
x=linspace(0,1,n);
y=zeros(1,n);
a=x(1:m+1);
y(1:m+1) = sin(2*pi*a);
y(2*m+1:-1:m+2)= y(1:m);
y(2*m+2:n)=-y(2:2*m+1);
y=sin(2*pi*linspace(0,1,21))
% SineTable,Prints a short table of sine evaluations,
clc
n = 21;
x = linspace(0,1,n);
y = sin(2*pi*x);
disp(' ')
disp(' k x(k) sin(x(k))')
disp('------------------------')
for k=1:21
degrees = (k-1)*360/(n-1);
disp(sprintf(' %2.0f %3.0f %6.3f ',k,degrees,y(k)));
end
disp( ' ');
disp('x(k) is given in degrees.')
disp(sprintf('One Degree = %5.3e Radians',pi/180))
Displaying Tables--sprintf
Sinetable.m
Help sinetable
A Simple Plot
n=20; %n=200;
x=linspace(0,1,n);
y=sin(2*pi*x);
hold on
plot(x,y) %plot(x,y,'r')
title('The function y=sin(2*pi*x)')
xlabel('x (in radians)')
ylabel('y') 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
The function y=sin(2*pi*x)
x (in radians)
y
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
The function y=sin(2*pi*x)
x (in radians)
y
n=
% Script File,SinePlot
% Displays increasingly smooth plots of sin(2*pi*x),
close all
pause(2)
for n = [4 8 12 16 20 50 100 200 400]
x = linspace(0,1,n);
y = sin(2*pi*x);
plot(x,y)
title(sprintf('Plot of sin(2*pi*x) based …
upon n = %3.0f points.',n))
pause(2)
end
homework
? P1.1.1
? P1.1.2
Vectorizing Function Evaluations
n=200;
x=linspace(0,1,n);
y=zeros(1,n);
for k=1:n
y(k)=((1-x(k)/24)/(1+x(k)/24+(x(k)/384)*x(k)))^8
end
plot(x,y)
x
e
xx
x
xf
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
?
?
8
2
3 8 424
1
24
1
)(
% Script File,ExpPlot
% Examines the function
% f(x) = ((1 - x/24)/(1 + x/24 + x^2/384))^8
% as an approximation to exp(-x) across [0,1].
x=linspace(0,1,200);
y=((1-x/24)./(1+x/24+(x/384).*x)).^8
plot(x,y)
title('(1-x/24)/(1+x/24+x^2/384))^8')
.^ pointwise vector exponentiation
./ pointwise vector divide
.* pointwise vector multiple
Scaling and superpositioning
x=linspace(-pi/2,11*pi/2,200);
y=tan(x);
plot(x,y)
axis([-pi/2 9*pi/2 -10 10])
-2 0 2 4 6 8 10 12 14 16 18
- 1 8
- 1 6
- 1 4
- 1 2
- 1 0
-8
-6
-4
-2
0
2
x 1 0
15
0 2 4 6 8 10 12 14
- 1 0
-8
-6
-4
-2
0
2
4
6
8
10
autoscaling
Axis([xmin xmax ymin ymax])
% Plots the function tan(x),-pi/2 <= x <= 9pi/2
close all;ymax = 10;
x = linspace(-pi/2,pi/2,40);
y = tan(x); plot(x,y)
axis([-pi/2 9*pi/2 -ymax ymax])
title('The Tangent Function');xlabel('x'); ylabel('tan(x)')
hold on
for k=1:4
xnew = x+ k*pi;
plot(xnew,y);
End
hold off 0 2 4 6 8 10 12 14-10
-8
-6
-4
-2
0
2
4
6
8
10
The Tangent Function
x
tan
(x
)
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
x=linspace(0,1,200);
y1=sin(2*pi*x); y2=cos(2*pi*x);
plot(x,y1); hold on;
plot(x,y2,'--');
plot([1/8 5/8],[1/sqrt(2) -1/sqrt(2)],'*')
hold off
x=linspace(0,1,200);
y1=sin(2*pi*x); y2=cos(2*pi*x);
plot(x,y1,x,y2,'--',[1/8 5/8],[1/sqrt(2) -1/sqrt(2)],'*')
Plot(….)
% Plots selected regular polygons.
close all;clc
theta = linspace(0,2*pi,361);
c = cos(theta);s = sin(theta);k=0;
for sides = [3 4 5 6 8 10 12 18 24]
stride = 360/sides;
k=k+1;
subplot(3,3,k)
plot(c(1:stride:361),s(1:stride:361))
axis([-1.2 1.2 -1.2 1.2])
axis equal
end
interesting
Subplot(m,n,k):Split the current
figure into an m-row by n
column array of subwindows
that are indexed left to right,top
to bottom
-1 0 1
-1
0
1
-1 0 1
-1
0
1
-1 0 1
-1
0
1
-1 0 1
-1
0
1
-1 0 1
-1
0
1
-1 0 1
-1
0
1
-1 0 1
-1
0
1
-1 0 1
-1
0
1
-1 0 1
-1
0
1
-10 -8 -6 -4 -2 0 2 4 6 8 10-20
-15
-10
-5
0
5
10
15
20
f(x)=2sin(x)+3sin(2x)+7sin(3x)+5sin(4x)
Some matrix computations
n=200;
x=linspace(-10,10,n)';
y=zeros(n,1);
for k=1:n
y(k)=2*sin(x(k))+3*sin(2*x(k))+7*sin(3*x(k))+5*sin(4*x(k));
end
plot(x,y)
title('f(x)=2sin(x)+3sin(2x)+7sin(3x)+5sin(4x)')
Vectorization
n=200;
x=linspace(-10,10,n)';
y(k)=2*sin(x(k))+3*sin(2*x(k))+7*sin(3*x(k))+5*sin(4*x(k));
plot(x,y)
title('f(x)=2sin(x)+3sin(2x)+7sin(3x)+5sin(4x)')
-10 -8 -6 -4 -2 0 2 4 6 8 10-20
-15
-10
-5
0
5
10
15
20
f(x)=2sin(x)+3sin(2x)+7sin(3x)+5sin(4x)
?
?
?
?
?
?
?
?
?
?
?
?
5
7
3
2
]4s i n,3s i n,2s i n,[ s in xxxx
)4s i n (5)3s i n (7)2s i n (3)s i n (2)( xxxxxf ????
By vector
close all
x = linspace(-10,10,200)';
A = [sin(x) sin(2*x) sin(3*x) sin(4*x)];
y = A*[2;3;7;5];
plot(x,y)
title('f(x) = 2sin(x) + 3sin(2x) + 7sin(3x) +
5sin(4x)')
for j=1:m
for k=1:n
A(k,j)=sin(j*x(k))
end
end
for j=1:m
A(:,j)=sin(j*x)
end
for k=1:n
A(k,:)=sin(k*(1:m))
end
)4s in (9)3s in (6)2s in (2)s in (8)(
)4s in (5)3s in (7)2s in (3)s in (2)(
xxxxxg
xxxxxf
????
????
n=200;x = linspace(-10,10,n)';
A = [sin(x) sin(2*x) sin(3*x) sin(4*x)];
y1 = A*[2;3;7;5]; y2=A*[8; 2; 6; 9];
plot(x,y1,x,y2)
-10 -8 -6 -4 -2 0 2 4 6 8 10-20
-15
-10
-5
0
5
10
15
20
y=A*[2 8;3 2;7 6;5 9];
plot(x,y(:,1),x,y(:,2))
plot(x,y)
By matrix
,,,
5)4(,1)3(,0)2(,2)1(
)4s in ()3s in ()2s in ()s in ()(
4321
4321
?
?????
????
????
????
ffff
xxxxxf
X=[1 2 3 4;2 4 6 8; 3 6 9 12;4 8 12 16];
Z=sin(x);
F=[-2;0;1;5];
alpha=Z\F
backslash
homework
? P1.2.1-P1.2.10
Building Exploratory
Environments
Random Process
-1 -0.5 0 0.5 1 1.5 20
10
20
30
40
50
60
Distribution of Values in rand(1000,1)
Mean = 0.502,Median = 0.510.
-3 -2 -1 0 1 2 30
10
20
30
40
Distribution of Values in randn(1000,1)
Mean = -0.043,Standard Deviation = 0.943
close all
subplot(2,1,1); x = rand(1000,1);
hist(x,30); axis([-1 2 0 60])
title('Distribution of Values in rand(1000,1)')
xlabel(sprintf('Mean = %5.3f,Median = …
%5.3f.',mean(x),median(x)))
subplot(2,1,2); x = randn(1000,1);
hist(x,linspace(-2.9,2.9,100));
title('Distribution of Values in randn(1000,1)')
xlabel(sprintf('Mean = %5.3f,Standard Deviation =…
%5.3f',mean(x),std(x)))
Clouds.m
close all; Points = rand(1000,2);
subplot(1,2,1)
plot(Points(:,1),Points(:,2),'.')
title('Uniform Distribution.')
axis([0 1 0 1])
axis square
Points = randn(1000,2);
subplot(1,2,2)
plot(Points(:,1),Points(:,2),'.')
title('Normal Distribution.')
axis([-3 3 -3 3])
axis square
2 3 4 5 6 7 8 9 10 11 12
0
20
40
60
80
100
120
140
160
180
Outcome of 1000 Dice Rolls.
% Script File,Dice
% Simulates 1000 rollings of a pair of dice.
close all
First = 1 + floor(6*rand(1000,1));
Second = 1 + floor(6*rand(1000,1));
Throws = First + Second;
hist(Throws,linspace(2,12,11));
title('Outcome of 1000 Dice Rolls.')
First = ceil(6*rand(1000,1));
-1.5 -1 -0.5 0 0.5 1 1.5
-1.5
-1
-0.5
0
0.5
1
1.5
p's estimate by Monte Carlo
Hundreds of Trials
0 50 100 150 200 250 300 350 400 450 500
3.08
3.1
3.12
3.14
3.16
3.18
3.2
3.22
3.24
3.26
Monte Carlo Estimate of Pi = 3.157
Hundreds of Trials
Darts.m
close all; rand('seed',.123456)
NumberInside = 0;
PiEstimate = zeros(500,1);
for k=1:500
x = -1+2*rand(100,1); y = -1+2*rand(100,1);
NumberInside = NumberInside + sum(x.^2 + y.^2 <= 1);
PiEstimate(k) = (NumberInside/(k*100))*4;
end
plot(PiEstimate)
title(sprintf('Monte Carlo Estimate of Pi = %5.3f',PiEstimate(500)));
xlabel('Hundreds of Trials')
Monte Carlo method
Polygon Smoothing
Points (xi,yi ); i=1:n+1
Plot(x,y,x,y,’*’)
xnew = [(x(1:n)+x(2:n+1))/2; (x(1)+x(2))/2];
ynew = [(y(1:n)+y(2:n+1))/2;(y(1)+y(2))/2];
plot(xnew,ynew)
A new polygon is displayed that is obtained by connecting
the side midpoints of the original polygon.
close all
n = input('Enter the number of edges:');figure; axis([0 1 0 1])
axis square;hold on
x = zeros(n,1); y = zeros(n,1);
for k=1:n
title(sprintf('Click in %2.0f more points.',n-k+1))
[x(k) y(k)] = ginput(1);
plot(x(1:k),y(1:k),x(1:k),y(1:k),'*')
end
x0 = x; y0 = y;x = [x;x(1)];y = [y;y(1)];plot(x,y,x,y,'*')
title('The Original Polygon') ;hold off;k=0;
xlabel('Click inside window to smooth,outside window to quit.')
[a,b] = ginput(1); v = axis;
while (v(1)<=a) & (a<=v(2)) & (v(3)<=b) & (b<=v(4));
k = k+1;
x = [(x(1:n)+x(2:n+1))/2;(x(1)+x(2))/2];
y = [(y(1:n)+y(2:n+1))/2;(y(1)+y(2))/2];
clf; plot(x,y,x,y,'*') ;axis([0 1 0 1]);
axis square
title(sprintf('Number of Smoothings = %1.0f',k))
xlabel('Click inside window to smooth,outside window to quit.')
v = axis;
[a,b] = ginput(1);
end
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
The Original Polygon
Click inside window to smooth,outside window to quit.
Smooth.m
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Number of Smoothings = 1
Click inside window to smooth,outside window to quit.
Error
? Focus on the mathematical errors that arise
through discretization and the rounding
errors that arise due to finite precision
arithmetic
Absolute and Relative Error
? Absolute error
? Relative error
? If the relative error is 10-d,then has
approximately d correct significant digits
|~| xx ?
||/|~| xxx ?
x~
% Script File,Stirling
% Prints a table showing error in Stirling's formula for n!
disp(' Stirling Absolute Relative')
disp(' n n! Approximation Error Error')
disp('----------------------------------------------------------------')
e=exp(1);nfact=1;
for n=1:13
nfact = n*nfact;
s = sqrt(2*pi*n)*((n/e)^n);
abserror = abs(nfact - s);
relerror = abserror/nfact;
s1 = sprintf(' %2.0f %10.0f %13.2f',n,nfact,s);
s2 = sprintf(' %13.2f %5.2e',abserror,relerror);
disp([s1 s2])
end
!2 n
e
nnS n
n ??
?
??
?
?? p
Stirling
Stirling Absolute Relative
n n! Approximation Error Error
1 1 0.92 0.08 7.79e-002
2 2 1.92 0.08 4.05e-002
3 6 5.84 0.16 2.73e-002
4 24 23.51 0.49 2.06e-002
5 120 118.02 1.98 1.65e-002
6 720 710.08 9.92 1.38e-002
7 5040 4980.40 59.60 1.18e-002
8 40320 39902.40 417.60 1.04e-002
9 362880 359536.87 3343.13 9.21e-003
10 3628800 3598695.62 30104.38 8.30e-003
11 39916800 39615625.05 301174.95 7.55e-003
12 479001600 475687486.47 3314113.53 6.92e-003
13 6227020800 6187239475.19 39781324.81 6.39e-003
xx
n
e
k
xe n
k
n
k
x ??
?
?? ?
?
? ?
?
0,
)!1(!0
1
% Script File,ExpTaylor Plots,as a function of n,the relative error in
%the Taylor approximation 1 + x + x^2/2! +...+ x^n/n! to exp(x).
nTerms = 50;
for x=[10 5 1 -1 -5 -10]
figure; f = exp(x)*ones(nTerms,1); s = 1; term = 1;
for k=1:50
term = x.*term/k; s = s+ term; err(k) = abs(f(k) - s);
end
relerr = err/exp(x); semilogy(1:nTerms,relerr)
ylabel('Relative Error in Partial Sum.');
xlabel('Order of Partial Sum.') title(sprintf('x = %5.2f',x))
end
0 5 10 15 20 25 30 35 40 45 50
10
-16
10
-14
10
-12
10
-10
10
-8
10
-6
10
-4
10
-2
10
0
R
e
la
t
iv
e
E
r
r
o
r
in
P
a
r
t
ia
l
S
u
m
.
O r d e r o f P a r t i a l S u m,
x = 1 0, 0 0
0 5 10 15 20 25 30 35 40 45 50
10
-10
10
-5
10
0
10
5
10
10
R
e
la
t
iv
e
E
r
r
o
r
in
P
a
r
t
ia
l
S
u
m
.
O r d e r o f P a r t i a l S u m,
x = - 1 0, 0 0
0 5 10 15 20 25 30 35 40 45 50
10
-16
10
-14
10
-12
10
-10
10
-8
10
-6
10
-4
10
-2
10
0
R
e
la
t
iv
e
E
r
r
o
r
in
P
a
r
t
ia
l
S
u
m
.
O r d e r o f P a r t i a l S u m,
x = 5, 0 0
0 2 4 6 8 10 12 14 16
10
-16
10
-14
10
-12
10
-10
10
-8
10
-6
10
-4
10
-2
10
0
R
e
la
t
iv
e
E
r
r
o
r
in
P
a
r
t
ia
l
S
u
m
.
O r d e r o f P a r t i a l S u m,
x = 1, 0 0
bottom
out
10-16
10-10
10-16 10-16
x=10
Relative error
x=-10
x=5 x=1
Rounding Errors
? The plots reveal that the mathematical
convergence theory does not quit apply
? The errors do not go to zero as the number of
terms in the series increases.
? The relative error for x=-10 is much worse that for
x=10
? Floating point arithmetic!
mathematics Development of algorithms
Inexact computer
arithmetic system
% Script File,Zoom
% Plots (x-1)^6 near x=1 with increasingly refined scale.
% Evaluation via x^6 - 6x^5 + 15x^4 - 20x^3 + 15x^2 - 6x +1
% leads to severe cancellation.
close all;k=0;n=100;
for delta = [.1,01,008,007,005,003 ]
x = linspace(1-delta,1+delta,n)';
y = x.^6 - 6*x.^5 + 15*x.^4 - 20*x.^3 + 15*x.^2 - 6*x +
ones(n,1);
k = k+1;
subplot(2,3,k); plot(x,y,x,zeros(1,n))
axis([1-delta 1+delta -max(abs(y)) max(abs(y))])
end
0.9 1 1.1
-5
0
5
x 10 -7
0.99 1 1.01
-1
-0.5
0
0.5
1 x 10
-12
0.995 1 1.005
-2
-1
0
1
2
x 10 -13
0.995 1 1.005
-1
-0.5
0
0.5
1
x 10 -13
0.995 1 1.005
-1
0
1
x 10 -14
0.998 1 1.002
-5
0
5
x 10 -15
y = x.^6 - 6*x.^5 + 15*x.^4 - 20*x.^3 + 15*x.^2 -
6*x + ones(n,1);
One
zero?
0.9 1 1.1
-1
-0.5
0
0.5
1 x 10
-6
0.99 1 1.01
-1
-0.5
0
0.5
1 x 10
-12
0.995 1 1.005
-2
-1
0
1
2
x 10 -13
0.995 1 1.005
-1
-0.5
0
0.5
1
x 10 -13
0.995 1 1.005
-1
0
1
x 10 -14
0.998 1 1.002
-5
0
5
x 10 -16
y = (x-1).^6
Algorithms that are equivalent mathematically may
behave very differently numerically.
The floating point number
x=1; p=0; y=1; z=x+y;
while x~=z
y=y/2;
p=p+1;
z=x+y;
end
disp(sprintf('p = %2.0f is the smallest
positive integer so 1+1/2^p = 1.',p))
12/11 ?? p
P=53
% q = smallest positive integer so 1/2^q = 0.
x=1; q=0;
while x>0
x=x/2;
q=q+1;
end;
02/1 ?q
q=1075
% r = smallest positive integer so 2^r = inf.
x=1; r=0;
while x~=inf
x=2*x; r=r+1;
end
in f2 ?r
r=1024
p<q?
underflow
overflow
% The Kahan example,In exact arithmetic f/e = 0/0:
disp(' ')
h = 1./2.;
x = 2./3,- h;
y = 3./5,- h;
e = (x+x+x) - h;
f = (y+y+y+y+y)-h;
z = f/e;
disp(sprintf('\nz = %10.5f',z))
z=1
Designing function
? An ability to write good Matlab function is
crucial.
? Three examples are used to clarify the
essential ideas,
– Taylor series
– Numerical differentiation
– Three-digit Arithmetic(see book)
Four ways to compute the
exponential of a vector
function y=MyExpF(x,n)
%Pre:
% x=scalar
% n=positive integer
%Post:
% y=n-th order Taylor approximation to exp(x)
y=1;
term=1;
for k=1:n
term=x*term/k;
y=y+term;
end
?
?
?
n
k
k
n k
xxT
0 !
)(
Help MyExpF
All variables inside the
function are local and
are not part of the
Matlab workspace
The input and output parameters are
formal parameters
close all;m = 50;x = linspace(-1,1,m);
exact = exp(x);figure;k=0;
y = zeros(1,m);
for n = [4 8 16 20]
for i=1:m
y(i) = MyExpF(x(i),n);
end
RelErr = abs(exact - y)./exact;
k=k+1;
subplot(2,2,k); plot(x,RelErr)
title(sprintf('n = %2.0f',n))
end
-1 -0.5 0 0.5 10
0.005
0.01
0.015
0.02
n = 4
-1 -0.5 0 0.5 10
2
4
6
8x 10
-6 n = 8
-1 -0.5 0 0.5 10
2
4
6
8x 10
-15 n = 16
-1 -0.5 0 0.5 10
2
4
6x 10
-16 n = 20
10-8
10-15 10-16
function y=MyExp1(x)
n=17;
p=length(x);
y=ones(p,1);
for i=1:p
y(i)=MyExpF(x(i),n);
end
function y=MyExp2(x)
y=ones(size(x));
n=17;
term=ones(size(x));
for k=1:n
term=x.*term/k;
y=y+term;
end
Length(x) MyExp1(x) MyExp2(x)
100 0.02 0.00
200 0.03 0.00
300 0.03 0.00
Run time
Faster than
result in book
disp(' Length(x) MyExp1(x) MyExp2(x) ')
disp(' Time Time')
disp('-----------------------------------------')
for L=50:50:300
xL=linspace(-1,1,L);
tstart=clock; y=MyExp1(xL); tfinish=clock;
t1=etime(tfinish,tstart);
tstart=clock; y=MyExp2(xL); tfinish-clock;
t2=etime(tfinish,tstart);
disp(sprintf('%6.0f %13.2f %13.2f ',L,t1,t2))
end
0 500 1000 1500 2000 2500 30000
0.05
0.1
0.15
0.2
0.25
0.3
0.35Length(x) MyExp1(x) MyExp2(x)
450 0.0470 0.0000000000000
500 0.0470 0.0160000000000
550 0.0620 0.0000000000000
600 0.0630 0.0000000000000
650 0.0620 0.0160000000000
700 0.0780 0.0000000000000
750 0.0780 0.0000000000000
800 0.0940 0.0000000000000
850 0.0940 0.0000000000000
900 0.0930 0.0000000000000
950 0.1100 0.0000000000000
function y = MyExpW(x)
% y = MyExpW(x)
% x is a scalar and y is a Taylor approximation to exp(x).
y = 0;
term = 1;
k=0;
while abs(term) > eps*abs(y)
k = k + 1;
y = y + term;
term = x*term/k;
end
function y = MyExp3(x)
% y = MyExp3(x)
% x is a column n-vector and y is a column n-vector
with the property that
% y(i) is a Taylor approximation to exp(x(i)) for i=1:n.
n = length(x);
y = ones(n,1);
for i=1:n
y(i) = MyExpW(x(i));
end
function y = MyExp4(x)
% y = MyExp4(x)
% x is an n-vector and y is an n-vector with the same shape and
%the property that y(i) is a Taylor approximation to exp(x(i)) for
i=1:n.
y = zeros(size(x));
term = ones(size(x));
k = 0;
while any(abs(term) > eps*abs(y))
y = y + term;
k = k+1;
term = x.*term/k;
end
Numerical Differentiatioin
2
2
)('')(')()( hfhafafhaf ?????
h
afhafD
h
)()( ???
2)('')('
hfafD
h ???
a=input('Enter a,');
h=logspace(-1,-16,16);
Dh=(sin(a+h)-sin(a))./h;
err=abs(Dh-cos(a));
0 2 4 6 8 10 12 14 16-20
-18
-16
-14
-12
-10
-8
-6
-4
-2
0
1.0000e-001
1.0000e-002
1.0000e-003
1.0000e-004
1.0000e-005
1.0000e-006
1.0000e-007
1.0000e-008
1.0000e-009
1.0000e-010
1.0000e-011
1.0000e-012
1.0000e-013
1.0000e-014
1.0000e-015
1.0000e-016
0.04293855333275
0.00421632485627
0.00042082550781
0.00004207444952
0.00000420736228
0.00000042074681
0.00000004182769
0.00000000296989
0.00000005254127
0.00000005848104
0.00000116870406
0.00004324021692
0.00073391590031
0.00370697619819
0.01480920644444
0.54030230586814
loss of accuracy
Explanation
|)(''|
2
|)('| ?fhafD h ??
h
fhafD h e p s2|)(''|
2
|)('| ??? ?
Truncation
error
|)(''|e p s /2 ?fh ?
minimize
The choice of h
2|)(''| Mxf ?
h
hMhe r r D ?2
2
)( 2 ??
2
2
M
h opt
?
? 22)( Mhe r r D opt ??
function d=Derivative(fname,a,delta,M2)
if nargin<=3
M2=10;
end
if nargin==2
delta=eps;
end
hopt=2*sqrt(delta/M2);
d=(feval(fname,a+hopt)-feval(fname,a))/hopt;
NARGIN returns the number of input
arguments that were used to call the function,
FEVAL(F,x1,...,xn) evaluates
the function
da=derivative(fname,a,delta,M2)
da=derivative(fname,a,delta)
da=derivative(fname,a)
da=derivative('sin',1,eps,1)
da=derivative('sin',1,eps)
da=derivative('sin',1)
Computing
-- A Matrix Vector Approach Using Matlab
Written by Charles F.Van Loan
陈 文 斌
Wbchen@fudan.edu.cn
复旦大学
网上的材料
? http://www.cs.cornell.edu/cv
? http://www.netlib.org/index.html
? http://www.scicomput.com
Chapter 1Power tools of the trade
1,Vectors and Plotting
2,More Vectors,More Plotting and New
Matrices
3,Building Exploratory Environments
4,Error
5,Designing Functions
Six
windows
intro
demo
Vectors and Plotting
Setting Up Vectors
Setting Up Vectors
? X=[10.1 20.2 30.3]
? X=[10.1; 20.2; 30.3]
? X=[10.1; 20.2; 30.3];
? X=[10.1 20.2 30.3]'
? X=[0,05,10,15,20,25,30,35,40,45,50 …
.55,60,65,70,75,80,85,90,95 1.0]
n=21;
h=1/(n-1);
for k=1:n
x(k)=(k-1)*h;
end
n=21;
h=1/(n-1);
x=zeros(1,n);
for k=1:n
x(k)=(k-1)*h;
end
1 "Dimension
mismatch" error
2,Efficient(memory)
In Matlab,variables are not declared by the
user but are created on a need-to-use basis
by a memory manager,Moreover,from
Matlab's point of view,every variable is a
complex matrix indexed from unity.
length
? u=[10 20 30];
? n=length(u);
? v=[10;20;30;40];
? m=length(v);
? u=[50 60];
? p=length(u);
help length
LENGTH Length of vector.
LENGTH(X) returns the length
of vector X,It is equivalent to
MAX(SIZE(X)) for non-empty
arrays and 0 for empty ones.
help
? help who
? help whos
? help clear
? help for
? help zeros
? help ;
? help []
help help
Regular vectors
? x=20:24
? x=20:2:29;
? x=10:-1:1
? x=3:2
colon notation
<Starting Index>:<stride>:<Bounding Index>
linspace
? x=linspace(a,b,n)
? xk=a+(k-1)*(b-a)/(n-1)
? x=linspace(0,1,10) 0 0.1111 0.2222 0.3333
0.4444 0.5556 0.6667
0.7778 0.8889 1.0000
Linspace(Left Endpoint,Right Endpoint,Number of Points)
logspace
? x=logspace(-2,2,6)
? x=logspace(a,b,n)
0.0100 0.0631 0.3981
2.5119 15.8489 100.0000
m=linspace(a,b,n)
for k=1:n
x(k)=10^m(k);
end
nkx nkabak,.,,,1,10 )1/()1)(( ?? ????
linspace(a,b) linspace(a,b,100)
logspace(a,b) logspace(a,b,50)
format
? x=.123456789012345*logspace(1,5,5)'
Format
short
Format
long
Format
short e
Format
long e
x?
Evaluating Functions
n=21;
x=linspace(0,1,n);
y=zeros(1,n);
for k=1:n
y(k)=sin(x(k));
end
n=21;
x=linspace(0,1,n);
y=sin(x);
vectorization
help elfun
Vectorization
? Speed
? Clarity
? Education
0 5 10 15 20 25-1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
Example
m=5;
n=4*m+1;
x=linspace(0,1,n);
y=zeros(1,n);
a=x(1:m+1);
y(1:m+1) = sin(2*pi*a);
y(2*m+1:-1:m+2)= y(1:m);
y(2*m+2:n)=-y(2:2*m+1);
y=sin(2*pi*linspace(0,1,21))
% SineTable,Prints a short table of sine evaluations,
clc
n = 21;
x = linspace(0,1,n);
y = sin(2*pi*x);
disp(' ')
disp(' k x(k) sin(x(k))')
disp('------------------------')
for k=1:21
degrees = (k-1)*360/(n-1);
disp(sprintf(' %2.0f %3.0f %6.3f ',k,degrees,y(k)));
end
disp( ' ');
disp('x(k) is given in degrees.')
disp(sprintf('One Degree = %5.3e Radians',pi/180))
Displaying Tables--sprintf
Sinetable.m
Help sinetable
A Simple Plot
n=20; %n=200;
x=linspace(0,1,n);
y=sin(2*pi*x);
hold on
plot(x,y) %plot(x,y,'r')
title('The function y=sin(2*pi*x)')
xlabel('x (in radians)')
ylabel('y') 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
The function y=sin(2*pi*x)
x (in radians)
y
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
The function y=sin(2*pi*x)
x (in radians)
y
n=
% Script File,SinePlot
% Displays increasingly smooth plots of sin(2*pi*x),
close all
pause(2)
for n = [4 8 12 16 20 50 100 200 400]
x = linspace(0,1,n);
y = sin(2*pi*x);
plot(x,y)
title(sprintf('Plot of sin(2*pi*x) based …
upon n = %3.0f points.',n))
pause(2)
end
homework
? P1.1.1
? P1.1.2
Vectorizing Function Evaluations
n=200;
x=linspace(0,1,n);
y=zeros(1,n);
for k=1:n
y(k)=((1-x(k)/24)/(1+x(k)/24+(x(k)/384)*x(k)))^8
end
plot(x,y)
x
e
xx
x
xf
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
?
?
8
2
3 8 424
1
24
1
)(
% Script File,ExpPlot
% Examines the function
% f(x) = ((1 - x/24)/(1 + x/24 + x^2/384))^8
% as an approximation to exp(-x) across [0,1].
x=linspace(0,1,200);
y=((1-x/24)./(1+x/24+(x/384).*x)).^8
plot(x,y)
title('(1-x/24)/(1+x/24+x^2/384))^8')
.^ pointwise vector exponentiation
./ pointwise vector divide
.* pointwise vector multiple
Scaling and superpositioning
x=linspace(-pi/2,11*pi/2,200);
y=tan(x);
plot(x,y)
axis([-pi/2 9*pi/2 -10 10])
-2 0 2 4 6 8 10 12 14 16 18
- 1 8
- 1 6
- 1 4
- 1 2
- 1 0
-8
-6
-4
-2
0
2
x 1 0
15
0 2 4 6 8 10 12 14
- 1 0
-8
-6
-4
-2
0
2
4
6
8
10
autoscaling
Axis([xmin xmax ymin ymax])
% Plots the function tan(x),-pi/2 <= x <= 9pi/2
close all;ymax = 10;
x = linspace(-pi/2,pi/2,40);
y = tan(x); plot(x,y)
axis([-pi/2 9*pi/2 -ymax ymax])
title('The Tangent Function');xlabel('x'); ylabel('tan(x)')
hold on
for k=1:4
xnew = x+ k*pi;
plot(xnew,y);
End
hold off 0 2 4 6 8 10 12 14-10
-8
-6
-4
-2
0
2
4
6
8
10
The Tangent Function
x
tan
(x
)
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1-1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
x=linspace(0,1,200);
y1=sin(2*pi*x); y2=cos(2*pi*x);
plot(x,y1); hold on;
plot(x,y2,'--');
plot([1/8 5/8],[1/sqrt(2) -1/sqrt(2)],'*')
hold off
x=linspace(0,1,200);
y1=sin(2*pi*x); y2=cos(2*pi*x);
plot(x,y1,x,y2,'--',[1/8 5/8],[1/sqrt(2) -1/sqrt(2)],'*')
Plot(….)
% Plots selected regular polygons.
close all;clc
theta = linspace(0,2*pi,361);
c = cos(theta);s = sin(theta);k=0;
for sides = [3 4 5 6 8 10 12 18 24]
stride = 360/sides;
k=k+1;
subplot(3,3,k)
plot(c(1:stride:361),s(1:stride:361))
axis([-1.2 1.2 -1.2 1.2])
axis equal
end
interesting
Subplot(m,n,k):Split the current
figure into an m-row by n
column array of subwindows
that are indexed left to right,top
to bottom
-1 0 1
-1
0
1
-1 0 1
-1
0
1
-1 0 1
-1
0
1
-1 0 1
-1
0
1
-1 0 1
-1
0
1
-1 0 1
-1
0
1
-1 0 1
-1
0
1
-1 0 1
-1
0
1
-1 0 1
-1
0
1
-10 -8 -6 -4 -2 0 2 4 6 8 10-20
-15
-10
-5
0
5
10
15
20
f(x)=2sin(x)+3sin(2x)+7sin(3x)+5sin(4x)
Some matrix computations
n=200;
x=linspace(-10,10,n)';
y=zeros(n,1);
for k=1:n
y(k)=2*sin(x(k))+3*sin(2*x(k))+7*sin(3*x(k))+5*sin(4*x(k));
end
plot(x,y)
title('f(x)=2sin(x)+3sin(2x)+7sin(3x)+5sin(4x)')
Vectorization
n=200;
x=linspace(-10,10,n)';
y(k)=2*sin(x(k))+3*sin(2*x(k))+7*sin(3*x(k))+5*sin(4*x(k));
plot(x,y)
title('f(x)=2sin(x)+3sin(2x)+7sin(3x)+5sin(4x)')
-10 -8 -6 -4 -2 0 2 4 6 8 10-20
-15
-10
-5
0
5
10
15
20
f(x)=2sin(x)+3sin(2x)+7sin(3x)+5sin(4x)
?
?
?
?
?
?
?
?
?
?
?
?
5
7
3
2
]4s i n,3s i n,2s i n,[ s in xxxx
)4s i n (5)3s i n (7)2s i n (3)s i n (2)( xxxxxf ????
By vector
close all
x = linspace(-10,10,200)';
A = [sin(x) sin(2*x) sin(3*x) sin(4*x)];
y = A*[2;3;7;5];
plot(x,y)
title('f(x) = 2sin(x) + 3sin(2x) + 7sin(3x) +
5sin(4x)')
for j=1:m
for k=1:n
A(k,j)=sin(j*x(k))
end
end
for j=1:m
A(:,j)=sin(j*x)
end
for k=1:n
A(k,:)=sin(k*(1:m))
end
)4s in (9)3s in (6)2s in (2)s in (8)(
)4s in (5)3s in (7)2s in (3)s in (2)(
xxxxxg
xxxxxf
????
????
n=200;x = linspace(-10,10,n)';
A = [sin(x) sin(2*x) sin(3*x) sin(4*x)];
y1 = A*[2;3;7;5]; y2=A*[8; 2; 6; 9];
plot(x,y1,x,y2)
-10 -8 -6 -4 -2 0 2 4 6 8 10-20
-15
-10
-5
0
5
10
15
20
y=A*[2 8;3 2;7 6;5 9];
plot(x,y(:,1),x,y(:,2))
plot(x,y)
By matrix
,,,
5)4(,1)3(,0)2(,2)1(
)4s in ()3s in ()2s in ()s in ()(
4321
4321
?
?????
????
????
????
ffff
xxxxxf
X=[1 2 3 4;2 4 6 8; 3 6 9 12;4 8 12 16];
Z=sin(x);
F=[-2;0;1;5];
alpha=Z\F
backslash
homework
? P1.2.1-P1.2.10
Building Exploratory
Environments
Random Process
-1 -0.5 0 0.5 1 1.5 20
10
20
30
40
50
60
Distribution of Values in rand(1000,1)
Mean = 0.502,Median = 0.510.
-3 -2 -1 0 1 2 30
10
20
30
40
Distribution of Values in randn(1000,1)
Mean = -0.043,Standard Deviation = 0.943
close all
subplot(2,1,1); x = rand(1000,1);
hist(x,30); axis([-1 2 0 60])
title('Distribution of Values in rand(1000,1)')
xlabel(sprintf('Mean = %5.3f,Median = …
%5.3f.',mean(x),median(x)))
subplot(2,1,2); x = randn(1000,1);
hist(x,linspace(-2.9,2.9,100));
title('Distribution of Values in randn(1000,1)')
xlabel(sprintf('Mean = %5.3f,Standard Deviation =…
%5.3f',mean(x),std(x)))
Clouds.m
close all; Points = rand(1000,2);
subplot(1,2,1)
plot(Points(:,1),Points(:,2),'.')
title('Uniform Distribution.')
axis([0 1 0 1])
axis square
Points = randn(1000,2);
subplot(1,2,2)
plot(Points(:,1),Points(:,2),'.')
title('Normal Distribution.')
axis([-3 3 -3 3])
axis square
2 3 4 5 6 7 8 9 10 11 12
0
20
40
60
80
100
120
140
160
180
Outcome of 1000 Dice Rolls.
% Script File,Dice
% Simulates 1000 rollings of a pair of dice.
close all
First = 1 + floor(6*rand(1000,1));
Second = 1 + floor(6*rand(1000,1));
Throws = First + Second;
hist(Throws,linspace(2,12,11));
title('Outcome of 1000 Dice Rolls.')
First = ceil(6*rand(1000,1));
-1.5 -1 -0.5 0 0.5 1 1.5
-1.5
-1
-0.5
0
0.5
1
1.5
p's estimate by Monte Carlo
Hundreds of Trials
0 50 100 150 200 250 300 350 400 450 500
3.08
3.1
3.12
3.14
3.16
3.18
3.2
3.22
3.24
3.26
Monte Carlo Estimate of Pi = 3.157
Hundreds of Trials
Darts.m
close all; rand('seed',.123456)
NumberInside = 0;
PiEstimate = zeros(500,1);
for k=1:500
x = -1+2*rand(100,1); y = -1+2*rand(100,1);
NumberInside = NumberInside + sum(x.^2 + y.^2 <= 1);
PiEstimate(k) = (NumberInside/(k*100))*4;
end
plot(PiEstimate)
title(sprintf('Monte Carlo Estimate of Pi = %5.3f',PiEstimate(500)));
xlabel('Hundreds of Trials')
Monte Carlo method
Polygon Smoothing
Points (xi,yi ); i=1:n+1
Plot(x,y,x,y,’*’)
xnew = [(x(1:n)+x(2:n+1))/2; (x(1)+x(2))/2];
ynew = [(y(1:n)+y(2:n+1))/2;(y(1)+y(2))/2];
plot(xnew,ynew)
A new polygon is displayed that is obtained by connecting
the side midpoints of the original polygon.
close all
n = input('Enter the number of edges:');figure; axis([0 1 0 1])
axis square;hold on
x = zeros(n,1); y = zeros(n,1);
for k=1:n
title(sprintf('Click in %2.0f more points.',n-k+1))
[x(k) y(k)] = ginput(1);
plot(x(1:k),y(1:k),x(1:k),y(1:k),'*')
end
x0 = x; y0 = y;x = [x;x(1)];y = [y;y(1)];plot(x,y,x,y,'*')
title('The Original Polygon') ;hold off;k=0;
xlabel('Click inside window to smooth,outside window to quit.')
[a,b] = ginput(1); v = axis;
while (v(1)<=a) & (a<=v(2)) & (v(3)<=b) & (b<=v(4));
k = k+1;
x = [(x(1:n)+x(2:n+1))/2;(x(1)+x(2))/2];
y = [(y(1:n)+y(2:n+1))/2;(y(1)+y(2))/2];
clf; plot(x,y,x,y,'*') ;axis([0 1 0 1]);
axis square
title(sprintf('Number of Smoothings = %1.0f',k))
xlabel('Click inside window to smooth,outside window to quit.')
v = axis;
[a,b] = ginput(1);
end
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
The Original Polygon
Click inside window to smooth,outside window to quit.
Smooth.m
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Number of Smoothings = 1
Click inside window to smooth,outside window to quit.
Error
? Focus on the mathematical errors that arise
through discretization and the rounding
errors that arise due to finite precision
arithmetic
Absolute and Relative Error
? Absolute error
? Relative error
? If the relative error is 10-d,then has
approximately d correct significant digits
|~| xx ?
||/|~| xxx ?
x~
% Script File,Stirling
% Prints a table showing error in Stirling's formula for n!
disp(' Stirling Absolute Relative')
disp(' n n! Approximation Error Error')
disp('----------------------------------------------------------------')
e=exp(1);nfact=1;
for n=1:13
nfact = n*nfact;
s = sqrt(2*pi*n)*((n/e)^n);
abserror = abs(nfact - s);
relerror = abserror/nfact;
s1 = sprintf(' %2.0f %10.0f %13.2f',n,nfact,s);
s2 = sprintf(' %13.2f %5.2e',abserror,relerror);
disp([s1 s2])
end
!2 n
e
nnS n
n ??
?
??
?
?? p
Stirling
Stirling Absolute Relative
n n! Approximation Error Error
1 1 0.92 0.08 7.79e-002
2 2 1.92 0.08 4.05e-002
3 6 5.84 0.16 2.73e-002
4 24 23.51 0.49 2.06e-002
5 120 118.02 1.98 1.65e-002
6 720 710.08 9.92 1.38e-002
7 5040 4980.40 59.60 1.18e-002
8 40320 39902.40 417.60 1.04e-002
9 362880 359536.87 3343.13 9.21e-003
10 3628800 3598695.62 30104.38 8.30e-003
11 39916800 39615625.05 301174.95 7.55e-003
12 479001600 475687486.47 3314113.53 6.92e-003
13 6227020800 6187239475.19 39781324.81 6.39e-003
xx
n
e
k
xe n
k
n
k
x ??
?
?? ?
?
? ?
?
0,
)!1(!0
1
% Script File,ExpTaylor Plots,as a function of n,the relative error in
%the Taylor approximation 1 + x + x^2/2! +...+ x^n/n! to exp(x).
nTerms = 50;
for x=[10 5 1 -1 -5 -10]
figure; f = exp(x)*ones(nTerms,1); s = 1; term = 1;
for k=1:50
term = x.*term/k; s = s+ term; err(k) = abs(f(k) - s);
end
relerr = err/exp(x); semilogy(1:nTerms,relerr)
ylabel('Relative Error in Partial Sum.');
xlabel('Order of Partial Sum.') title(sprintf('x = %5.2f',x))
end
0 5 10 15 20 25 30 35 40 45 50
10
-16
10
-14
10
-12
10
-10
10
-8
10
-6
10
-4
10
-2
10
0
R
e
la
t
iv
e
E
r
r
o
r
in
P
a
r
t
ia
l
S
u
m
.
O r d e r o f P a r t i a l S u m,
x = 1 0, 0 0
0 5 10 15 20 25 30 35 40 45 50
10
-10
10
-5
10
0
10
5
10
10
R
e
la
t
iv
e
E
r
r
o
r
in
P
a
r
t
ia
l
S
u
m
.
O r d e r o f P a r t i a l S u m,
x = - 1 0, 0 0
0 5 10 15 20 25 30 35 40 45 50
10
-16
10
-14
10
-12
10
-10
10
-8
10
-6
10
-4
10
-2
10
0
R
e
la
t
iv
e
E
r
r
o
r
in
P
a
r
t
ia
l
S
u
m
.
O r d e r o f P a r t i a l S u m,
x = 5, 0 0
0 2 4 6 8 10 12 14 16
10
-16
10
-14
10
-12
10
-10
10
-8
10
-6
10
-4
10
-2
10
0
R
e
la
t
iv
e
E
r
r
o
r
in
P
a
r
t
ia
l
S
u
m
.
O r d e r o f P a r t i a l S u m,
x = 1, 0 0
bottom
out
10-16
10-10
10-16 10-16
x=10
Relative error
x=-10
x=5 x=1
Rounding Errors
? The plots reveal that the mathematical
convergence theory does not quit apply
? The errors do not go to zero as the number of
terms in the series increases.
? The relative error for x=-10 is much worse that for
x=10
? Floating point arithmetic!
mathematics Development of algorithms
Inexact computer
arithmetic system
% Script File,Zoom
% Plots (x-1)^6 near x=1 with increasingly refined scale.
% Evaluation via x^6 - 6x^5 + 15x^4 - 20x^3 + 15x^2 - 6x +1
% leads to severe cancellation.
close all;k=0;n=100;
for delta = [.1,01,008,007,005,003 ]
x = linspace(1-delta,1+delta,n)';
y = x.^6 - 6*x.^5 + 15*x.^4 - 20*x.^3 + 15*x.^2 - 6*x +
ones(n,1);
k = k+1;
subplot(2,3,k); plot(x,y,x,zeros(1,n))
axis([1-delta 1+delta -max(abs(y)) max(abs(y))])
end
0.9 1 1.1
-5
0
5
x 10 -7
0.99 1 1.01
-1
-0.5
0
0.5
1 x 10
-12
0.995 1 1.005
-2
-1
0
1
2
x 10 -13
0.995 1 1.005
-1
-0.5
0
0.5
1
x 10 -13
0.995 1 1.005
-1
0
1
x 10 -14
0.998 1 1.002
-5
0
5
x 10 -15
y = x.^6 - 6*x.^5 + 15*x.^4 - 20*x.^3 + 15*x.^2 -
6*x + ones(n,1);
One
zero?
0.9 1 1.1
-1
-0.5
0
0.5
1 x 10
-6
0.99 1 1.01
-1
-0.5
0
0.5
1 x 10
-12
0.995 1 1.005
-2
-1
0
1
2
x 10 -13
0.995 1 1.005
-1
-0.5
0
0.5
1
x 10 -13
0.995 1 1.005
-1
0
1
x 10 -14
0.998 1 1.002
-5
0
5
x 10 -16
y = (x-1).^6
Algorithms that are equivalent mathematically may
behave very differently numerically.
The floating point number
x=1; p=0; y=1; z=x+y;
while x~=z
y=y/2;
p=p+1;
z=x+y;
end
disp(sprintf('p = %2.0f is the smallest
positive integer so 1+1/2^p = 1.',p))
12/11 ?? p
P=53
% q = smallest positive integer so 1/2^q = 0.
x=1; q=0;
while x>0
x=x/2;
q=q+1;
end;
02/1 ?q
q=1075
% r = smallest positive integer so 2^r = inf.
x=1; r=0;
while x~=inf
x=2*x; r=r+1;
end
in f2 ?r
r=1024
p<q?
underflow
overflow
% The Kahan example,In exact arithmetic f/e = 0/0:
disp(' ')
h = 1./2.;
x = 2./3,- h;
y = 3./5,- h;
e = (x+x+x) - h;
f = (y+y+y+y+y)-h;
z = f/e;
disp(sprintf('\nz = %10.5f',z))
z=1
Designing function
? An ability to write good Matlab function is
crucial.
? Three examples are used to clarify the
essential ideas,
– Taylor series
– Numerical differentiation
– Three-digit Arithmetic(see book)
Four ways to compute the
exponential of a vector
function y=MyExpF(x,n)
%Pre:
% x=scalar
% n=positive integer
%Post:
% y=n-th order Taylor approximation to exp(x)
y=1;
term=1;
for k=1:n
term=x*term/k;
y=y+term;
end
?
?
?
n
k
k
n k
xxT
0 !
)(
Help MyExpF
All variables inside the
function are local and
are not part of the
Matlab workspace
The input and output parameters are
formal parameters
close all;m = 50;x = linspace(-1,1,m);
exact = exp(x);figure;k=0;
y = zeros(1,m);
for n = [4 8 16 20]
for i=1:m
y(i) = MyExpF(x(i),n);
end
RelErr = abs(exact - y)./exact;
k=k+1;
subplot(2,2,k); plot(x,RelErr)
title(sprintf('n = %2.0f',n))
end
-1 -0.5 0 0.5 10
0.005
0.01
0.015
0.02
n = 4
-1 -0.5 0 0.5 10
2
4
6
8x 10
-6 n = 8
-1 -0.5 0 0.5 10
2
4
6
8x 10
-15 n = 16
-1 -0.5 0 0.5 10
2
4
6x 10
-16 n = 20
10-8
10-15 10-16
function y=MyExp1(x)
n=17;
p=length(x);
y=ones(p,1);
for i=1:p
y(i)=MyExpF(x(i),n);
end
function y=MyExp2(x)
y=ones(size(x));
n=17;
term=ones(size(x));
for k=1:n
term=x.*term/k;
y=y+term;
end
Length(x) MyExp1(x) MyExp2(x)
100 0.02 0.00
200 0.03 0.00
300 0.03 0.00
Run time
Faster than
result in book
disp(' Length(x) MyExp1(x) MyExp2(x) ')
disp(' Time Time')
disp('-----------------------------------------')
for L=50:50:300
xL=linspace(-1,1,L);
tstart=clock; y=MyExp1(xL); tfinish=clock;
t1=etime(tfinish,tstart);
tstart=clock; y=MyExp2(xL); tfinish-clock;
t2=etime(tfinish,tstart);
disp(sprintf('%6.0f %13.2f %13.2f ',L,t1,t2))
end
0 500 1000 1500 2000 2500 30000
0.05
0.1
0.15
0.2
0.25
0.3
0.35Length(x) MyExp1(x) MyExp2(x)
450 0.0470 0.0000000000000
500 0.0470 0.0160000000000
550 0.0620 0.0000000000000
600 0.0630 0.0000000000000
650 0.0620 0.0160000000000
700 0.0780 0.0000000000000
750 0.0780 0.0000000000000
800 0.0940 0.0000000000000
850 0.0940 0.0000000000000
900 0.0930 0.0000000000000
950 0.1100 0.0000000000000
function y = MyExpW(x)
% y = MyExpW(x)
% x is a scalar and y is a Taylor approximation to exp(x).
y = 0;
term = 1;
k=0;
while abs(term) > eps*abs(y)
k = k + 1;
y = y + term;
term = x*term/k;
end
function y = MyExp3(x)
% y = MyExp3(x)
% x is a column n-vector and y is a column n-vector
with the property that
% y(i) is a Taylor approximation to exp(x(i)) for i=1:n.
n = length(x);
y = ones(n,1);
for i=1:n
y(i) = MyExpW(x(i));
end
function y = MyExp4(x)
% y = MyExp4(x)
% x is an n-vector and y is an n-vector with the same shape and
%the property that y(i) is a Taylor approximation to exp(x(i)) for
i=1:n.
y = zeros(size(x));
term = ones(size(x));
k = 0;
while any(abs(term) > eps*abs(y))
y = y + term;
k = k+1;
term = x.*term/k;
end
Numerical Differentiatioin
2
2
)('')(')()( hfhafafhaf ?????
h
afhafD
h
)()( ???
2)('')('
hfafD
h ???
a=input('Enter a,');
h=logspace(-1,-16,16);
Dh=(sin(a+h)-sin(a))./h;
err=abs(Dh-cos(a));
0 2 4 6 8 10 12 14 16-20
-18
-16
-14
-12
-10
-8
-6
-4
-2
0
1.0000e-001
1.0000e-002
1.0000e-003
1.0000e-004
1.0000e-005
1.0000e-006
1.0000e-007
1.0000e-008
1.0000e-009
1.0000e-010
1.0000e-011
1.0000e-012
1.0000e-013
1.0000e-014
1.0000e-015
1.0000e-016
0.04293855333275
0.00421632485627
0.00042082550781
0.00004207444952
0.00000420736228
0.00000042074681
0.00000004182769
0.00000000296989
0.00000005254127
0.00000005848104
0.00000116870406
0.00004324021692
0.00073391590031
0.00370697619819
0.01480920644444
0.54030230586814
loss of accuracy
Explanation
|)(''|
2
|)('| ?fhafD h ??
h
fhafD h e p s2|)(''|
2
|)('| ??? ?
Truncation
error
|)(''|e p s /2 ?fh ?
minimize
The choice of h
2|)(''| Mxf ?
h
hMhe r r D ?2
2
)( 2 ??
2
2
M
h opt
?
? 22)( Mhe r r D opt ??
function d=Derivative(fname,a,delta,M2)
if nargin<=3
M2=10;
end
if nargin==2
delta=eps;
end
hopt=2*sqrt(delta/M2);
d=(feval(fname,a+hopt)-feval(fname,a))/hopt;
NARGIN returns the number of input
arguments that were used to call the function,
FEVAL(F,x1,...,xn) evaluates
the function
da=derivative(fname,a,delta,M2)
da=derivative(fname,a,delta)
da=derivative(fname,a)
da=derivative('sin',1,eps,1)
da=derivative('sin',1,eps)
da=derivative('sin',1)
