Exercise1
1.
1 ./点除
如果a、b是矩阵,a./b就是a、b中对应的每个元素相除,得到一个新的矩阵;
如果a、b是两个数,那么a./b就是普通的除法
2 /
如果a、b是矩阵, 比如aX=b,求出X。那么X=a\b(代表a的逆乘以b);或者Xa=b,X=b/a(代表b乘以a的逆)。
如果a、b是两个数,那么a/b就是普通的除法。
代码如下:
a=[3 8.5 8.5 9 9 768]
b=[5 11 14 12 16 1024]
GR0=a./b
GR=abs(GR0-0.618)
输出:
''' assignment01
a =
3.0000 8.5000 8.5000 9.0000 9.0000 768.0000
b =
5 11 14 12 16 1024
GR0 =
0.6000 0.7727 0.6071 0.7500 0.5625 0.7500
GR =
0.0180 0.1547 0.0109 0.1320 0.0555 0.1320
故8.5X14英寸的法律文书最接近黄金矩形
2.
补充:
在Matlab中,Inf为无穷大量+∞,-inf为无穷小量-∞,在Matlab程序执行时,即使遇到了以0为除数的运算,也不会终止程序的运行,而只给出一个“除0”警告,并将结果赋成inf,继续执行
NaN非数值量则是因为0/0,或者Inf/Inf型的非正常运算。另:有时候sup/inf也叫上确界/下确界
输入:
%-------02-------%
A=[1 1;1 0];
X=[1 0;0 1];
i=0;
while X<inf,
i=i+1;X=X*A;
end
i
X
输出:
i =
1476
X =
1.0e+308 *
Inf 1.3070
1.3070 0.8078
故需要执行1476此X=X*A操作
Tips:
如果变量后边加了;则不会输出双精度浮点数(double)是计算机使用的一种数据类型,使用 64 位(8字节) 来存储一个浮点数。 它可以表示十进制的15或16位有效数字,其可以表示的数字的绝对值范围大约是:-1.7E+308~1.7E+308
(E是指数的意思,比如7.823E5=782300 这里E5表示10的5次方,E代表的英文是exponent,有时也可用index number来表示。)
单精度浮点数(float) | 双精度浮点数(double) | |
---|---|---|
在内存中占有的字节数不同 | 4个字节(1byte=1B=8bit) | 8个字节 |
有效数字位数不同 | 7 | 16 |
所能表示数的范围不同 | -3.40E+38 ~ +3.40E+38 | -1.79E+308 ~ +1.79E+308 |
在程序中处理速度不同 | 更快 | 更慢 |
3.
0301
输入:
%-------0301-------%
y3=log((fibonacci(18)'));x3=[1:18];polyfit(x3,y3,1)
输出
ans =0.4834 -0.3505
当取到第200位斐波那契数时,拟合一次函数斜率为0.4812
ans =0.4812 -0.3259
0302
输入:
clc;clear;clf;
t=[];%储存运行时间
for i=20:40
t1=clock;
fibnum(i);
t2=clock;
t3=etime(t2,t1)
i
t=[t,t3];
end
figure
%取对数
ylog=log(t*1000);
%用一次函数拟合对数函数
pfitlog=polyfit([20:40],log(t*1000),1)
yfitlog=polyval(pfitlog,[20:40]);
%红*:第n个fibnum与前两项运行时间之差,蓝色圆圈:对数,青线:一次函数
plot([20:40],t,'r*',[20:40],ylog,'bo',[20:40],yfitlog)
t50=1/1000*exp(polyval(pfitlog,50))
输出:
t3 =
0.0020
i =
20
t3 =
0.0010
i =
21
t3 =
1.0000e-03
i =
22
t3 =
0.0020
i =
23
t3 =
0.0040
i =
24
t3 =
0.0060
i =
25
t3 =
0.0100
i =
26
t3 =
0.0170
i =
27
t3 =
0.0390
i =
28
t3 =
0.0440
i =
29
t3 =
0.0690
i =
30
t3 =
0.1080
i =
31
t3 =
0.1790
i =
32
t3 =
0.2670
i =
33
t3 =
0.4540
i =
34
t3 =
0.7110
i =
35
t3 =
1.1590
i =
36
t3 =
1.8540
i =
37
t3 =
3.0380
i =
38
t3 =
4.9330
i =
39
t3 =
7.9250
i =
40
pfitlog =
0.4633 -9.6310
t50 =
473.9673
图像:
可以看到,随着n的增大,函数拟合效果越好(当n<20,因运行时间过短,误差太大,故舍弃)
当n=50时,用时t50 = 473.9673
注:多次运行函数后数值pfitlog稳定在0.46左右
其他计时函数:
tic;fun;toc;t=toc-tic%matlab推荐这个 原因: in the case of a Pentium 4 with hyperthreading running Windows,there can be a significant difference between the values returnedby cputime versus tic and toc.
cputime 是CPU每个内核的运行时间,未必是程序执行时间,所以cputime一般都不准确,tic toc一般能精确测定代码运行时间,所以,一般都用tic toc
t0=cputime;fun;t=cputime-t0
t1=clock;fun;t2=clock;t=etime(t2,t1)
profile on;profile clear;fun1;fun2;fun3;profile viewer%显示各个fun运行的时间
如果运算难以量化,可以计算复杂度
disp示例:
disp([‘toc计算第’,num2str(i),‘次循环运行时间:’,num2str(toc)]);
toc计算第1次循环运行时间:2.5628
04
0401
输入:
function F = finitefern(varargin)
%FINITEFERN MATLAB implementation of the Fractal Fern.
% Michael Barnsley, "Fractals Everywhere", Academic Press, 1993.
%
% FINITEFERN with no arguments plots 100000 points.
% FINITEFERN(N) plots N points.
% FINITEFERN(N,'s') shows each step.
%
% F = FINITEFERN(N,r,c) returns an r-by-c sparse logical
% bit map array that can be viewed with
% spy(F)
% or
% image(F)
% colormap([1 1 1; 0 2/3 0])
% F can be saved in PNG (Portable Network Graphics) format with
% imwrite(full(F),'myfern.png','png','bitdepth',1)
%
% See also: FERN.
% Copyright 2014 Cleve Moler
% Copyright 2014 The MathWorks, Inc.
showstep = (nargin >= 1) && ischar(varargin{end});
if showstep || (nargout == 0)
clf
shg
set(gcf,'menubar','none','numbertitle','off','name','Finite Fern', ...
'color','black')%黑色背景
pink = [0 2/3 0];
darkred = [2/3 0 0];
pink=[242/255,156/255,177/255]
end
if showstep
finish = uicontrol('style','toggle','string','finish', ...
'value',0,'background','black');%黑色背景
end
if (nargin >= 1) && ~ischar(varargin{1})
n = varargin{1};
else
n = 100000;
end
p = [ .85 .92 .99 1.00];
A1 = [ .85 .04; -.04 .85]; b1 = [0; 1.6];
A2 = [ .20 -.26; .23 .22]; b2 = [0; 1.6];
A3 = [-.15 .28; .26 .24]; b3 = [0; .44];
A4 = [ 0 0 ; 0 .16];
x = [.5; .5];
xs = zeros(2,n);
xs(:,1) = x;
for j = 2:n
r = rand;
if r < p(1)
x = A1*x + b1;
elseif r < p(2)
x = A2*x + b2;
elseif r < p(3)
x = A3*x + b3;
else
x = A4*x;
end
xs(:,j) = x;
if showstep
h = plot(xs(1,1:n-1),xs(2,1:n-1),'.',x(1),x(2),'o');
set(h(1),'markersize',6,'color',pink);
set(h(2),'color',darkred);
axis([-3 3 0 10])
axis off
showstep = get(finish,'value') == 0;
if ~showstep, delete(finish), end
pause(.01)
end
end
if nargout == 0
plot(xs(1,:),xs(2,:),'.','markersize',1,'color',pink);
axis([-3 3 0 10])
axis off
else
if nargin < 3
r = 768; c = 1024;
else
r = varargin{2}; c = varargin{3};
end
j = round((xs(1,:)+3)/6*c);
i = round((9-.9*xs(2,:)+.5)/10*r);
F = sparse(i,j,1,r,c)~=0;
end
输出:
0402
输入:
function F = finitefern(varargin)
%FINITEFERN MATLAB implementation of the Fractal Fern.
% Michael Barnsley, "Fractals Everywhere", Academic Press, 1993.
%
% FINITEFERN with no arguments plots 100000 points.
% FINITEFERN(N) plots N points.
% FINITEFERN(N,'s') shows each step.
%
% F = FINITEFERN(N,r,c) returns an r-by-c sparse logical
% bit map array that can be viewed with
% spy(F)
% or
% image(F)
% colormap([1 1 1; 0 2/3 0])
% F can be saved in PNG (Portable Network Graphics) format with
% imwrite(full(F),'myfern.png','png','bitdepth',1)
%
% See also: FERN.
% Copyright 2014 Cleve Moler
% Copyright 2014 The MathWorks, Inc.
showstep = (nargin >= 1) && ischar(varargin{end});
if showstep || (nargout == 0)
clf
shg
set(gcf,'menubar','none','numbertitle','off','name','Finite Fern', ...
'color','black')%黑色背景
pink = [0 2/3 0];
darkred = [2/3 0 0];
pink=[242/255,156/255,177/255]
end
if showstep
finish = uicontrol('style','toggle','string','finish', ...
'value',0,'background','black');%黑色背景
end
if (nargin >= 1) && ~ischar(varargin{1})
n = varargin{1};
else
n = 100000;
end
p = [ .85 .92 .99 1.00];
A1 = [ .85 .04; -.04 .85]; b1 = [0; 1.6];
A2 = [ .20 -.26; .23 .22]; b2 = [0; 1.6];
A3 = [-.15 .28; .26 .24]; b3 = [0; .44];
A4 = [ 0 0 ; 0 .16];
x = [.5; .5];
xs = zeros(2,n);
xs(:,1) = x;
for j = 2:n
r = rand;
if r < p(1)
x = A1*x + b1;
elseif r < p(2)
x = A2*x + b2;
elseif r < p(3)
x = A3*x + b3;
else
x = A4*x;
end
xs(:,j) = x;
if showstep
h = plot(xs(2,1:n-1),xs(1,1:n-1),'.',x(2),x(1),'o');%交换x,y
set(h(1),'markersize',6,'color',pink);
set(h(2),'color',darkred);
axis([-3 3 0 10])
axis off
showstep = get(finish,'value') == 0;
if ~showstep, delete(finish), end
pause(.01)
end
end
if nargout == 0
plot(xs(2,:),xs(1,:),'.','markersize',1,'color',pink);%交换x,y
axis([10 0 3 -3])%[-3 3 0 10]
axis off
else
if nargin < 3
r = 768; c = 1024;
else
r = varargin{2}; c = varargin{3};
end
j = round((xs(1,:)+3)/6*c);
i = round((9-.9*xs(2,:)+.5)/10*r);
F = sparse(i,j,1,r,c)~=0;
end
输出:
0403
输入:
clear;clf;clc;
bg=[0 0 85];%深蓝背景
fg=[255 255 255];%白色点
sz=get(0,'screensize');
rand('state',0)
X=finitefern(500000,sz(4),sz(3));
d=fg-bg;
R=uint8(bg(1)+d(1)*X);
G=uint8(bg(2)+d(2)*X);
B=uint8(bg(3)+d(3)*X);
F=cat(3,R,G,B);
figure
imwrite(F,'myfern.png','png','bitdepth',8)
输出:
myfern.png像素点比fern多,前者为1368x912,后者为1024x768
05
输入
clc;clear;A=magic(4)
null(A);
null(A,'r')
null(sym(A))
rref(A)
输出
A =
16 2 3 13
5 11 10 8
9 7 6 12
4 14 15 1
ans =%null(A);
-0.2236
-0.6708
0.6708
0.2236
ans =%null(A,'r')
-1
-3
3
1
ans = %null(sym(A))
-1
-3
3
1
ans =%rref(A)
1 0 0 1
0 1 0 3
0 0 1 -3
0 0 0 0
说明:
rref(A)最后一行为0,说明A的列线性无关
null(A),null(sym(A)),null(A,‘r’)的结果说明Ax=0有解,即列线性无关
Exercise2
相关函数
Jocobi
:function y=Jacobi(A,b,x) D = diag(diag(A)); L = -tril(A,-1); U = -triu(A,1); B = D\(L+U); f = D\b; y = B*x+f; iter=1; while norm(y-x)>1.0e-6 x = y; y = B*x+f; iter=iter+1 y end fprintf('Number of iterations(Jacobi):%d',iter); end
Gauss-Seidel
:function y=Seidel(A,b,x) D=diag(diag(A)); L=-tril(A,-1); U=-triu(A,1); G=(D-L)\U; f=(D-L)\b; y=G*x+f; iter=1; while norm(y-x)>1.0e-6 iter=iter+1 x=y; y=G*x+f end fprintf('Number of iterations(Seidel):%d',iter); end
1.
输入:
A=[10 -1 0;-1 10 -2;0 -2 10];
b=[9 7 6]';
x=[0 0 0]';
Jacobi(A,b,x)
seidel(A,b,x)
结果:
Number of iterations(Jacobi):11
ans =
0.995789443750000
0.957894656250000
0.791578887500000
%--------------------%
Number of iterations(Seidel):7
ans =
0.995789468437500
0.957894734218750
0.791578946843750
2.
输入:
b=[9 7 6]';
x=[0 0 0]';
A1=[1 2 -2;1 1 1;2 2 1];
Jacobi(A1,b,x)
seidel(A1,b,x)
结果:
Number of iterations(Jacobi):4
ans =
-27
26
8
%-------------------%
Number of iterations(Seidel):1012
ans =
NaN
NaN
NaN
迭代矩阵谱半径比较:
B =
0 -2 2
-1 0 -1
-2 -2 0
ans =
1.233221836101528e-05
%----------------------%
G =
0 -2 2
0 2 -3
0 0 2
ans =
2
即,用Jacobi更快。原因:使用高斯赛德尔迭代矩阵G谱半径大于1,发散
3.
0301
输入与输出:
>> A=[1 2 3;4 5 6;7 8 9];b=[1 3 5]';
>> rank(A)
ans =
2
>> rref(A)
ans =
1 0 -1
0 1 2
0 0 0
因为A秩亏,所以Ax=b有无穷多解,解的集合:k*[1 -2 1]’+[2 0 1]‘ (k∈R)
0302
输入与输出:
>> A=magic(4);b=[1:4]';f=A\b
警告: 矩阵接近奇异值,或者缩放错误。结果可能不准确。RCOND = 4.625929e-18。
f =
1.0e+15 *
1.5895
4.7685
-4.7685
-1.5895
>> seidel(A,b,[0 0 0 0]')
Number of iterations(Seidel):272
ans =
Inf
Inf
Inf
-Inf
因为magic(4)是一个秩亏矩阵,Ax=b由无穷多组解,左除运算给出符合条件的一组解
高斯赛德尔结果:多次迭代后溢出
0303
输入与输出:
K>> bslashtx(A,b)
ans =
1.0e+15 *
-0.562949953421312
-1.688849860263935
1.688849860263935
0.562949953421312
解的对比:
>> norm(A*(A\b)-b)
警告: 矩阵接近奇异值,或者缩放错误。结果可能不准确。RCOND = 4.625929e-18。
ans =
10.8397
>> norm(A*bslashtx(A,b)-b)
ans =
6.3689
从舍入误差角度,使用bslashtx计算舍入误差更小,所得解更好
0304
通过lugui(A)可知:
如图,使用左除运算则第三主元为-0.8913,小于其他系数
使用bslashtx函数则经过选主元后人为减少了舍入误差
4.
输入:
J=[];
GS=[];%J,GS用来储存次数
for x=1:20
A=randn(100,100)+diag(diag(100*ones(100,100)))';%生成随机矩阵
b=ones(100,1);
x=zeros(100,1);
%声明:为方便输出 编辑Jacobi和seidel函数使其输出结果变为循环次数
J=[J,Jacobi(A,b,x)];
GS=[GS,seidel(A,b,x)];
end
x=1:20;
plot(x,J,'r.-',x,GS,'bo-')%绘制曲线
输出:
当A为2000*2000的矩阵:
输入:
J2=[];
GS2=[];%J,GS用来储存次数
for x=1:20
A2=randn(2000,2000)+diag(diag(2000*ones(2000,2000)))';%生成随机矩阵
b=ones(2000,1);
x=zeros(2000,1);
J2=[J2,Jacobi(A2,b,x)];
GS2=[GS2,seidel(A2,b,x)];
end
x=1:20;
plot(x,J2,'r.-',x,GS2,'bo-')
输出:
综上:在可迭代的情形下,高斯赛德尔迭代方法收敛速度比雅可比更快或二者速度相当
5.
输入:
a=sqrt(2)/2
A=[0 1 0 0 0 -1 0 0 0 0 0 0 0;...
0 0 1 0 0 0 0 0 0 0 0 0 0;...
a 0 0 -1 -a 0 0 0 0 0 0 0 0;...
a 0 1 0 a 0 0 0 0 0 0 0 0 ;...
0 0 0 1 0 0 0 -1 0 0 0 0 0;...
0 0 0 0 0 0 1 0 0 0 0 0 0;...
0 0 0 0 a 1 0 0 -a -1 0 0 0;...
0 0 0 0 a 0 1 0 a 0 0 0 0;...
0 0 0 0 0 0 0 0 0 1 0 0 -1;...
0 0 0 0 0 0 0 0 0 0 1 0 0;...
0 0 0 0 0 0 0 1 a 0 0 -a 0;...
0 0 0 0 0 0 0 0 a 0 1 a 0;...
0 0 0 0 0 0 0 0 0 0 0 a 1];%可以先用Excel单元格输入矩阵各分量的值再复制到matlab
b=[0 10 0 0 0 0 0 15 0 20 0 0 0 ]';
f=A\b
输出:
f =
-28.284271247461898
20.000000000000000
10.000000000000000
-29.999999999999996
14.142135623730947
20.000000000000000
0
-29.999999999999996
7.071067811865476
24.999999999999996
20.000000000000000
-35.355339059327370
24.999999999999996