MATLAB符号运算

符号运算

定义符号变量

• 简单符号变量创建

syms x  % 观察工作区， sym是symbolic的缩写
syms a b c % 一次多定义几个符号
1
2

• 符号方程创建，两种方法

syms a x
y = a*x+x^2
% y = str2sym('a*x+x^2')  %Matlab 2017b 版本后推出
1
2
3

• 符号矩阵：矩阵中元素是符号

% 符号矩阵
syms alpha
M = [cos(alpha)  -sin(alpha);
       sin(alpha)  cos(alpha)]
1
2
3
4

简单运算

%% 简单运算 下面为得到的结果
syms a b c d e
y = a + b
% a + b
x = c - d
% c - d
y1 = x*y
% (a + b)*(c - d)
y2 = y1/y
% c - d
y3 = y1^3
% (a + b)^3*(c - d)^3
y4 = sqrt(y3)   
% ((a + b)^3*(c - d)^3)^(1/2)
y5 = exp(y4)
% exp(((a + b)^3*(c - d)^3)^(1/2))
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

符号方程/表达式的整理

• 化简
  • 不是所有的都能化简，太复杂的这个函数可能化简不了

% 化简
syms a
y=(cot(a/2)-tan(a/2))*(1+tan(a)*tan(a/2))
simplify(y)
% 2/sin(a)
1
2
3
4
5

• 因式分解

%% 因式分解
factor(12) % 对常数进行因式分解 
%      2     2     3
syms m n x
y = -24*m^2*x-16*n^2*x
factor(y)
% [ -8, x, 3*m^2 + 2*n^2]
y1=m^3-n^3
factor(y1)
% [ m - n, m^2 + m*n + n^2]
1
2
3
4
5
6
7
8
9
10

• 多项式展开

%% 多项式展开
syms a x
y = a*(x^2-a)^2+(x-2)
expand(y)
% a^3 - 2*a^2*x^2 + a*x^4 + x - 2
1
2
3
4
5

• 合并
  • 指定对哪个符号合并

%% 合并
syms x y
z = (x+y)^2*y+5*y*x-2*x^3
% expand(z)   - 2*x^3 + x^2*y + 2*x*y^2 + 5*x*y + y^3
collect(z,x)
% y*x^2 - 2*x^3 + (2*y^2 + 5*y)*x + y^3
collect(z,y)
% y^3 + 2*x*y^2 + (x^2 + 5*x)*y - 2*x^3
1
2
3
4
5
6
7
8

• 计算分子和分母
  • 函数的参数，必须是符号变量而不是数值变量
  • 对方程计算分子分母

% [z1,z2] = numden(2.5)  % 会报错，因为numden的输入变量不能是数值，只能是符号变量
% ans = sym(2.5);  % sym函数可以将数值2.5转换为符号
[z1,z2] = numden(sym(2.5)) % 对常数计算分子与分母
% z1 = 5
% z2 = 2
syms x y
z = 1/x*y+x/(x^2-2*y)
[z1,z2] = numden(z)  %z1分子，z2分母
% z1 = - x^2*y - x^2 + 2*y^2
% z2 = x*(- x^2 + 2*y)
1
2
3
4
5
6
7
8
9
10

• 让结果显示的更加自然
  • mupad代数工具箱
  • 未来的版本里面没有这个工具箱，可以新建一个实时脚本
  • 将复杂的式子复制进去，点击全部运行，就会生成公式，可以复制为图片也可以是latex

syms x y
M = (1/x*y+x/(x^2-2*y)-x^2/(3+y)^2)^2;
expand(M)  
% y^2/x^2 + x^4/(y^4 + 12*y^3 + 54*y^2 + 108*y + 81) + (2*x^3)/(- x^2*y^2 - 6*x^2*y - 9*x^2 + 2*y^3 + 12*y^2 + 18*y) - (2*y)/(- x^2 + 2*y) + x^2/(x^4 - 4*x^2*y + 4*y^2) - (2*x*y)/(y^2 + 6*y + 9)
mupad % 未来的版本可能会移除这个工具箱，可以点击Matlab的主页，新建实时脚本
1
2
3
4
5

