MATLAB符号计算
文章目录
一、符号对象
1、符号对象的建立
(1)sym函数:用于建立单个符号对象,调用格式:
符号对象名=sym(A)
将由A来建立符号对象。其中,A可以是一个数值常量、数值矩阵或者数值表达式(不加单引号)。此时符号对象为一个符号常量;A也可以是一个变量名(加单引号),这时符号对象为一个符号变量。
对比一下程序:
由以上例子可以看出:
- 符号计算的结果是一个精确的数学表达式。
- 数值计算的结果是一个数值。
(2)syms命令
syms命令可以一次定义多个符号变量,其一般调用格式如下:
syms 符号变量名1 符号变量名2 … 符号变量名n
例如,要同时定义四个符号变量a、b、c、d,则命令如下:
.>>syms a b c d
2、符号对象的运算
(1)四则运算
符号表达式的四则运算与数值运算一样,用+、-、*、/、^运算符实现,其运算结果依然是一个符号表达式。
(2)关系运算
- 6种关系运算符:<、<=、>、>=、==、~=
- 对应的6个函数:lt()、le()、gt()、ge()、eq()、ne()
若参与运算的是符号表达式,其结果是一个符号表达式;若参与运算的是符号矩阵,其结果是由符号关系表达式组成的矩阵。
在进行符号对象的运算前,可用assume函数对符号对象设置值域,调用格式为:
assume(condition)
asume(expr,set)
第一种格式指定变量满足条件condition,第二种格式指定表达式expr属于集合set
(3)逻辑运算
- 3种逻辑运算符:&(与)、|(或)、~(非)
- 4个逻辑运算函数:and(a,b)、or(a,b)、not(a,b)、xor(a,b)
(4)因式分解与展开运算
- factor(s):对符号表达式s分解因式。
- expand(s):对符号表达式s进行展开。
- collect(s):对符号表达式s合并同类项。
- collect(s,v):对符号表达式s按变量v合并同类项
例子:梅森素数的验证问题
梅森数是指形如2p-1的一类整数,其中p是素数。将梅森数记为Mp。如果梅森数是素数,就称为梅森素数。
验证M19、M23、M29、M31是否为梅森素数。
(5)其他运算
- 提取有理分式的分子分母:[n,d]=numden(s)
- 提取符号表达式的系数:c=coeffs(s,x)
- 符号表达式化简:simplify(s)
- 符号多项式与多项式系数向量之间的转化:
– 符号多项式转换为多项式系数向量:p=sym2poly(s)
– 多项式系数向量转换为符号多项式:s=poly2sym( p)
例子:求方程
a
x
2
+
b
x
+
c
=
0
ax^2+bx+c=0
ax2+bx+c=0的根
注意此处不能调用sym2poly函数
(6)符号运算中变量的确定
- 如果没有明确指定自变量,MATLAB将按以下原则确定主变量并对其进行相应运算:
– 寻找除i、j之外,在字母顺序上最接近x的小写字母。
– 若表达式中有两个符号变量与x的距离相等,则ASCII码大的优先。 - symvar()函数可以用于查找一个符号表达式中的符号变量:symvar(s,n)
函数返回符号表达式s中的n个符号变量。因此,可以用symvar(s,1)查找表达式s的主变量。
3、符号矩阵
符号矩阵也是一种符号表达式,所以符号表达式运算都可以在矩阵意义下进行。注意这些函数作用于符号矩阵时,是分别作用于矩阵的每一个元素。
例如:建立符号矩阵并化简
[
a
3
−
b
3
s
i
n
2
α
+
c
o
s
2
α
15
x
y
−
3
x
2
x
−
5
y
78
]
\begin{bmatrix} a^3-b^3 &sin^2α+cos^2α\\ {15xy-3x^2}\over {x-5y} & 78 \end{bmatrix}
[a3−b3x−5y15xy−3x2sin2α+cos2α78]
例子:当λ取何值时,以下齐次线性方程组有非零解
{
(
1
−
λ
)
x
1
−
2
x
2
+
4
x
3
=
0
2
x
1
+
(
3
−
λ
)
x
2
+
x
3
=
0
x
1
+
x
2
+
(
1
−
λ
)
x
3
=
0
\left\{ \begin{aligned} (1-λ)x_1-2x_2+4x_3=0\\ 2x_1+(3-λ)x_2+x_3=0\\ x_1+x_2+(1-λ)x_3=0 \end{aligned} \right.
⎩⎪⎨⎪⎧(1−λ)x1−2x2+4x3=02x1+(3−λ)x2+x3=0x1+x2+(1−λ)x3=0
分析:对于齐次线性方程组,若有非零解的充要条件是其系数矩阵的行列式等于0或者其秩小于n。
因此,当lamda=0、2或者3时,系数矩阵行列式等于0。此时方程组有非零解。
二、符号微积分
1、符号函数的极限
(1)求符号函数极限的命令为limit,其调用格式为:
limit(f,x,a)
即求函数f关于变量x在a点的极限
(2)limit函数的另一种功能是求单边极限:
limit(f,x,a,‘right’)
limit(f,x,a,‘left’)
(3)例子:求下列极限
-
lim
x
→
a
x
m
−
a
m
x
−
a
\displaystyle\lim_{x \rightarrow a}{{\sqrt[m]{x}-\sqrt[m]{a}}\over{x-a}}
x→alimx−amx−ma
-
lim
n
→
∞
(
1
+
1
n
)
n
\displaystyle\lim_{n \rightarrow ∞}{(1+{1\over n})}^n
n→∞lim(1+n1)n
2、符号函数的导数
(1)diff(f,x,n)
即求函数f关于变量x的n阶导数
(2)例子:求下列函数的导数
-
y
=
1
+
e
x
,
求
y
′
y=\sqrt{1+e^x},求y'
y=1+ex,求y′
-
z
=
x
e
y
y
2
,
求
z
x
′
,
z
y
′
z={{xe^y}\over y^2},求z_x',z_y'
z=y2xey,求zx′,zy′
3、符号函数的积分
(1)不定积分
int(f,x),求函数f对变量x的不定积分
例子:①
∫
(
3
−
x
2
)
3
d
x
\int_{} (3-x^2)^3dx
∫(3−x2)3dx
②
∫
5
x
t
1
+
x
2
d
t
\int_{}{5xt\over {1+x^2}}dt
∫1+x25xtdt
(2)定积分
int(f,x,a,b),其中a、b分别表示定积分的下限和上限。
- 当函数f关于变量x在闭区间[a,b]可积,函数返回一个定积分结果。
- 当a、b中有一个是inf时,函数返回一个广义积分。
- 当a、b中有一个符号表达式,函数返回一个符号函数。
例子
①
∫
1
2
∣
1
−
x
∣
d
x
\int_{1}^2|1-x|dx
∫12∣1−x∣dx
②
∫
−
∞
+
∞
1
1
+
x
2
d
x
\int_{-∞}^{+∞}{1\over {1+x^2}}dx
∫−∞+∞1+x21dx
③
∫
2
s
i
n
x
4
x
t
d
t
\int_{2}^{sinx}{4x\over t}dt
∫2sinxt4xdt
4、例子
(1)第一问
xi=0:50:600;
yi=[4.4,4.5,4.6,4.8,4.9,5.1,5.4,5.2,5.5,5.2,4.9,4.8,4.7];
p=polyfit(xi,yi,3);
plot(xi,yi,'o',xi,polyval(p,xi));
syms x y;
y=poly2sym(p,x);
s=int(y,x,0,600);
v=s*0.6;
eval(v)
求得流量
v
=
1787.4
(
m
3
/
s
)
v=1787.4(m^3/s)
v=1787.4(m3/s)
(2)第二问
xi=0:50:600;
yi=[4.4,4.5,4.6,4.8,4.9,5.1,5.4,5.2,5.5,5.2,4.9,4.8,4.7];
p=polyfit(xi,yi,3);
syms y x;
y=poly2sym(p,x);
yii=diff(y,x);
x=50:60;
k=eval(yii);
all(abs(k)<1/1.5)
结果为1
三、级数
1、级数求和
(1)求无穷级数的和需要符号表达式求和函数symsum(s,v,n,m),其中s表示一个级数的通项,是一个符号表达式。v是求和变量,省略时使用系统的默认变量。n和m是求和变量v的初值和末值。
(2)例子:
- 求下列级数之和
① s 1 = 1 + 4 + 9 + 16 + … … + 10000 s_1=1+4+9+16+……+10000 s1=1+4+9+16+……+10000
② s 2 = 1 − 1 2 + 1 3 − 1 4 + … … + ( − 1 ) n + 1 1 n + … … s_2=1-{1\over 2}+{1\over 3}-{1\over 4}+……+(-1)^{n+1}{1\over n}+…… s2=1−21+31−41+……+(−1)n+1n1+……
③ s 3 = 1 − 1 3 + 1 5 − 1 7 + … … + ( − 1 ) n + 1 1 2 n − 1 + … … s_3=1-{1\over 3}+{1\over 5}-{1\over 7}+……+(-1)^{n+1}{1\over 2n-1}+…… s3=1−31+51−71+……+(−1)n+12n−11+……
- 银行利率计算问题
syms k r;
p2=eval(symsum(50000*(1+0.045/k)^k,k,2,2))%半年计息
p4=eval(symsum(50000*(1+0.045/k)^k,k,4,4))%按季度计息
p12=eval(symsum(50000*(1+0.045/k)^k,k,12,12))%按月计息
p_inf=limit(50000*(1+r/k)^k,k,inf)%用极限求无限期数
从结果知,当期数趋于无限大时,年利率确定,总金额也确定。因此,当年利率为4.5%,计息期趋于无穷时,一年后总金额为:
50000*exp(0.045)=5.2301e+04
(3)需要注意的是,符号计算中因为小数都是表示为有理分数的形式,随着计算次数的增加,容易导致分子或分母出现极大整数从而无法计算。
2、泰勒级数
(1)taylor(f,v,a,Name,Value)
该函数将函数f按变量v在a点展开为泰勒级数,v的默认值与diff函数相同,a的默认值为0。Name和Value为选项设置,前者为选项名,后者为该选项的值。
Name有3个可取字符串:
①’ExpansionPoint’:指定展开点,对应值可以是标量或向量,未设置时,展开点是0.
②’Order’:指定截断参数,对应值为一个正整数。未设置时,截断参数为6,即展开式的最高阶为5。
③’OrderMode’:指定展开式采用绝对阶或相对阶,对应值为’Absolute’或’Relative’。未设置时取’Absolute’。
(2)例子
①求函数f在x=1处的5阶泰勒级数展开式
f
(
x
)
=
1
+
x
+
x
2
1
−
x
+
x
2
f(x)={1+x+x^2\over1-x+x^2}
f(x)=1−x+x21+x+x2
(3)使用泰勒级数展开,可以让计算机只需要通过四则运算便能计算其他复杂函数。
例子:利用泰勒展开式计算三角函数cosx的值
四、符号方程的求解
1、代数方程符号求解
(1)MATLAB中求解符号表达式的代数方程可以用solve()实现:
- solve(s):求解符号表达式s的代数方程,求解变量为默认变量。
- solve(s,v):求解符号表达式s的代数方程,求解变量为v。
- solve(s1,s2,……,sn,v1,v2,……,vn):求解符号表达式s1,s2,……,sn组成的代数方程组,求解变量分别为v1,v2,……,vn。
(2)例子:解方程
a
x
2
+
b
x
+
c
=
0
ax^2+bx+c=0
ax2+bx+c=0
以下给出四种形式
(3)需要注意的是,solve()函数求解方程时不一定准确,有时候方程有解,solve却给出无解。
例如下面方程组:
{
x
+
y
=
98
x
3
+
y
3
=
2
\left\{ \begin{aligned} x+y=98\\ \sqrt[3]x+\sqrt[3]y=2 \end{aligned} \right.
{x+y=983x+3y=2
但是稍作变换,令
y
=
(
2
−
x
3
)
3
y=(2-\sqrt[3]x)^3
y=(2−3x)3,得到
x
+
(
2
−
x
3
)
3
=
98
x+(2-\sqrt[3]x)^3=98
x+(2−3x)3=98
因此,利用solve求解方程组需要尝试不同方法。
2、常微分方程符号求解
- MATLAB中,大写字母D表示导数。如,Dy表示y’,D2y表示y’’,Dy(0)=5表示y’(0)=5。但注意,此方法未来新版本替换,而是使用syms声明变量。例如: d 2 y d t 2 d^2y\over dt^2 dt2d2y,定义代码为:syms y(t);diff(y,t,2)
- 符号常微分方程求解可以通过函数dsolve()来实现,格式:dsolve(e,c,v)
用于求解常微分方程e在初值条件c下得特解。参数v是方阵中的自变量,省略时按默认处理,若没给出初值条件c,则求方程的通解。 - dsolve在求常微分方程组时格式为:
dsolve(e1,e2,……,en,c1,c2,……,cn,v)
(1)例子
①求下列微分方程或方程组的解
-
d
y
d
x
=
x
2
+
y
2
2
x
2
{dy\over dx}={x^2+y^2\over 2x^2}
dxdy=2x2x2+y2
-
{
d
y
d
x
=
4
x
−
2
y
d
y
d
x
=
2
x
−
y
\left\{ \begin{aligned} {dy\over dx}=4x-2y\\ {dy\over dx}=2x-y \end{aligned} \right.
⎩⎪⎪⎨⎪⎪⎧dxdy=4x−2ydxdy=2x−y
②疾病传染问题
(1)当λ=μ=0.3,b=0.7
syms y(t) t;
e1=diff(y,t)==0.3*y*(1-y)-0.3*y;
e2=y(0)==0.7;
y=dsolve(e1,e2,t)
fplot(y,[0,100])
xlabel('时间t');
ylabel('病人比例y')
当λ=μ时,病人比例随着时间增加而减少
(2)当λ=0.3,μ=0.15,b=0.15时
syms y(t) t;
e1=diff(y,t)==0.3*y*(1-y)-0.15*y;
e2=y(0)==0.3;
y=dsolve(e1,e2,t)
fplot(y,[0,100])
xlabel('时间t');
ylabel('病人比例y')
m=limit(y,t,inf)
此时,病人比例随着时间增加而增加,并趋于0.5