Posted 2020-07-29Updated 2020-07-29笔记2 minutes read (About 303 words)

MATLAB绘制三维图形

平面网格数据的生成

利用矩阵运算生成

x = 2:6;
y = (3:8)';
X = ones(size(y))*x;
Y = y*ones(size(x));

利用meshgrid函数生成

[X, Y] = meshgrid(x, y);

其中，参数x、y为向量，存储网格点坐标的X、Y为矩阵。

x = 2:6;
y = (3:8)';
[X, Y] = meshgrid(x, y);
Z = randn(size(X));
plot3(X, Y, Z)
grid on

绘制三维曲面的函数

mesh(x, y, z, c)用于绘制三维网格图

surf(x, y, z, c)用于绘制三维曲面图

其中，x、y是网格坐标矩阵，z是网格点上的高度矩阵，c用于指定在不同高度下的曲面颜色。

mesh(z, c)

surf(z, c)

当x、y省略时，z矩阵的第2维下标当作x轴坐标，z矩阵的第1维下标当作y轴坐标。

1
2
3

t = 1:5;
z = [0.5*t; 2*t; 3*t];
mesh(z);

meshc带等高线的三维网格曲面函数

meshz带底座的三维网格曲面函数

surfc带等高线的曲面函数

surfl带光照效果的曲面函数

标准三维曲面

sphere

[x, y, z] = sphere(n)

cylinder

[x, y, z] = cylinder(R, n)

参数n决定曲面的圆滑程度，默认为20。

Posted 2020-07-26Updated 2020-07-26笔记a few seconds read (About 39 words)

数值微积分与方程求解

数值微分与数值积分

线性方程组求解

非线性方程求解与函数极值计算

常微分方程数值求解

Posted 2020-07-26Updated 2020-07-26笔记11 minutes read (About 1648 words)

数据分析与多项式计算

数据统计分析

max()：求向量或矩阵元素的最大值。

min()：求向量或矩阵元素的最小值。

当参数为向量时：

y = max(X)：返回向量 X 的最大值存入 y，如果 X 中包含复数元素，则按模取最大值。

[y, k] = max(X)：返回向量 X 的最大值存入 y，最大值元素的序号存入 k，如果 X 中包含复数元素，则按模取最大值。

当参数为矩阵时：

max(A)：返回一个行向量，向量的第 i 个元素是矩阵 A 的第 i 列上的最大值。

[Y, U] = max(A)：返回行向量 Y 和 U，Y 向量记录 A 的每列的最大值，U 向量记录每列最大值的行号。

max(A, [], dim)：dim 取 1 或 2。dim 取 1 时，该函数的功能与 max(A) 完全相同；dim 取 2 时，该函数返回一个列向量，其第 i 个元素是 A 矩阵的第 i 行上的最大值。

求矩阵的最大值：

1	>> max(max(A))

1	>> max(A(:))

mean()：求算术平均值。

median()：求中值。

sum()：求和函数。

prod()：求积函数。

cumsum()：累加和函数。

cumprod()：累乘积函数。

标准差：

std()：计算标准差函数。

调用格式：

std(X)：计算向量 X 的标准差。

std(A)：计算矩阵 A 是各列的标准差。

std(A, flag, dim)：flag 取 0 或 1，当 flag = 0 时，按 S1 公式计算样本标准差；当 flag = 1 时，按 S2 公式计算总体标准差。默认情况下，flag = 0，dim = 1。

生成满足正态分布的 50000 * 4 随机矩阵，用不同的形式求其各列之间的标准差。

x = randn(50000, 4);
y1 = std(x, 0, 1);
y2 = std(x, 1, 1);
x1 = x';
y3 = std(x1, 0, 2);
y3 = y3';
y4 = std(x1, 1, 2);
y4 = y4';

相关系数：

corrcoef()

调用格式：

corrcoef(A)：返回的相关系数矩阵的第 i 行第 j 列元素表示矩阵 A 中第 i 列和第 j 列的相关系数。

corrcoef(X, Y)：这里 X，Y是向量。

排序：

sort()

调用格式：

sort(X)：对向量 X 按升序排序。

[Y, I] = sort(A, dim, mode)：其中 dim 指明对 A 的列还是行进行排序，mode 指明升序还是降序排序，若取 “ascend”，则按升序；若取 “descend”，则按降序。默认为升序。输出参数中，Y 是排序后的矩阵，I 记录 Y 中元素在 A 中的位置。

多项式计算

多项式乘法

conv(P1, P2)：P1，P2是多项式系数向量。

多项式除法

[Q, r] = deconv(P1, P2)：Q 返回多项式 P1 除以 P2 的商式，r 返回 P1 除以 P2 的余式。

多项式的求导

polyder()

p = polyder(P)：求多项式 P 的导函数。

p = polyder(P, Q)：求 P*Q 的导函数。

[p, q] = polyder(P, Q)：求 P/Q 的导函数。导函数的分子存入 p，分母存入 q。

多项式的求值

polyval(p, x)：代数多项式求值。p 为多项式系数向量，x 可以是标量、向量或矩阵。

polyvalm(p, x)：矩阵多项式求值。x 为方阵。

多项式求根

roots(p)：p 为多项式系数向量。

若已知多项式的全部根，则可以用 poly 函数建立起该多项式：p = poly(x)。

求多项式 - 38.89 t^2 + 126.11 t - 3.42 的极大值：

p = [-38.89, 126.11, -3.42];
q = polyder(p);
a = roots(q);
b = polyval(p, a);
x = 0:0.1:2;
plot(x, polyval(p, x), a, b, 'rp')

数据插值

数据插值可以根据有限个点的取值状况，合理估算出附近其他点的取值，从而节约大量的实验和测试资源。数据插值是一种函数逼近的方法。

interp1()：一维插值函数。

Y1 = interp1(X, Y, X1, method)

根据 X，Y 的值，计算函数在 X1 处的值。其中，X，Y是两个等长的已知向量，分别表示采样点和采样值。X1 是一个向量或标量，表示要插值的点。

