matlab基础语法

任意符号矩阵生成

利用命令

1
A = sym('a',[n,m]) , B = sym('b%d%d',[n,m])

其中，$n,m$表示矩阵的行数和列数，若$n,m$其中一个值为1，那么我们定义的是一个向量（行/列）。

若你想使用方阵，但无需给出具体的值，上述两个命令都可以实现，不同的是，前者的形式为ai_j，后者为bij。注意的是，单下标的列向量无需指定下标形式。

例子：

1
A = sym('a%d%d',[3,4]), f = sym('f',[3,1])

结果：

1
2
3
4
5
6
7
8
9
10
11
12
13
>> A
A =
 
[a11, a12, a13, a14]
[a21, a22, a23, a24]
[a31, a32, a33, a34]

>>f
f =
 
f1
f2
f3

其他数据结构

讲一些常用的例子

基本赋值结构

给定函数$f(x)$，我们想计算一些嵌套函数，比如$f(f(f(x)))$，那么我们可以用如下的方法解决：

例子：若$f(x) = x^2-x-1$，试着求$F(x) = f(f(f(f(f(f(x))))))$，若结果是多项式，则$\deg f$是多少？

使用如下命令，可以在MATLAB得到结果

1
2
3
4
syms x
f(x) = x^2-x-1
F(x) = f(f(f(f(f(x)))))
F1 = expand(F)

前者F(x)会得到一个没有展开的多项式，后者expand命令则表示将F全部展开。

1
2
3
4
5
>> F(x) = f(f(f(f(f(x)))))
 
F(x) =
 
(x + (x + (- x^2 + x + 1)^2 - x^2 - (x + (- x^2 + x + 1)^2 - x^2)^2 + 1)^2 + (- x^2 + x + 1)^2 - x^2 - (x + (- x^2 + x + 1)^2 - x^2)^2)^2 - (x + (- x^2 + x + 1)^2 - x^2 - (x + (- x^2 + x + 1)^2 - x^2)^2 + 1)^2 - x - (- x^2 + x + 1)^2 + x^2 + (x + (- x^2 + x + 1)^2 - x^2)^2 - 1
1
2
3
4
5
>> F1 = expand(F)
F1(x) =
 
x^32 - 16*x^31 + 96*x^30 - 200*x^29 - 444*x^28 + 2968*x^27 - 3052*x^26 - 11804*x^25 + 30944*x^24 + 9832*x^23 - 112076*x^22 + 68760*x^21 + 216048*x^20 - 279328*x^19 - 218500*x^18 + 525046*x^17 + 47653*x^16 - 591744*x^15 + 164520*x^14 + 416864*x^13 - 232442*x^12 - 175256*x^11 + 154806*x^10 + 35134*x^9 - 58854*x^8 + 1496*x^7 + 12614*x^6 - 2044*x^5 - 1382*x^4 + 300*x^3 + 69*x^2 - 9*x - 1

我们可以简单的看出，在expand(F)后，最高次是32。

函数调用

函数调用的基本结构为

1
[返回元列表] = fun_name(输入变元列表)

其中fun_name指的是调用的函数名字，一般函数名对应于一个MATLAB路径下的文件，比如：函数名my_fun应该对应于my_fun.m文件，一些函数名需要对应于MATLAB内核中的内核函数，比如inv()函数之类。

返回变元列表，输入变元列表指的是由若干个变量名组成，变量名之间用逗号分隔，或者空格分隔。例如 [U S V] = svd(X)，该函数对给定的X矩阵进行奇异值分解，所得结果由 U、S、V三个变量返回。

当然，格式不是死的，我们可以根据需要进行调用。譬如内核函数d = eig(A)，则直接返回A矩阵的特征向量d，若调用格式[V,D] = eig(A)，则除了返回特征值D之外，还返回特征向量矩阵V，若调用格式为eig(A,B) 则将求解广义特征值问题。

冒号表达式

现在，我们想生成一些数据，可以使用冒号表达式，它的格式为v = s1:s2:s3，其中s1是起点，s2是步距，s3是终点。该函数v生成一个从s1开始，隔s2取一个点，直到不超过s3的最大值的元素组成的向量。若省略s2，则默认s2 = 1.

子矩阵提取

提取子矩阵的具体方法为：B = A(v1,v2)，其中v1表示保留的行号构成的向量，而v2表示要保留的列号构成的向量。