符号函数求导与差分

diff()函数用作符号函数上是求导，用在矩阵上是差分

一元函数求导

syms x
y = x^4-5*x^2+6
diff(y) %求一阶导数
% 4*x^3 - 10*x
diff(y,2) %求二阶导数
% 12*x^2 - 10
1
2
3
4
5
6

• 一个现象
  • matlab求导有一套自己的规则，求出来的结果可能十分繁琐，所以这时候需要simplify化简
  • 例如，求cos(x)*tan(x)十阶导数，结果很长，但化简后只有sinx，很简单

y = cos(x)*tan(x)
dy = diff(y,10)  %求十阶导数
simplify(dy)
y = sin(x)*tan(x)
dy = diff(y,10)  %求十阶导数
simplify(dy)
1
2
3
4
5
6

多元函数求导（求偏导）

syms x1 x2 x3
y1 = x1^5*x2+x2*x3-x1^2*x3
py1 = diff(y1,x1,1) % 对x1求一阶偏导
% 5*x2*x1^4 - 2*x3*x1
py2 = diff(y1,x1,2) % 对x1求二阶偏导
% 20*x2*x1^3 - 2*x3
py3 = diff(y1,x1,x2) % 先对x1求偏导，再对x2求偏导
% 5*x1^4
py4 = diff(y1,x2,x1) % 先对x2求偏导，再对x1求偏导
% 5*x1^4


1
2
3
4
5
6
7
8
9
10
11
12

求差分

注意，如果diff函数作用的对象不是符号函数，而是矩阵，那么对应的功能是求差分

• 一阶差分二阶差分

A=[4 5 6 3 2 1];
diff(A)  % 求向量A的一阶差分     1     1    -3    -1    -1
diff(A,2)  % 在一阶差分的基础上再差分一次     0    -4     2     0
1
2
3

• 矩阵差分
  • 默认标识符是按行（列的方向上）差分，标识符为2时按列差分

A=[4 5 6; 
     7 4 2;
     5 6 2]
A1=diff(A)  % 下一行减去上一行求一阶差分
%      3    -1    -4
%     -2     2     0
A2=diff(A,2) % 下一行减去上一行求二阶差分（在一阶差分的基础上再差分一次）
%     -5     3     4
1
2
3
4
5
6
7
8

• 按行/列差分

A3=diff(A,2,1) % 最后面的1表示在行上进行差分（在列的方向上进行差分）
%     -5     3     4
A4=diff(A,1,2)  % 后一列减去前一列求一阶差分， 最后面的2表示在列上进行差分（在行的方向上进行差分）
%      1     1
%     -3    -2
%      1    -4
A4=diff(A,2,2) % 后一列减去前一列求二阶差分
%      0
%      1
%     -5
1
2
3
4
5
6
7
8
9
10

计算不定积分与定积分

计算不定积分

int(y,x)

积分英语：integral

y是被积函数，x是要对x进行积分

注意，Matlab计算时不会给我们加上常数C

syms x
y = x^2
int(y,x)
% x^3/3
1
2
3
4

syms x
y = 1/x
int(y,x)
% log(x)  

1
2
3
4
5

注意，Matlab计算1/x形式的不定积分时不会给我们加上绝对值~

• 关于对数

常见的有以2、e、10为底的求对数函数，如果说相求其他的，比如以6为底的，可以用换底公式

比如${\log }_{6}36=\frac{{\ln }_{}36}{{\ln }_6}$

syms x
y = x^2 / (1+x^2)
int(y,x)
% x - atan(x)

syms x
y = 1/(exp(x)+1)
int(y,x)
% x - log(exp(x) + 1)

syms x a
y = 1/sqrt(x^2-a^2)
int(y,x)
% log(x + (x^2 - a^2)^(1/2))
1
2
3
4
5
6
7
8
9
10
11
12
13
14

求不定积分的结果可能不同，可以求个导数验证一下

计算定积分

int(y,x,a,b)