x = [0,3,5,7,9,11,12,13,14,15];
y = [0,1,2,1.7,2.0,2.1,2.0,1.8,1.2,1.0,1.6];
x1 = 0:0.1:15;
y1 = interp1(x, y, x1, 'spline');
plot(x1, y1)

method 用于指定插值方法，常用的取值有以下四种：

linear：线性插值，默认方法。将插值点与靠近的两个数据点用直线连接，然后在直线上选取对应插值点的数据。

nearest：最近插值点。选择最近样本点的值作为插值数据。(如果是中间点，则取后一个数据点的值)

pchip：分段 3 次埃尔米特插值。采用分段三次多项式，除满足插值条件，还需满足在若干节点处相邻段插值函数的一阶导数相等，使得曲线光滑的同时，还具有保形性。

spline：3 次样条插值。每个分段内构造一个三次多项式，使其插值函数除满足插值条件外，还要求在各节点处具有连续的一阶和二阶导数。

为什么这两种插值方法都用 3 次多项式而不用更高次的？

多项式次数并非越高越好。次数越高，越容易产生震荡而偏离原函数，这种现象称为龙格(Runge)现象。

interp2()：二维插值函数

Z1 = interp2(X, Y, Z, X1, Y1, method)

其中，X，Y 是两个向量，表示两个参数的采样点，Z 是采样点对应的函数值。X1，Y1 是两个标量或向量，表示要插值的点。

曲线拟合

x = 1790:10:2020;
y = [...];    # 此处数据省略。。。
plot(x, y, '*')
p = polyfit(x, y, 3);
polyval(p, 2020)
plot(x, y, '*', polyval(p, x))

与数据插值类似，曲线拟合也是一种函数逼近的方法，但是拟合出的曲线不一定经过所有样本点。

曲线拟合的原理：最小二乘法

polyfit()：多项式拟合函数。可求得最小二乘拟合多项式系数。

调用格式：

P = polyfit(X, Y, m)

[P, S] = polyfit(X, Y, m)

[P, S, mu] = polyfit(X, Y, m)

根据样本数据 X 和 Y，产生一个 m 次多项式 P 及其在采样点误差数据 S，mu 是一个二元向量，mu(1) 是 mean(X)，而 mu(2) 是 std(X)。

Posted 2020-07-16Updated 2020-07-19笔记9 minutes read (About 1355 words)

MATLAB绘图

二维曲线

plot函数

plot函数的基本用法

plot(x, y)

其中，x、y 分别存储 x 坐标和 y 坐标。

绘制一条折线：

1
2
3

x = [1,2,3,4];
y = [7,3,9,1];
plot(x,y)

最简单的plot函数调用格式

plot(x)

当 x 是实向量时，则以该向量元素的下标为横坐标，元素的值为纵坐标，绘制曲线。

当 x 是复向量时，则以该向量元素的实部为横坐标，虚部为纵坐标，绘制曲线。

x = [1,2,4,6];
y = [6,1,8,3];
c = complex(x, y);
plot(c)

plot(x, y)函数参数的变化形式

一般地，x、y 是长度相等的向量。

当 x 是向量，y 是矩阵时：

如果矩阵 y 的列数等于 x 的长度，则以向量 x 为横坐标，以 y 的每个行向量为纵坐标绘制曲线，曲线的条数等于 y 的行数。

如果矩阵 y 的行数等于 x 的长度，则以向量 x 为横坐标，以 y 的每个列向量为纵坐标绘制曲线，曲线的条数等于 y 的列数。

1
2
3

x = linspace(0,2*pi,100);
y = [sin(x);sin(2*x);sin(x/2)];  # y是一个矩阵
plot(x,y)

t = 0:0.01:2*pi;
t1 = t';
x = [t1,t1,t1];
y = [sin(t1),sin(2*t1),sin(0.5*t1)];
plot(x,y)

含多个输入参数的plot函数

ploy(x1, y1, x2, y2, … xn, yn)

其中，每一向量对构成一组数据的横、纵坐标，绘制一条曲线。

t1 = linspace(0,2*pi,10);
t2 = linspace(0,2*pi,20);
t3 = linspace(0,2*pi,100);
plot(t1,sin(t1),t2,sin(t2)+1,t3,sin(t3)+2)

含选项的plot函数

plot(x, y, 选项)

1
2
3

线型："-"(实线), ":"(虚线), "-."(点画线), "--"(双画线)
颜色："r"(红), "g"(绿), "b"(蓝), "w"(白), "k"(黑)
数据点标记："*"(星号), "o"(圆圈), "s"(方块), "p"(五角星), "^"(朝上三角符号)

x = (0:pi/50:2*pi)';
y1 = 2*exp(-0.5*x)*[1,-1];
y2 = 2*exp(-0.5*x).*sin(2*pi*x);
x1 = 0:0.5:6;
y3 = 2*exp(-0.5*x1).*sin(2*pi*x1);
plot(x,y1,'k:',x,y2,'b--',x1,y3,'rp')

fplot函数

fplot函数的基本用法

fplot(f, lims, 选项)

其中，f 代表一个函数，通常采用函数句柄的形式。lims 为 x 轴的取值范围，用二元向量 [xmin, xmax] 描述，默认值为 [-5, 5]。选项定义与 plot 函数相同。

1	fplot(@(x)sin(1./x),[0,0.2],'b')

双输入函数参数的用法

fplot(funx, funy, tlims, 选项)

其中，funx、funy 代表函数，通常采用函数句柄的形式。tlims 为参数函数 funx 和 funy 的自变量的取值范围，用二元向量 [tmin, tmax] 描述。

1	fplot(@(t)t.sin(t),@(t)t.cos(t),[0,10*pi],'r')

绘制图形的辅助操作

图形标注

title(图形标题)

xlabel(x 轴说明)

ylabel(y轴说明)

text(x, y, 图形说明)

legend(图例1, 图例2, …)

title函数

title函数的基本用法