若想删除矩阵的第i行元素，则直接给出简单的命令A(i,:) = []，i也可以是向量。其中:在v1表示提取全部行，v2表示提取全部列。

我们用如下的命令定义一个矩阵A

1
2
A = [1,2,3;4,5,6;7,8,9]
B1 = A(1:2:end, :)

第一行建立一个3x3矩阵，第二行提取奇数行，:表示带有全部列

结果：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
>> A = [1:3;4:6;7:9]

A =

     1     2     3
     4     5     6
     7     8     9

>> B1 = A(1:2:end,:)

B1 =

     1     2     3
     7     8     9

B2 = A([3,2,1],[1,1,1])将提取A矩阵的3,2,1行，并且由3个首列构成矩阵。改成[1,3,1]则会吧第二列变成A的第三列倒转。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
>> B2 = A([3,2,1],[1,1,1])

B2 =

     7     7     7
     4     4     4
     1     1     1

>> B2 = A([3,2,1],[1,3,1])

B2 =

     7     9     7
     4     6     4
     1     3     1

使用命令B3 = A(:,end:-1:1)可以让A矩阵左右翻转。

1
2
3
4
5
6
7
>> B3 = A(:,end:-1:1)

B3 =

     3     2     1
     6     5     4
     9     8     7

A(2,:) = []删除第二行, A(:,3) = []删除第三列

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

A =

     1     2
     7     8

接下来我们讲讲一些函数运算

MATLAB的运算

符号表达式的处理

我们可以用符号运算工具推导数学公式，但结果可能不是最简形式，或者是我们期待的。因此我们要做一些化简处理。

调用函数simplify() 会调用所有化简函数对表达式进行简化。除了这个还有一些其他的函数可以用

格式 函数说明
[n,d] = numden(s) 提取符号表达式s的分子n和分母d
expand(s,fun) 表达式展开，可用于三角函数或者多项式
collect(s,s1) 按照s1表达式对s做同类项合并
v = factor(s) 对s进行因式分解
subs(s,s1,s2) 对s中的子表达式s1换成s2，并且可以拓展到更多个替换.subs(s,s1,s2,…,sn)
rewrite(s,fun) 对表达式s进行自定义改写。其中fun是自选的改写方式，例如’sin’表示化简结果中只含正弦函数,’sincos’表示只有正弦和余弦函数，也可以用一些其他的函数替代，例如exp,sqrt等等

例子：我们尝试吧$\tan(x+y)$改写为只有$\sin x$的表达式

命令：

1
2
3
4
5
6
7
syms x,y; F = tan(x+y)
F1 = rewrite(expand(F),'sin')
 
F1 =
 
(sin(x)/(2*sin(x/2)^2 - 1) + sin(y)/(2*sin(y/2)^2 - 1))/((sin(x)*sin(y))/((2*sin(x/2)^2 - 1)*(2*sin(y/2)^2 - 1)) - 1)

现在，给定$x\in R$，试着把表达式$y = \mid 2x^2-3\mid + \mid 4x-5\mid$表示为分段函数

MATLAB中，分段函数用piecewise表示，那么输入如下命令即可

1
2
3
syms x real; y(x) = abs(2*x^2-3)+abs(4*x-5)
y1 = rewrite(y,'piecewise')

1
2
3
4
5
6
>> y1 = rewrite(y,'piecewise')
 
y1(x) =
 
piecewise(5/4 <= x, 2*x^2 + 4*x - 8, x <= 5/4 & 0 <= 2*x^2 - 3, 2*x^2 - 4*x + 2, x <= 5/4 & 2*x^2 - 3 <= 0, - 2*x^2 - 4*x + 8)

也就是

循环结构

for循环

我们用一个简单的例子讲for的使用

求解$S = \sum^{100}_{i=1} i$可以这样做：

1
2
3
4
5
6

s = 0
for i=1:100
   s = s+i
   end
s

即可。

那么一些递推序列怎么办呢，比如假设$a_1 = 3$，第二项为$a_2 = \sqrt{1+a_1}$，此后的每一项都是$a_{k+1} = \sqrt{1+a_k}$，我们如何计算$a_{32}$呢

在计算机中，对于一般递推式我们使用称为迭代的方法，即先写出递推式的代码，然后循环变量k执行31次。

1
2
3
4
5
6

a =3; format long
for k=1:31
a=sqrt(1+a)
disp(a) %显示每一步的结果
end

结果：