a为下界，b为上界

syms x
y = sin(x)
int(y,x,0,pi/2) 
% 1

syms x a b
y = exp(x)
int(y,x,a,b)
% exp(b) - exp(a)
1
2
3
4
5
6
7
8
9

• 边界是无穷，用inf

syms x
y = (sin(x))^2 / x^2
b=int(y,x,0,+inf)
% pi/2

1
2
3
4
5

• 注意，不是所有的函数都可以利用int函数计算出最后的结果
  • 因为有些函数算不出来积分

syms x
y = 1 / exp(x) * log(x+2*x^2+sin(x))
int(y,x,0,4)

1
2
3
4

• 于是可以用数值积分

我们可以计算数值积分：数值积分可用于求定积分的近似值。在数值分析中，数值积分是计算定积分数值的方法和理论。
在数学分析中，给定函数的定积分的计算不总是可行的，许多定积分不能用已知的积分公式得到精确值。

integral()求数值积分

注意要用函数句柄

函数句柄与匿名函数知识

y = @(x) 1 ./ exp(x) .* log(x+2.*x.^2+sin(x))  % 注意，写成函数句柄时，要用点乘或者点除
integral(y,0,4)

% 主意这里不能再用y，因为y已经是一个函数句柄，指向了一个函数
xx = 0:0.1:4;
yy = 1 ./ exp(xx) .* log(xx+2*xx.^2+sin(xx));
plot(xx,yy,'-')
1
2
3
4
5
6
7

求解方程和方程组

不同的MATLAB版本之间的语法存在不兼容的情况：要会用帮助文档
% 视频里面用到的是Matlab2017a版本，如果低版本版本可能会报错。
% 更多关于Matlab求方程的介绍可看这个博客：求解方程方法

solve函数

求解单变量方程

要写两个等号：这里的等号一定要有两个，一个等号表示赋值，两个等号才表示左右两边相等

第一个参数是方程，后面是未知数

clear;clc
syms x
answ = solve(sin(x) == 1, x) 
answ = solve(sin(x) == 1)  % 只有一个符号变量x，所以可以不指定未知数
1
2
3
4

• 当方程很长的时候，可以先将方程赋给一个变量，然后再放入solve中

% 也可以这样写
clear;clc
syms x
eqn = (sin(x) == 1);  % eqn = sin(x) == 1;  
answ = solve(eqn, x)
% 因为三角函数是周期函数，如果要得到所有的解，则需要加上条件
[answ, params, condtions] = solve(eqn, x, 'ReturnConditions', true)
1
2
3
4
5
6
7

‘ReturnConditions’, true

这个参数意思是返回详细信息

三个返回值

answ是解

params是解里的变量

condtions是解里变量的条件，比如：整数

求解多变量方程

里面有不止1个未知数，在solve中指定要求解的未知数

clear;clc
syms a b c x
eqn = (a*x^2 + b*x + c == 0);
answ1 = solve(eqn, x)  % 将x视为未知数求解 
%  -(b + (b^2 - 4*a*c)^(1/2))/(2*a)
%  -(b - (b^2 - 4*a*c)^(1/2))/(2*a)
answ2 = solve(eqn, a) % 将a视为未知数求解
% -(c + b*x)/x^2
1
2
3
4
5
6
7
8

求解方程组

将几个方程放入矩阵组成方程组

求解返回一个结构体

用.取出各个变量的解

当然也可以用个数相等的[answ_u, answ_v]来接住答案，就不需要结构体了

clear;clc
syms u v a
eqn = [2*u + v == a, u - v == 1];
answ = solve(eqn, [u, v])
answ.u
answ.v
[answ_u, answ_v] = solve(eqn, [u, v])
1
2
3
4
5
6
7

solve 可能会警告

syms x
eqn = (sin(x) == x^2 - 1);
solve(eqn, x)  % 警告: Cannot solve symbolically. Returning a numeric approximation instead. 

1
2
3
4

• 告诉我们求解的结果可能不是很好

• fplot画表达式的图形，观察一下