x = -2*pi:0.05:2*pi;
y = sin(x);
plot(x, y);
title({'MATLAB', 'y=sin x'})   # 有多行标题时要用{}

在图形标题中使用LaTeX格式控制符

title('y=cos{\omega}t')
title('y=e^{axt}')
title('X_{1}{\geq}X_{2}')
title('{\bf y=cos{\omega}t + {\beta}}')

1
2
3

\bf 加粗
\it 斜体
\rm 正体

含属性设置的title函数

title(图形标题，属性名，属性值)

Color属性：

1	title('y=sin x','Color','r')

FontSize属性：

1	title('y=sin x','FontSize',24)

xlabel和ylabel

x = -2*pi:0.05:2*pi;
y = sin(x);
plot(x, y)
title('y = sin(x)')
xlabel('-2\pi \leq x \leq 2\pi')

text函数和gtext函数

text(x, y, 说明)

gtext(说明)

1
2
3

>> text(-2*pi, 0, '-2{\pi}')
>> text(3, 0.28, '\leftarrow sin(x)')
>> gtext('\leftarrow sin(x)')

legend函数

legend(图例1, 图例2, …)

1
2
3

x = linspace(0,2*pi,100);
plot(x,[sin(x);sin(2*x);sin(3*x)])
legend('sin(x)','sin(2x)','sin(3x)')

坐标控制

axis函数

axis([xmin, xmax, ymin, ymax, zmin, zmax])

1	axis([-pi,pi,-10,10])

axis equal：纵、横坐标轴采用等长刻度

axis square：产生正方形坐标系(默认为矩形)

axis auto：使用默认设置

axis off：取消坐标轴

axis on：显示坐标轴

给坐标系加网格和边框

网格

grid on：显示网格

grid off：不显示网格

grid：在两种状态下切换

边框

box on

box off

box

一般来说，绘图命令每执行一次，就刷新一次图形窗口，原有图形不会保持。要想保持原有图形，可以用图形保持命令：

hold on

hold off

hold

t = linspace(0,2*pi,100);
x = sin(t); y = cos(t);
plot(x, y, 'b')
hold on;
plot(2*x, 2*y, 'r--')
grid on
axis([-2.2,2.2,-2.2,2.2])
axis equal   # 显示正圆而不是椭圆

图形窗口的分割

子图：同一图形窗口中的不同坐标系下的图形称为子图。

subplot函数

subplot(m, n, p)：其中，m 和 n 指定将图形窗口分成 m * n 个绘图区，p 指定当前活动区。

x = linspace(0,2*pi,60);
subplot(2,2,1)
plot(x,sin(x)-1)
title('sin(x)-1'); axis([0,2*pi,-2,0])
subplot(2,1,2)
plot(x,cos(x)+1)
title('cos(x)+1'); axis([0,2*pi,0,2])
subplot(4,4,3)
plot(x,tan(x))
title('tan(x)');
subplot(4,4,8)
plot(x,cot(x))
title('cot(x)')

其他坐标系下的二维曲线图

对数坐标图

semilogx(x1, y1, 选项1, x2, y2, 选项2, …)

semilogy(x1, y1, 选项1, x2, y2, 选项2, …)

loglog(x1, y1, 选项1, x2, y2, 选项2, …)

x = 0:0.1:10;
y = 1./x;
subplot(2, 2, 1);
plot(x, y)
title('plot(x,y)'); subplot(2, 2, 2);
semilogx(x, y)
title('semilogx(x,y)'); grid on
subplot(2, 2, 3);
semilogy(x, y)
title('semilogy(x,y)'); grid on
subplot(2, 2, 4);
loglog(x, y)
title('loglog(x,y)'); grid on

极坐标图

polar(theta, rho, 选项)

其中，theta 为极角，rho 为极径，选项的内容与 plot 函数相同。

绘制心形线：

t = 0:pi/100:2*pi;
r = 1 - sin(t);
subplot(1, 2, 1)
polar(t, r)
subplot(1, 2, 2)
t1 = t - pi/2;
r1 = 1 - sin(t1);
polar(t, r1)

Posted 2020-07-15Updated 2020-07-16笔记2 minutes read (About 337 words)

MATLAB程序流程控制

MATLAB 程序文件包含脚本文件和函数文件。

脚本文件是可以在命令行窗口直接执行的文件，也叫命令文件。

函数文件是定义一个函数，不能直接执行，而必须以函数调用的方式来调用它。

文件的建立

用命令按钮创建文件(主页—>新建脚本)；
用 edit 命令创建文件：

1	>> edit test

顺序结构

1）数据的输入：input()

>> A = input('Please enter the value of variable A: ')
Please enter the value of variable A: 100
A =
   100

2）数据的输出：disp()

1 2	>> A = [1,2,3;4,5,6] >> disp(A)

3）程序的暂停：pause(延迟秒数)

或者，Ctrl + C。

选择结构

当条件结果为标量时，非零表示条件成立；

当条件结果为矩阵时，非空且不包含零元素表示条件成立。

c = input('Please enter a char: ','s');
if c >= 'A' && c <= 'Z'
    disp(lower(c))
elseif c >= 'a' && c <= 'z'
    disp(upper(c))
elseif c >= '0' && c <= '9'
    disp(str2double(c))
else
    disp(c)
end

switch

x = input('x = ?');
switch fix(x)
    case 2
        disp(111)
    case 2
        disp(222)
    case {3,4,5}
        disp(333)
    otherwise
        disp(444)
end

输入2、4、8时，分别输出111、333、444。

函数参数的可调性

nargin：输入实参的个数

nargout：输出实参的个数

function font=fib(a,b,c)
if nargin==1
    font=a;
elseif nargin==2
    font=a+b;
elseif nargin==3
    font=a+b+c;
end

Posted 2020-07-15Updated 2020-07-15笔记44 minutes read (About 6657 words)

MATLAB矩阵处理-续

矩阵求值

向量和矩阵的范数

1）向量的 3 种常用范数

向量1—范数：向量元素的绝对值之和。

向量2—范数：向量元素平方和的平方根。