1
2
3
4
5
6
>> a

a =

   1.618033988749895

while表达式

while和for不同的地方在于，while是带条件的循环，一旦条件表达式不成立则自动终止循环。

一样，我们使用while循环求解$\sum^{100}_{i=1}i$

1
2
3
4
5
s = 0;i =1
while(i <= 100)
s = s+i
end
s

现在，我们求出满足$s = \sum^m_{i=1}i > 100000$ 的最小m

1
2
3
4
5
s = 0; m=0
while( s <= 10000)
m = m+1;
s = s+m
end

向量化编程

循环结构的执行较慢，在实际编程中若能对矩阵或者向量运算时，尽量不要使用循环结构。而是使用向量化方法优化。

假设一组圆，半径为r = 1,1.2,0.9,0.7,0.85。试着求出这些圆的面积。

面积公式为$S = \pi r^2$。

若用for 循环，则需要先定义数组

1
2
3
4
r = [1,1.2,0.9,0.7,0.85]
for i:length(r)
S(i) = pi *r(i)^2
end

若使用向量化，则代码为

1
S = pi.*r.^2

函数的基本结构

现在，我们尝试写一些函数用来实现一些特定的功能。

假设我想写一个函数生成$n\times m$阶希尔伯特矩阵，第i行j列元素$h_{i,j} = 1/(i+j-1)$。

那么我们可以这样做，首先，我们的函数名字叫myhilb，它有两种调用格式，即H =myhilb(n,m) 和 H = myhilb(n)

前者生成矩阵，后者生成方阵

M-函数由function语句引导，基本结构为

1
function [out1,out2,...,outm] = fun_name(in1,in2,...,inn)

函数的输入和返回由nargin和nargout两个保留变量给出。在进入函数的时候自动生成，我们演示如何创建一个自己的函数。首先新建脚本， 名字为myhilb.m。接着在里面输入如下内容

1
2
3
4
5

function = A=myhilb(n,m)
if nargin == 1,m=n,end %若输入只有1个变量，则定义m=n。
for i = 1:n, for j = 1:m, A(i,j) = 1/(i+j-1);end,end

保存后即可在命令行界面用如下方法调用该函数

1
2
A1 = myhilb(4,3)
A2 = myhilb((sym(4)))

结果为：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

>> A1 = myhilb(4,3)

A1 =

   1.000000000000000   0.500000000000000   0.333333333333333
   0.500000000000000   0.333333333333333   0.250000000000000
   0.333333333333333   0.250000000000000   0.200000000000000
   0.250000000000000   0.200000000000000   0.166666666666667

>> A2 = myhilb(sym(4))
 
A2 =
 
[  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]

函数的递归调用

现在，我们编写一个求解阶乘n!的函数

核心思想是不停的调用函数直到停止，那么第k个变量为 k - n*my_fact(n-1)

一样的，新建脚本，名字my_fact.m，并写入如下代码

1
2
3
4
5
6
7
function k = my_fact(n)
if nargin ~= 1, error('Error:Only one input variable accepted'); end %判断有几个输入
if abs(n-floor(n)) > eps || n<0 
    error('you son of beach get me mtf negative integer');
end %判断输入是否是负数
if n> 1, k=n*my_fact(n-1); %递归调用
elseif any([0 1]==n),k=1;end %函数出口，即n=1的时候。

调用结果：

1
2
3
4
5
6
>> k = my_fact(5)

k =

   120

更多例子：若输入参数不确定个数，可以用varargin变量表示，使用varargin{i}提取第i个输入。

现在，我们编写一个加法mat_add，它接收任意多矩阵作为变量，并且输出加后的结果

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
function result = mat_add(varargin)
    if nargin == 1
       result = varargin{1};
    end
    if nargin == 0
        error('至少一个输入');
    end
    [rows,cols] = size(varargin{1});
    for i = 2:nargin
        [r,c]=size(varargin{i});
        if r~=rows || c ~= cols
            error('矩阵需要有相同的维度');
        end
    end
    result = zeros(rows,cols);
    for i = 1:nargin
        result = result +varargin{i};
    end
end




结果:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
A = [1,2,3;4,5,6;7,8,9]
B = [7,9,7;4,6,4;1,3,1]
>> S = mat_add(A)

S =

     1     2     3
     4     5     6
     7     8     9

>> S = mat_add(A,B)

S =

     8    11    10
     8    11    10
     8    11    10