% 画图看看 
fplot(sin(x), [-2 2])  % fplot函数可绘制表达式的图形
hold on
fplot(x^2 - 1, [-2 2]) 

1
2
3
4
5

• 可以看出结果有两个交点，而solve函数只给出了一个结果，同时还有警告

• 使用vpasolve函数

vpasolve函数

求解相比solve更加强大一下，同时可以定义求解的区间

% 用vpasolve函数指定求[0 2]上的解
syms x
eqn = sin(x) == x^2 - 1;
vpasolve(eqn, x, [0 2])
vpasolve(eqn, x, [-1 0])
vpasolve(eqn, x, [-10 10])
1
2
3
4
5
6

帮助文档

vpasolve returns all solutions only for polynomial equations. 
% For nonpolynomial equations, there is no general method of finding all solutions.
% When you look for numerical solutions of a nonpolynomial equation or system that has several solutions,
% then, by default, vpasolve returns only one solution, if any. 
% To find more than just one solution, set random to true. 
% Now, calling vpasolve repeatedly might return several different solutions.
1
2
3
4
5
6

vpasolve函数的帮助文档上面说，如果是多项式方程，那么该函数可以得到所有的解。如果是非多项式方程，可能只会返回一个解

我们画图知道这个方程有两个解，如果将区间定位[-10,10]（太大），那么可能就搜索不到解。按理来说有两个解，这里只返回了一个

如果要找到更多解的话，对参数’random’,设置true，多次重复调用，可能就会找到不同的解

vpasolve(eqn, x, 'random', true) 
vpasolve(eqn, x, -5)   % 给定搜索的起始点
1
2

因为指定了不同的初始搜索点，所以更可能获得不同的结果

一个复杂的例子

syms x y
eqn = [x^2 - 2*x - 3*x*y == 10, y^4 == exp(-2*x/3*y)]
[answ_x, answ_y] = vpasolve(eqn, [x, y], 'random', true)
1
2
3

• 多次运行可以获得多个结果

• 还是需要画图来查看有多少个解，大概在什么区间，可以缩小搜索范围

• fimplicit()函数画二元隐函数

  • 后面加个3可以绘制三维隐函数

fimplicit(x^2 - 2*x - 3*x*y == 10, [-10 10],'r')  % R2016b版本之后才有
hold on
fimplicit(y^4 == exp(-2*x/3*y), [-10 10],'b')  % R2016b版本之后才有
1
2
3

• 根据图指定搜索的范围，求解
• 同时可以在图上将结果标点

[answ_x, answ_y] = vpasolve(eqn, [x, y],[-4 -1;1 5])  % 指定搜索的范围：x位于[-4 -1], y位于[1 5]
hold on
plot(answ_x, answ_y,'ko', 'MarkerSize',10) 

1
2
3
4

• 这里answ_x, answ_y是符号变量
• 如果后续要参与数值运算的话，可以使用double()函数将其从符号转换为数值

fslove函数

最强大的

fsolve是Matlab优化工具箱中的一个函数，可专门用来求解特别复杂的方程和方程组

x0 = [0,0];  % 初始值
result_x = fsolve(@my_fun,x0)
1
2

• my_fun函数是自己定义的方程组
  • x指未知数，是一个向量，比如这里包括x1和x2
  • F是返回的方程组

function F = my_fun(x)
    F(1) = exp(-exp(-(x(1)+x(2)))) - x(2)*(1+x(1)^2);
    F(2) = x(1)*cos(x(2)) + x(2)*sin(x(1)) - 0.5;
end
1
2
3
4

• 当然可以用vpasolve函数试试

clear; clc
syms x1 x2
eqn =  [exp(-exp(-(x1+x2))) - x2*(1+x1^2) == 0, x1*cos(x2) + x2*sin(x1) - 0.5 == 0]
[answ_x1, answ_x2] = vpasolve(eqn, [x1, x2], [0 0])

1
2
3
4
5

内容来源：数学建模清风第四次直播：
利用matlab快速实现机器学习
配套的讲义和后续视频
欢迎大家学习