向量∞—范数：所有向量元素绝对值中的最大值。

norm(V) 或 norm(V, 2)：计算向量 V 的 2—范数。

norm(V, 1)：计算向量 V 的 1—范数。

norm(V, inf)：计算向量 V 的 ∞—范数。

2）矩阵的范数

矩阵 A 的 1—范数：矩阵列元素绝对值之和的最大值。

矩阵 A 的 2—范数：A’A 矩阵的最大特征值的平方根。

矩阵 A 的 ∞—范数：所有矩阵行元素绝对值之和的最大值。

>> x = [2,0,1; -1,1,0; -3,3,0]
>> n = norm(x)
n =
    4.7234
>> n = norm(x, 1)
n =
     6

矩阵的条件数

矩阵 A 的条件数等于 A 的范数与 A 的逆矩阵的范数的乘积。

条件数越接近 1，矩阵的性能越好。

cond(A, 1)：计算 A 的 1—范数下的条件数。

cond(A) 或 cond(A, 2)：计算 A 的 2—范数下的条件数。

cond(A, inf)：计算 A 的 ∞—范数下的条件数。

例如，求 2~10 阶希尔伯特矩阵的条件数。

for n = 2:10
    c(n) = cond(hilb(n));
end
format long
c'

ans =

   1.0e+13 *

                   0
   0.000000000001928
   0.000000000052406
   0.000000001551374
   0.000000047660725
   0.000001495105864
   0.000047536735691
   0.001525757556663
   0.049315340455101
   1.602502816811318

矩阵的特征值与特征向量

E = eig(A)：求矩阵 A 的全部特征值，构成向量 E。

[X, D] = eig(A)：求矩阵 A 的全部特征值，构成对角阵 D，并产生矩阵 X，X 各列是相应的特征向量。

>> A = randi([1,10],3,3)
A =
     2     3     2
     8     5    10
     4     1    10
>> B = randi([1,10],2,2)
B =
     6     3
     1     4
>> A1 = zeros(3,2)
A1 =
     0     0
     0     0
     0     0
>> B1 = zeros(2,3)
B1 =
     0     0     0
     0     0     0
>> C = [A,A1;B1,B]      # 或者直接C = [A,zeros(3,2); zeros(2,3),B]
C =
     2     3     2     0     0
     8     5    10     0     0
     4     1    10     0     0
     0     0     0     6     3
     0     0     0     1     4
>> x1 = eig(A)
x1 =
   13.9462
   -1.1038
    4.1576
>> x2 = eig(B)
x2 =
     7
     3
>> x = eig(C)
x =
   13.9462
   -1.1038
    4.1576
    7.0000
    3.0000
>> [D1,E1]=eig(A)
D1 =
   -0.2887   -0.7464   -0.5276
   -0.8167    0.6308   -0.7002
   -0.4996    0.2121    0.4811
E1 =
   13.9462         0         0
         0   -1.1038         0
         0         0    4.1576
>> [D2,E2]=eig(B)
D2 =
    0.9487   -0.7071
    0.3162    0.7071
E2 =
     7     0
     0     3
>> [D,E]=eig(C)
D =
   -0.2887   -0.7464   -0.5276         0         0
   -0.8167    0.6308   -0.7002         0         0
   -0.4996    0.2121    0.4811         0         0
         0         0         0    0.9487   -0.7071
         0         0         0    0.3162    0.7071
E =
   13.9462         0         0         0         0
         0   -1.1038         0         0         0
         0         0    4.1576         0         0
         0         0         0    7.0000         0
         0         0         0         0    3.0000

MATLAB 提供了一个 eigshow 函数，可以演示单位圆上的向量 x 和 Ax 之间的关系。

稀疏矩阵

稀疏存储方式只存储矩阵的非零元素的值及其位置，即行号和列号。

A = sparse(S)：将矩阵 S 转化为稀疏存储方式的矩阵 A。

S = full(A)：将矩阵 A 转化为完全存储方式的矩阵 S。

>> A = sparse(eye(5))
A =
   (1,1)        1
   (2,2)        1
   (3,3)        1
   (4,4)        1
   (5,5)        1
>> B = full(A)
B =
     1     0     0     0     0
     0     1     0     0     0
     0     0     1     0     0
     0     0     0     1     0
     0     0     0     0     1
>> whos
  Name      Size            Bytes  Class     Attributes

  A         5x5               128  double    sparse    
  B         5x5               200  double

直接建立稀疏存储矩阵

sparse 函数的其他调用格式：

sparse(m, n)：生成一个 m*n 的所有元素都是零的稀疏矩阵。

sparse(u, v, S)：其中 u, v, S 是 3 个等长的向量。S 是要建立的稀疏存储矩阵的非零元素，u(i)、v(i) 分别是 S(i) 的行和列下标。

用 spconvert 函数直接建立稀疏存储矩阵，其调用格式为：B = spconvert(A)

A 为一个 m 3 或 m 4 的矩阵，其每行表示一个非零元素，m 是非零元素的个数。

A (i, 1) 表示第 i 个非零元素所在的行；

A (i, 2) 表示第 i 个非零元素所在的列；

A (i, 3) 表示第 i 个非零元素所在的实部；

A (i, 4) 表示第 i 个非零元素所在的虚部；

若矩阵的全部元素都是实数，则无须第 4 列。

>> A = [2,2,1;2,1,-1;2,4,3]
A =
     2     2     1
     2     1    -1
     2     4     3
>> B = spconvert(A)
B =
   (2,1)       -1
   (2,2)        1
   (2,4)        3

带状稀疏矩阵的稀疏存储

稀疏矩阵有两种基本类型：无规则结构的稀疏矩阵与有规则结构的稀疏矩阵。

带状稀疏矩阵是指所有非零元素集中在对角线上的矩阵。

[B, d] = spdiags(A)：从带状稀疏矩阵 A 中提取全部非零对角线元素赋给矩阵 B 及其位置向量 d。

