数学实验练习

最新推荐文章于 2021-03-19 06:26:24 发布

陈小孩

最新推荐文章于 2021-03-19 06:26:24 发布

阅读量709

点赞数 2

分类专栏：数学实验

本文链接：https://blog.csdn.net/weixin_44226560/article/details/103130791

版权

数学实验专栏收录该内容

3 篇文章 0 订阅

订阅专栏

文章目录

Exercise1
- 1.
- 2.
- 3.
- - 0301
  - 0302
- 04
- 05
Exercise2
- `相关函数`
- 1.
- 2.
- 3.
- - 0301
  - 0302
  - 0303
  - 0304
- 4.
- 5.
- 6.

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

6.

陈小孩

关注

2
点赞
踩
1

收藏

觉得还不错? 一键收藏
0
评论
数学实验练习

Exercise11.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 7...
复制链接

扫一扫

专栏目录