第 6章 MATLAB符号计算
6.1 符号计算基础
6.2 符号导数及其应用
6.3 符号积分
6.4 级数
6.5 代数方程的符号求解
6.6 常微分方程的符号求解
6.1 符号计算基础
6.1.1 符号对象
1,建立符号变量和符号常数
(1)sym函数
sym函数用来建立单个符号量,例如,a=sym('a')建立符号变量 a,此后,用户可以在表达式中使用变量 a进行各种运算。
例 6.1考察符号变量和数值变量的差别。
在 MATLAB命令窗口,输入命令:
a=sym('a');b=sym('b');c=sym('c');d=sym('d'); %定义 4个符号变量
w=10;x=5;y=-8;z=11; %定义 4个数值变量
A=[a,b;c,d] %建立符号矩阵 A
B=[w,x;y,z] %建立数值矩阵 B
det(A) %计算符号矩阵 A的行列式
det(B) %计算数值矩阵 B的行列式例 6.2比较符号常数与数值在代数运算时的差别。
在 MATLAB命令窗口,输入命令:
pi1=sym('pi');k1=sym('8');k2=sym('2');k3=sym('3'); % 定义符号变量
pi2=pi;r1=8;r2=2;r3=3; % 定义数值变量
sin(pi1/3) % 计算符号表达式值
sin(pi2/3) % 计算数值表达式值
sqrt(k1) % 计算符号表达式值
sqrt(r1) % 计算数值表达式值
sqrt(k3+sqrt(k2)) % 计算符号表达式值
sqrt(r3+sqrt(r2)) % 计算数值表达式值
(2)syms函数
syms函数的一般调用格式为:
syms var1 var2 … varn
函数定义符号变量 var1,var2,…,varn 等。用这种格式定义符号变量时不要在变量名上加字符分界符 ('),变量间用空格而不要用逗号分隔。
2,建立符号表达式例 6.3用两种方法建立符号表达式。
在 MATLAB窗口,输入命令:
U=sym('3*x^2+5*y+2*x*y+6') %定义符号表达式 U
syms x y; %建立符号变量 x,y
V=3*x^2+5*y+2*x*y+6 %定义符号表达式 V
2*U-V+6 %求符号表达式的值例 6.4计算 3阶范得蒙矩阵行列式的值。设 A是一个由符号变量 a,b,c确定的范得蒙矩阵。
命令如下:
syms a b c;
U=[a,b,c];
A=[[1,1,1];U;U.^2] %建立范得蒙符号矩阵
det(A) %计算 A的行列式值例 6.5建立 x,y的一般二元函数。
在 MATLAB命令窗口,输入命令:
syms x y;
f=sym('f(x,y)');
6.1.2 基本的符号运算
1,符号表达式运算
(1)符号表达式的四则运算例 6.6符号表达式的四则运算示例。
在 MATLAB命令窗口,输入命令:
syms x y z;
f=2*x+x^2*x-5*x+x^3 %符号表达式的结果为最简形式
f=2*x/(5*x) %符号表达式的结果为最简形式
f=(x+y)*(x-y) %符号表达式的结果不是 x^2-y^2,而是 (x+y)*(x-y)
(2)因式分解与展开
factor(S) 对 S分解因式,S是符号表达式或符号矩阵。
expand(S) 对 S进行展开,S是符号表达式或符号矩阵。
collect(S) 对 S合并同类项,S是符号表达式或符号矩阵。
collect(S,v) 对 S按变量 v合并同类项,S是符号表达式或符号矩阵。
例 6.7 对符号矩阵 A的每个元素分解因式。
命令如下:
syms a b x y;
A=[2*a^2*b^3*x^2-4*a*b^4*x^3+10*a*b^6*x^4,3*x*y-
5*x^2;4,a^3-b^3];
factor(A) %对 A的每个元素分解因式例 6.8 计算表达式 S的值。
命令如下:
syms x y;
s=(-7*x^2-8*y^2)*(-x^2+3*y^2);
expand(s) %对 s展开
collect(s,x) %对 s按变量 x合并同类项 (无同类项 )
factor(ans) % 对 ans分解因式
(3)表达式化简
MATLAB提供的对符号表达式化简的函数有:
simplify(S) 应用函数规则对 S进行化简。
simple(S) 调用 MATLAB的其他函数对表达式进行综合化简,并显示化简过程。
例 6.9化简命令如下:
syms x y;
s=(x^2+y^2)^2+(x^2-y^2)^2;
simple(s) %MATLAB自动调用多种函数对 s进行化简,并显示每步结果
2,符号矩阵运算
transpose(S) 返回 S矩阵的转置矩阵。
determ(S) 返回 S矩阵的行列式值。
colspace(S) 返回 S矩阵列空间的基。
[Q,D]=eigensys(S) Q返回 S矩阵的特征向量,D返回 S矩阵的特征值。
6.1.3 符号表达式中变量的确定
MATLAB中的符号可以表示符号变量和符号常数。
findsym可以帮助用户查找一个符号表达式中的的符号变量。该函数的调用格式为:
findsym(S,n)
函数返回符号表达式 S中的 n个符号变量,若没有指定 n,则返回 S中的全部符号变量。
在求函数的极限、导数和积分时,如果用户没有明确指定自变量,MATLAB将按缺省原则确定主变量并对其进行相应微积分运算。可用
findsym(S,1)查找系统的缺省变量,事实上,
MATLAB按离字符 'x'最近原则确定缺省变量。
6.2 符号导数及其应用
6.2.1函数的极限
limit函数的调用格式为:
limit(f,x,a)
limit函数的另一种功能是求单边极限,其调用格式为:
limit(f,x,a,'right') 或 limit(f,x,a,'left')
例 6.10求极限。
在 MATLAB命令窗口,输入命令:
syms a m x;
f=(x^(1/m)-a^(1/m))/(x-a);
limit(f,x,a) %求极限 (1)
f=(sin(a+x)-sin(a-x))/x;
limit(f) %求极限 (2)
limit(f,inf) %求 f函数在 x→∞( 包括 +∞
和 -∞)处的极限
limit(f,x,inf,'left') %求极限 (3)
f=(sqrt(x)-sqrt(a)-sqrt(x-a))/sqrt(x*x-a*a);
limit(f,x,a,'right') %求极限 (4)
6.2.2 符号函数求导及其应用
MATLAB中的求导的函数为:
diff(f,x,n)
diff函数求函数 f对变量 x的 n阶导数。参数 x的用法同求极限函数 limit,可以缺省,缺省值与 limit相同,n的缺省值是 1。
例 6.11求函数的导数。
命令如下:
syms a b t x y z;
f=sqrt(1+exp(x));
diff(f) %求 (1)。未指定求导变量和阶数,按缺省规则处理
f=x*cos(x);
diff(f,x,2) %求 (2)。求 f对 x的二阶导数
diff(f,x,3) %求 (2)。求 f对 x的三阶导数
f1=a*cos(t);f2=b*sin(t);
diff(f2)/diff(f1) %求 (3)。按参数方程求导公式求 y对 x的导数
(diff(f1)*diff(f2,2)-diff(f1,2)*diff(f2))/(diff(f1))^3 %求 (3)。求 y对 x的二阶导数
f=x*exp(y)/y^2;
diff(f,x) %求 (4)。 z对 x的偏导数
diff(f,y) %求 (4)。 z对 y的偏导数
f=x^2+y^2+z^2-a^2;
zx=-diff(f,x)/diff(f,z) %求 (5)。按隐函数求导公式求 z对 x的偏导数
zy=-diff(f,y)/diff(f,z) %求 (5)。按隐函数求导公式求 z对 y的偏导数例 6.12在曲线 y=x3+3x-2上哪一点的切线与直线
y=4x-1平行。
命令如下:
x=sym('x');
y=x^3+3*x-2; %定义曲线函数
f=diff(y); %对曲线求导数
g=f-4;
solve(g) %求方程 f-4=0的根,即求曲线何处的导数为 4
6.3 符号积分
6.3.1不定积分在 MATLAB中,求不定积分的函数是 int,其调用格式为:
int(f,x)
int函数求函数 f对变量 x的不定积分。参数 x可以缺省,缺省原则与 diff函数相同。
例 6.13求不定积分。
命令如下:
x=sym('x');
f=(3-x^2)^3;
int(f) %求不定积分 (1)
f=sqrt(x^3+x^4);
int(f) %求不定积分 (2)
g=simple(ans) %调用 simple函数对结果化简
6.3.2 符号函数的定积分定积分在实际工作中有广泛的应用。在 MATLAB中,定积分的计算使用函数:
int(f,x,a,b)
例 6.14求定积分。
命令如下:
x=sym('x');t=sym('t');
int(abs(1-x),1,2) %求定积分 (1)
f=1/(1+x^2);
int(f,-inf,inf) %求定积分 (2)
int(4*t*x,x,2,sin(t)) %求定积分 (3)
f=x^3/(x-1)^100;
I=int(f,2,3) %用符号积分的方法求定积分 (4)
double(I) %将上述符号结果转换为数值例 6.15求椭球的体积。
命令如下:
syms a b c z;
f=pi*a*b*(c^2-z^2)/c^2;
V=int(f,z,-c,c)
V =
4/3*pi*a*b*c
例 6.16轴的长度为 10米,若该轴的线性密度计算公式是
f(x)=6+0.3x千克 /米 (其中 x为距轴的端点距离 ),求轴的质量。
(1)符号函数积分。在 MATLAB命令窗口,输入命令:
syms x;
f=6+0.3*x;
m=int(f,0,10)
(2)数值积分。
先建立一个函数文件 fx.m:
function fx=fx(x)
fx=6+0.3*x;
再在 MATLAB命令窗口,输入命令:
m=quad('fx',0,10,1e-6)
例 6.17求空间曲线 c从点 (0,0,0)到点 (3,3,2)的长度。
求曲线 c的长度是曲线一型命令如下:
syms t;
x=3*t;y=3*t^2;z=2*t^3;
f=diff([x,y,z],t) %求 x,y,z对参数 t的导数
g=sqrt(f*f') %计算一型积分公式中的根式部分
l=int(g,t,0,1) %计算曲线 c的长度
6.3.3 积分变换
1,傅立叶 (Fourier)变换在 MATLAB中,进行傅立叶变换的函数是:
fourier(fx,x,t) 求函数 f(x)的傅立叶像函数 F(t)。
ifourier(Fw,t,x) 求傅立叶像函数 F(t)的原函数 f(x)。
例 6.18求函数的傅立叶变换及其逆变换。
命令如下:
syms x t;
y=abs(x);
Ft=fourier(y,x,t) %求 y的傅立叶变换
fx=ifourier(Ft,t,x) %求 Ft的傅立叶逆变换
2,拉普拉斯 (Laplace)变换在 MATLAB中,进行拉普拉斯变换的函数是:
laplace(fx,x,t) 求函数 f(x)的拉普拉斯像函数 F(t)。
ilaplace(Fw,t,x) 求拉普拉斯像函数 F(t)的原函数 f(x)。
例 6.19计算 y=x2的拉普拉斯变换及其逆变换,
命令如下:
x=sym('x');y=x^2;
Ft=laplace(y,x,t) %对函数 y进行拉普拉斯变换
fx=ilaplace(Ft,t,x) %对函数 Ft进行拉普拉斯逆变换
3,Z变换对数列 f(n)进行 z变换的 MATLAB函数是:
ztrans(fn,n,z) 求 fn的 Z变换像函数 F(z)
iztrans(Fz,z,n) 求 Fz的 z变换原函数 f(n)
例 6.20求数列 fn=e-n的 Z变换及其逆变换。
命令如下:
syms n z
fn=exp(-n);
Fz=ztrans(fn,n,z) %求 fn的 Z变换
f=iztrans(Fz,z,n) %求 Fz的逆 Z变换
4,积分变换的应用例 6.21用拉普拉斯方法解微分方程。
6.4 级数
6.4.1 级数的符号求和级数符号求和函数 symsum,调用格式为:
symsum(a,n,n0,nn)
例 6.22求级数之和。
命令如下:
n=sym('n');
s1=symsum(1/n^2,n,1,inf) %求 s1
s2=symsum((-1)^(n+1)/n,1,inf) %求 s2。未指定求和变量,缺省为 n
s3=symsum(n*x^n,n,1,inf) %求 s3。此处的求和变量 n不能省略。
s4=symsum(n^2,1,100) %求 s4。计算有限级数的和
6.4.2 函数的泰勒级数
MATLAB中提供了将函数展开为幂级数的函数 taylor,其调用格式为:
taylor(f,v,n,a)
例 6.23求函数在指定点的泰勒展开式。
命令如下:
x=sym('x');
f1=(1+x+x^2)/(1-x+x^2);
f2=sqrt(1-2*x+x^3)-(1-3*x+x^2)^(1/3);
taylor(f1,x,5) %求 (1)。展开到 x的 4次幂时应选择 n=5
taylor(f2,6) %求 (2)。
例 6.24将多项式表示成 x+1的幂的多项式。
命令如下:
x=sym('x');
p=1+3*x+5*x^2-2*x^3;
f=taylor(p,x,-1,4)
例 6.25应用泰勒公式近似计算 。
命令如下:
x=sym('x');
f=(1-x)^(1/12); %定义函数,4000^(1/12)=2f(96/2^12)
g=taylor(f,4) %求 f的泰勒展开式 g,有 4000^(1/12)≈2g(96/2^12)
b=96/2^12;
a=1-b/12-11/288*b^2-253/10368*b^3 %计算 g(b)
2*a %求 4000^(1/12)的结果
4000^(1/12) %用 MATLAB的乘方运算直接计算
6.4.3 函数的傅立叶级数
MATLAB 5.x版中,尚未提供求函数傅立叶级数的内部函数。下面我们自己设计一个简化的求任意函数的傅立叶级数的函数文件。
function mfourier=mfourier(f,n)
syms x a b c;
mfourier=int(f,-pi,pi)/2; %计算 a0
for i=1:n
a(i)=int(f*cos(i*x),-pi,pi);
b(i)=int(f*sin(i*x),-pi,pi);
mfourier=mfourier+a(i)*cos(i*x)+b(i)*sin(i*x);
end
return
调用该函数时,需给出被展开的符号函数 f和展开项数 n,不可缺省。
例 6.26在 [-π,π]区间展开函数为傅立叶级数。
命令如下:
x=sym('x');a=sym('a');
f=x;
mfourier(f,5) %求 f(x)=x的傅立叶级数的前 5项
f=abs(x);
mfourier(f,5) %求 f(x)=|x|的傅立叶级数的前 5项
syms a;
f=cos(a*x);
mfourier(f,6) %求 f(x)=cos(ax)的傅立叶级数的前 6项
f=sin(a*x);
mfourier(f,4) %求 f(x)=sin(ax)的傅立叶级数的前 4项
6.5代数方程的符号求解
6.5.1线性方程组的符号求解
MATLAB中提供了一个求解线性代数方程组的函数
linsolve,其调用格式为:
linsolve(A,b)
例 6.27求线性方程组 AX=b的解。
解方程组 (1)的命令如下:
A=[34,8,4;3,34,3;3,6,8];
b=[4;6;2];
X=linsolve(A,b) %调用 linsolve函数求 (1)的解
A\b %用另一种方法求 (1)的解解方程组 (2)的命令如下:
syms a11 a12 a13 a21 a22 a23 a31 a32 a33 b1 b2 b3;
A=[a11,a12,a13;a21,a22,a23;a31,a32,a33];
b=[b1;b2;b3];
X=linsolve(A,b) %调用 linsolve函数求 (2)的解
XX=A\b %用左除运算求 (2)的解
6.5.2 非线性方程组的符号求解求解非线性方程组的函数是 solve,调用格式为:
solve('eqn1','eqn2',…,'eqnN','var1,var2,…,varN')
例 6.28 解方程。
命令如下:
x=solve('1/(x+2)+4*x/(x^2-4)=1+2/(x-2)','x') %解方程 (1)
f=sym('x-(x^3-4*x-7)^(1/3)=1');
x=solve(f) %解方程 (2)
x=solve('2*sin(3*x-pi/4)=1') %解方程 (3)
x=solve('x+x*exp(x)-10','x') %解方程 (4)。仅标出方程的左端例 6.29】 求方程组的解。
命令如下:
[x y]=solve('1/x^3+1/y^3=28','1/x+1/y=4','x,y') %解方程组 (1)
[x y]=solve('x+y-98','x^(1/3)+y^(1/3)-2','x,y') %解方程组 (2)
Warning,Explicit solution could not be found.
> In C:\MATLABR11\toolbox\symbolic\solve.m at line 136
x =
[ empty sym ]
y =
[]
对方程组 (2)MATLAB给出了无解的结论,显然错误,请看完全与其同构的方程组 (3)。输入命令如下:
[u,v]=solve('u^3+v^3-98','u+v-2','u,v') %解方程组 (3)
[x v]=solve('x^2+y^2-5','2*x^2-3*x*y-2*y^2') %解方程组 (4)
6.6常微分方程的符号求解
MATLAB的符号运算工具箱中提供了功能强大的求解常微分方程的函数 dsolve。该函数的调用格式为:
dsolve('eqn1','condition','var')
该函数求解微分方程 eqn1在初值条件 condition下的特解。参数 var描述方程中的自变量符号,省略时按缺省原则处理,若没有给出初值条件 condition,则求方程的通解。
dsolve在求微分方程组时的调用格式为:
dsolve('eqn1','eqn2',…,'eqnN','condition1',…,'conditionN','var1',…,'varN')
函数求解微分方程组 eqn1,…,eqnN在初值条件 conditoion1,…,
conditionN下的解,若不给出初值条件,则求方程组的通解,
var1,…,varN给出求解变量。
例 6.30 求微分方程的通解。
命令如下:
y=dsolve('Dy-(x^2+y^2)/x^2/2','x') %解 (1)。方程的右端为 0时可以不写
y=dsolve('Dy*x^2+2*x*y-exp(x)','x') %解 (2)
y=dsolve('Dy-x/y/sqrt(1-x^2)','x') %解 (3)
例 6.31 求微分方程的特解。
命令如下:
y=dsolve('Dy=2*x*y^2','y(0)=1','x') %解 (1)
y=dsolve('Dy-x^2/(1+y^2)','y(2)=1','x') %解 (2)
例 6.32用微分方程的数值解法和符号解法解方程,并对结果进行比较。
在 MATLAB命令窗口,输入命令:
y=dsolve('Dy+2*y/x-4*x','y(1)=2','x') %用符号方法得到方程的解析解为了求方程的数值解,需要按要求建立一个函数文件 fxyy.m:
function f=fxyy(x,y)
f=(4*x^2-2*y)/x; %只能是 y'=f(x,y)的形式,当不是这种形式时,要变形。
return
输入命令:
[t,w]=ode45('fxyy',[1,2],2); %得到区间 [1,2]中的数值解,以向量 t,w存储。
为了对两种结果进行比较,在同一个坐标系中作出两种结果的图形。输入命令:
x=linspace(1,2,100);
y=x.^2+1./x.^2; %为作图把符号解的结果离散化
plot(x,y,'b.',t,w,'r-');
6.6.3 常微分方程组求解例 6.33 求微分方程组的解。
命令如下:
[x,y]=dsolve('Dx=4*x-2*y','Dy=2*x-y','t') %解方程组 (1)
[x,y]=dsolve('D2x-y','D2y+x','t') %解方程组 (2)