A = spdiags(B, d, m, n)：产生带状稀疏矩阵的稀疏存储矩阵 A。其中 m，n 为原带状稀疏矩阵的行数和列数，矩阵 B 的第 i 列即为原带状稀疏矩阵的第 i 条非零对角线，向量 d 为原带状稀疏矩阵的所有非零对角线的位置。

>> A = [11,0,0,12,0,0;0,21,0,0,22,0;0,0,31,0,0,32;41,0,0,42,0,0;0,51,0,0,52,0]
A =
    11     0     0    12     0     0
     0    21     0     0    22     0
     0     0    31     0     0    32
    41     0     0    42     0     0
     0    51     0     0    52     0
>> [B, d] = spdiags(A)
B =
     0    11    12
     0    21    22
     0    31    32
    41    42     0
    51    52     0
d =
    -3
     0
     3
>> A = spdiags(B, d, 5, 6)
A =
   (1,1)       11
   (4,1)       41
   (2,2)       21
   (5,2)       51
   (3,3)       31
   (1,4)       12
   (4,4)       42
   (2,5)       22
   (5,5)       52
   (3,6)       32

单位矩阵的稀疏存储

speye(m, n) 返回一个 m * n 的稀疏存储单位矩阵。

>> speye(3)
ans =
   (1,1)        1
   (2,2)        1
   (3,3)        1

avatar

kf1 = [1,1,2,1,0];
k0 = [2,4,6,6,1];
k1 = [0,3,1,4,2];
B = [kf1,k0,k1];
d = [-1;0;1];
A = spdiags(B, d, 5, 5);
f = [0,3,2,1,5];
x = A \ f

注意：当参与运算的数据对象不全是稀疏存储矩阵时，所得结果是完全存储形式。

Posted 2020-07-15Updated 2020-07-15笔记3 hours read (About 27671 words)

MATLAB矩阵处理

特殊矩阵

zeros函数：产生全 0 矩阵

ones函数：产生全 1 矩阵

eye函数：产生对角线为 1 的矩阵。当矩阵是方阵时，得到一个单位矩阵。

rand函数：产生（0，1）区间均匀分布的随机矩阵。

randn函数：产生均值为 0，方差为 1 的标准正态分布随机矩阵。

randi函数：生成均匀分布的伪随机整数。

randi (iMax) 在开区间 (0, iMax) 生成均匀分布的伪随机整数；

randi (iMax, m, n) 在开区间 (0, iMax) 生成 m * n 型随机矩阵；

randi ([iMin, iMax], m, n) 在开区间 (iMin, iMax) 生成 m * n 型随机矩阵。

详细可参考：click here

1
2
3

zeros(m):产生 m * m 零矩阵
zeros(m, n):产生 m * n 零矩阵
zeros(size(A)): 产生与矩阵A同样大小的零矩阵

>> A = zeros(2, 3)
A =
     0     0     0
     0     0     0
>> A = zeros(size(reshape(A, 3, 2)))
A =
     0     0
     0     0
     0     0

通用的特殊矩阵

eg. 首先产生 5 阶两位随机整数矩阵 A，再产生均值为 0.6、方差为 0.1 的 5 阶正态分布随机矩阵 B，最后验证 (A + B) I = I A + B I，其中 I 为单位矩阵。

rand函数：产生（0，1）区间均匀分布的随机数 x。

fix(a + (b - a + 1) * x)：产生 [a, b] 区间上均匀分布的随机整数。

randn函数：产生均值为 0、方差为 1 的标准正态分布随机数 x。

μ + σx：得到均值为 μ、方差为 σ2 的随机数。

>> A = fix(10 + (99 - 10 + 1) * rand(5))    # 产生5行5列随机矩阵A
>> B = 0.6 + sqrt(0.1) * randn(5)           # 均值为0.6，方差为0.1的5行5列正态分布随机矩阵B
>> C = eye(5)
>> (A + B) * C == C * A + B * C

用于专门学科的特殊矩阵

魔方矩阵—Magic Square

是有相同的行数和列数，并在每行每列、对角线上的和都相等的矩阵。魔方矩阵中的每个元素不能相同。

>> M = magic(3)
M =
     8     1     6
     3     5     7
     4     9     2

n 阶魔方阵由 1，2，3，…，n2 共 n2 个整数组成，且每行、每列以及主、副对角线上各 n 个元素之和都相等。

n 阶魔方阵每行每列元素之和为(1 + 2 + 3 + … + n2) / n = (n + n3) / 2。

注：1 + 2 + 3 +…+ n2 = (n2 + n4) / 2，可用数学归纳法证明。

magic(n) 产生的是一个特定的魔方阵。

>> M = magic(8)
>> sum(M(1,:))
ans =
   260
>> sum(M(:,:))
ans =
   260   260   260   260   260   260   260   260
>> sum(M(:,1))
ans =
   260

范德蒙矩阵

vander(V) 生成以向量 V 为基础的 Vandermonde 矩阵。

>> A = vander(1:5)
A =
     1     1     1     1     1
    16     8     4     2     1
    81    27     9     3     1
   256    64    16     4     1
   625   125    25     5     1

Vandermonde 矩阵常用在各种通信系统的纠错编码中，例如，常用的 Reed-Solomon 编码就以其为基础。

希尔伯特矩阵

H (i, j) = 1 / (i + j - 1)，i、j分别为其行标和列标。

hilb(n) 生成n阶 Hilbert 矩阵。

invhilb(n) 求n阶 Hilbert 矩阵的逆，使用其他方法求逆会因为原始数据的微小扰动而产生不可靠的计算结果。

>> format rat       # 以有理数形式输出
>> H = hilb(4)
H =
       1              1/2            1/3            1/4     
       1/2            1/3            1/4            1/5     
       1/3            1/4            1/5            1/6     
       1/4            1/5            1/6            1/7

伴随矩阵

compan(p)，其中 p 是一个多项式的系数向量，高次幂系数在前，低次幂系数在后。

# x^3-2x^2-5x+6
>> p = [1, -2, -5, 6]
>> A = compan(p)
A =
     2     5    -6
     1     0     0
     0     1     0

帕斯卡矩阵

帕斯卡矩阵的第一行元素和第一列元素都为 1，其余位置的元素是该元素的左边元素与上面元素相加。

1 2	P(i, j) = P(i - 1, j) + P(i, j - 1) P(1, j) = P(i, 1) = 1

pascal(n)

>> format rat
>> P = pascal(5)
P =
       1              1              1              1              1
       1              2              3              4              5
       1              3              6             10             15
       1              4             10             20             35
       1              5             15             35             70
>> inv(P)       # 其逆矩阵所有元素都是整数
ans =

       5            -10             10             -5              1
     -10             30            -35             19             -4
      10            -35             46            -27              6
      -5             19            -27             17             -4
       1             -4              6             -4              1

矩阵变换

对角矩阵：只有对角线上有非零元素的矩阵。

数量矩阵：对角线上的元素相等的对角矩阵。

单位矩阵：对角线上的元素都为 1 的对角矩阵。

提取对角线元素

diag(A)：提取矩阵 A 的主对角线元素，产生一个列向量。

diag(A)：提取矩阵 A 第 k 条对角线的元素，产生一个列向量。

构造对角矩阵

diag(V)：以向量 V 为主对角线元素，产生对角矩阵。

diag(V, k)：以向量 V 为第 k 条对角线元素，产生对角矩阵。

上三角矩阵

triu(A)：提取矩阵 A 的主对角线及以上的元素。

triu(A, k)：提取矩阵 A 的第 k 条对角线及以上的元素。

>> triu(ones(4), -1)
ans =
     1     1     1     1
     1     1     1     1
     0     1     1     1
     0     0     1     1

下三角矩阵

tril(A)

矩阵的转置

转置运算符是小数点后面接单引号( .’ )。

共轭转置，其运算符是单引号( ‘ )，它在转置的基础上还要取每个数的复共轭。

>> A = [1, 3; 3+4i, 1-2i]
A =
   1.0000 + 0.0000i   3.0000 + 0.0000i
   3.0000 + 4.0000i   1.0000 - 2.0000i
>> A.'
ans =
   1.0000 + 0.0000i   3.0000 + 4.0000i
   3.0000 + 0.0000i   1.0000 - 2.0000i
>> A'
ans =
   1.0000 + 0.0000i   3.0000 - 4.0000i
   3.0000 + 0.0000i   1.0000 + 2.0000i

矩阵的翻转

fliplr(A)：对矩阵 A 实施左右翻转。

flipud(A)：对矩阵 A 实施上下翻转。

例如，验证魔方阵的主对角线、副对角线元素之和相等。

>> A = magic(5)
>> a1 = sum(diag(A))
>> B = flipud(A)
>> a2 = sum(diag(B))
>> a1 == a2

矩阵求值

方阵的行列式

det(A)：求方阵 A 所对应的行列式的值。

验证 det(A^-1) = 1 / det(A).

>> format rat
>> A = randi([1,10], 3, 3)
A =
       5             10              4
       1              5             10
      10              5              4
>> B = inv(A)
B =
      -1/21          -2/63           8/63
      16/105         -2/63         -23/315
      -1/14           5/42           1/42
>> a2 = det(B)
a2 =
       1/630
>> a1 = 1/det(A)
a1 =
       1/630

矩阵的秩

rank(A)

求 3~20 阶魔方阵的秩。

for n = 3:20
    r(n) = rank(magic(n));
end
bar(r);
grid on
axis([2,21,0,20])

avatar

矩阵的迹

等于矩阵的对角线元素之和，也等于特征值之和。

trace(A)

Posted 2020-07-14Updated 2020-07-15笔记18 minutes read (About 2737 words)

MATLAB学习记录

注：这篇文章主要是我学习基础语法时的一些零散记录。

MATLAB数值数据

数值数据类型

整型

无/有符号 8/16/32/64 位整数

数据类型转换函数

>> x = uint8(129)
129
>> x = int8(129)    # 因为int最大是127，129>127，所以返回最大值
127
>> x = int8(-129)   # int最小是-128
-128

浮点型

单精度：single 函数

双精度：double 函数

数值数据默认是双精度型。

>> class(4)
ans =
     'double'
>> class(single(4))
ans =
     'single'

复型

实部和虚部默认双精度型。虚数单位用 i 或 j 来表示。

1	6+5i 和 6+5j 表示的是一样的

real 函数求实部，imag 函数求虚部。

数值数据的输出格式

format 命令：>> format long（转换为 long 型输出）

转换回默认类型 short：>> format

format 命令只影响数据输出格式。

常用数学函数

abs 函数可以求实数的绝对值、复数的模、字符串的ASCII码值。

B = rem(A, 10) 是对 A 的每一个元素模10取余，得到一个布尔矩阵。

power(x, 2) 等价于 x^2。

取整函数：

fix：向0取整

floor：向下取整

ceil：向上取整

round：四舍五入

1 2	>> x = sin(pi/2) # 弧度制 >> x = sind(90) # 角度制

1 2	>> x = sqrt(7)-2i # 复数 >> y = exp(pi/2) # e的pi/2次幂

>> x = 1 : 100
>> k = isprime(x)   # 输出的是0/1的布尔值
>> k1 = find(k)     # 输出的是1:100之间的素数
>> p = x(k)         # 同上
>> p1 = x(k1)       # 同上

变量

变量名以字母开头，只能由字母、数字和下划线组成，最多63个字符。变量名区分字母的大小写。标准函数名以及命令名必须用小写字母。

预定义变量

ans 是默认赋值变量

i 和 j 代表虚数单位

pi 代表圆周率

NaN 代表非数

变量管理

用 who 和 whos 查看变量。

用于保存MATLAB工作区变量的文件叫内存变量文件，扩展名为.mat，也叫 MAT 文件。

save 命令：创建内存变量文件

load 命令：装入内存变量文件

1 2	>> save mydata a x # mydata.mat >> load mydata

MATLAB矩阵

矩阵的建立

1）直接输入法

将矩阵元素用 [ ] 括起来，按矩阵行的顺序输入各元素，同一行元素用逗号或空格分隔，不同行用分号分隔。

>> A = [1, 2, 3; 4, 5, 6]
A =
     1     2     3
     4     5     6

2）利用已建好的矩阵建立更大的矩阵

1
2
3

>> A = [1,2,3;4,5,6;7,8,9]
>> B = [-1,0,3;4,2,-5;9,-6,1]
>> C = [A, B; B, A]

还可以用实部矩阵和虚部矩阵构成复数矩阵

1
2
3

>> B = [1,2,3;4,5,6]
>> C = [6,7,8;9,10,11]
>> A = B + i*C

冒号表达式

1	e1 : e2 : e3

e1 为初始值，e2 为步长，e3 为终止值。

1
2
3

>> t = 0 : 1 : 5
t =
     0     1     2     3     4     5

e2 可以省略，默认为 1.

linspace函数

1	linspace(a, b, n)

a 为第一个元素，b 为最后一个元素，n 为元素总数。

1
2
3

>> x = linspace(0, pi, 6)
x =
     0    0.6283    1.2566    1.8850    2.5133    3.1416

结构矩阵和单元矩阵

结构矩阵

格式为：结构矩阵元素 . 成员名 = 表达式

1 2	>> a(1).x1 = 10; a(1).x2 = 'A'; a(1).x3 = [1, 2, 3]; >> a(2).x1 = 16; a(2).x2 = 'K'; a(2).x3 = [9, 2, 7];

单元矩阵

和一般矩阵相似，但是要用 { } 括起来。

>> b = {10,'A',[1,2,3]; 16,'K',[9,2,7]}   # b是 2×3 cell数组
b =
     {[10]}    {'A'}    {1×3 double}
     {[16]}    {'K'}    {1×3 double}

矩阵元素的引用

通过下标引用

1	>> A(4, 5) = 10 # 矩阵A的第4行第5列元素赋值为10

>> A = [1,2,3;4,5,6]   # A是2行3列
>> A(4, 5) = 10        # 超出矩阵范围，自动填充为0
A =
     1     2     3     0     0
     4     5     6     0     0
     0     0     0     0     0
     0     0     0     0    10

通过序号引用

在MATLAB中，矩阵元素按列存储。矩阵元素的序号就是其在内存中的排列顺序。

>> A = [1,2,3;4,5,6]
A =
     1     2     3
     4     5     6
>> A(3)
ans =
     2

序号与下标是一一对应的，m n 矩阵 A 的元素 A( i, j ) 的序号为 ( j - 1 ) m + i 。

矩阵元素的序号与下标可以通过 sub2ind 和 ind2sub 函数实现相互转换。

sub2ind

1	D = sub2ind(S, I, J) 下标转序号

S 是行数和列数组成的向量，通常用 size(A) 获取。I 是转换矩阵元素的行下标，J 是列下标。如果 I 和 J 是矩阵，则表示要将多个元素的行列下标转换为序号。

>> A = [1:3;4:6]
A =
     1     2     3
     4     5     6
>> B = sub2ind(size(A),[1,2;2,2],[1,1;3,2])
B =
     1     2
     6     4

size(A) 得到的是向量 [2, 3] 。然后依次得到：

行下标为 1，列下标为 1 的元素的序号是 1；

行下标为 2，列下标为 1 的元素的序号是 2；

行下标为 2，列下标为 3 的元素的序号是 6；

行下标为 2，列下标为 2 的元素的序号是 4。

ind2sub

1	[I, J] = ind2sub(S, D) 序号转下标

D 是序号。

>> [I, J] = ind2sub([3,3], [1,3,5])
I =
     1     3     2
J =
     1     1     2

利用冒号表达式获得子矩阵

A(i, :)              第i行的全部元素
A(:, j)              第j列的全部元素
A(i:i+m, k:k+n)      第i~i+m行且在第k~k+n列的全部元素
A(i:i+m, :)          第i~i+m行的全部元素

end运算符：表示某一维的末尾元素下标。（注意是某一维，可以是行或列）

1 2	>> A(end, :) # 最后一行 >> A([1,4], 3:end) # 第1、4行从第3列到最后一列的元素

删除矩阵中的元素

利用空矩阵删除矩阵的元素

空矩阵是指没有任何元素的矩阵。

>> x = []

>> A = [1,2,3,0,0;7,0,9,2,6;1,4,-1,1,8]
>> A(:, [2,4]) = []          # 删除了第2、4列元素
A =
     1     3     0
     7     9     6
     1    -1     8

改变矩阵的形状

reshape(A, m, n)：在矩阵总元素保持不变的前提下，将矩阵 A 重新排成 m * n 的二维矩阵。

reshape 函数只改变原矩阵的行数和列数，但并不改变原矩阵元素个数及其存储顺序。

A(:) 是将 A 的每一列元素堆叠起来，成为一个列向量。

>> A = [1,2,3;4,5,6]
>> B = A(:)
B =
     1
     4
     2
     5
     3
     6

所以这里 A(:) 等价于 reshape(A, 6, 1)。

MATLAB基本运算

算术运算

非奇异矩阵是行列式不为 0 的矩阵，也就是可逆矩阵。n 阶方阵 A 是非奇异方阵的充要条件是 A 为可逆矩阵，即A的行列式不为零。矩阵(方阵)A可逆与矩阵A非奇异是等价的概念。

除法运算：

在MATLAB中，有两种矩阵除法运算：右除/和左除\。

如果 A 矩阵是非奇异方阵，则 B/A等效于B inv(A)，A\B 等效于 inv(A) B。其中 inv(A) 是 A 的逆矩阵。

乘方运算：^

点运算：.*、./、. 和 .^。

两矩阵进行点运算是指它们的对应元素进行运算，所以两矩阵必须同型。

A . B 和 A B 是不一样的。前者是 A 中每个元素与 B 中对应位置的元素相乘，得到与 A 、B 都同型的矩阵。而 A * B 就是普通的矩阵乘法。

例如，当 x = 1, 4, 7 时，分别求 y = sinxcosx 的值：

1 2	>> x = 1:3:7 >> y = sin(x) .* cos(x)

关系运算

关系运算符：

1	< <= > >= == ~=

当参与比较的是两个同型的矩阵时，是对应位置的元素按照标量关系运算规则逐个进行，最后得到的结果是一个与原矩阵同型的矩阵，它的元素由 0 / 1 组成。

当参与比较的是一个标量和一个矩阵时，是这个标量与矩阵中的每一个元素进行比较，最后得到的结果是一个与原矩阵同型的矩阵，它的元素由 0 / 1 组成。

逻辑运算符：

1	&(与) \|(或) ~(非)

在算术运算、关系运算与逻辑运算中，算术运算的优先级最高，逻辑运算的优先级最低，但非(~)是单目运算符，优先级高于双目运算符。

例如判断水仙花数（一个三位数，其三位上的数字的立方之和等于这个数本身）：

m = 100 : 999;
m1 = rem(m, 10);
m2 = rem(fix(m/10), 10);
m3 = fix(m/100);
k = find(m == m1.*m1.*m1 + m2.*m2.*m2 + m3.*m3.*m3);
s = m(k);

>> k      # 序号
k =
    54   271   272   308
>> s      # 元素
s =
   153   370   371   407

字符串处理

在MATLAB中，字符串是用单引号括起来的字符序列。

>> s = 'preccrep workstation'
>> ss = s(1:8)
ss =
    'preccrep'

若字符串中的字符含有单引号，则该字符要用两个单引号来表示。

1
2
3

>> 'I''m a hacker.'
ans =
    'I'm a hacker.'

还可建立多行字符串，形成字符串矩阵：

>> ch = ['abcdef'; '123456']
>> ch(2, 3)
ans =
    '3'

例如：建立一个字符串向量，然后做如下处理：

取第1~5个字符作为子字符串；
将字符串倒序排列；
将字符串中的小写字母全部变为大写；
统计字符串中的小写字母数。

>> ch = 'ABc123dhUn'
>> subch = ch(1:5)
>> revch = ch(end:-1:1)
>> k = find(ch>='a' & ch<='z')
>> ch(k) = ch(k)-('a'-'A')
>> length(k)

字符串的执行

eval(s) 是将字符串 s 当作命令执行，也就是用来执行一个字符串表达式，并返回表达式的值。

1
2
3

t = pi;
m = '[t, sin(t), cos(t)]';
r = eval(m);

1
2
3

>> r
r =
    3.1416    0.0000   -1.0000

字符串与数值之间的转换

abs 和 double 都可以用来获取字符串矩阵所对应的 ASCII 码数值矩阵。

char 函数可以把 ACSII 码矩阵转换为字符串矩阵。

>> s1 = 'MATLAB'
>> a = ans(s1)
a =
    77    65    84    76    65    66
>> s2 = char(a+32)
s2 =
    'matlab'

字符串的比较

两种方法：利用关系运算符比较和字符串比较函数。

关系运算符比较：两个字符串长度相等时，每个字符逐个比较，返回结果是一个数值向量，其元素为0/1。

字符串比较函数：用来判断字符串是否相等。

strcmp() 函数的基本功能是比较两个字符串是否相等，相等返回 1，否则返回 0。其基本用法是：

1	ans = strcmp(s1, s2)

但是，如果我们要查找字符串数组中等于某字符串的索引时，该如果操作？strcmp() 函数也提供了这个功能，用法相同，其中，s1是字符串数组，s2是字符串，返回值为逻辑类型，大小与字符串数组s1相同。

strncmp(s1, s2, n) 比较前 n 个字符。
strcmpi(s1, s2) 在忽略字母大小写的情况下，比较两字符串是否相等。
strncmpi(s1, s2, n) 同理。

字符串的查找与替换

findstr(s1, s2) 返回短字符串在长字符串中的开始位置。

strrep(s1, s2, s3) 将字符串 s1 中的所有子串 s2 替换为字符串 s3。

strfind() 用于查找字符串数组中包含某字符串的记录：

1	k = strfind(str, pattern)

输出结果k表示pattern在str中出现的位置，若不出现，则返回[]。比如：

S = ‘Find the starting indices of the pattern string’;
k = strfind(S, ‘in’)
k = 
     2     15     19     45

平面网格数据的生成

利用矩阵运算生成

利用meshgrid函数生成

绘制三维曲面的函数

标准三维曲面

数值微分与数值积分

线性方程组求解

非线性方程求解与函数极值计算

常微分方程数值求解

数据统计分析

多项式计算

数据插值

曲线拟合

二维曲线

plot函数

plot函数的基本用法

最简单的plot函数调用格式

plot(x, y)函数参数的变化形式

含多个输入参数的plot函数

含选项的plot函数

fplot函数

fplot函数的基本用法

双输入函数参数的用法

绘制图形的辅助操作

图形标注

title函数

xlabel和ylabel

text函数和gtext函数

legend函数

坐标控制

axis函数

给坐标系加网格和边框

图形窗口的分割

其他坐标系下的二维曲线图

对数坐标图

极坐标图

矩阵求值

矩阵的特征值与特征向量

稀疏矩阵

特殊矩阵

矩阵变换

矩阵求值

MATLAB数值数据

数值数据类型

整型

浮点型

复型

数值数据的输出格式

常用数学函数

变量

预定义变量

变量管理

MATLAB矩阵

矩阵的建立

冒号表达式

结构矩阵和单元矩阵

结构矩阵

单元矩阵

矩阵元素的引用

通过下标引用

通过序号引用

sub2ind

ind2sub

利用冒号表达式获得子矩阵

删除矩阵中的元素

利用空矩阵删除矩阵的元素

改变矩阵的形状

MATLAB基本运算

算术运算

关系运算

字符串处理

Links

Categories

Recents

Archives

Tags

Subscribe for updates