Matlab符号計算

在科學研究和工程應用中，除了存在大量的數值計算外，還有對符号對象進行的運算，即直接對抽象的符号對象進行的計算，并将所得到的結果以标準的符号形式來表示。符号計算可以得到比數值計算更一般的結果。Matlab的符号計算是通過集成在Matlab中的符号運算工具箱（Symbolic Math Toolbox）來實現的。本章主要介紹符号計算基礎，符号微積分，級數的符号求和，代數方程和微分方程的符号求解等内容。

一．符号計算基礎

Matlab提供了一種符号數據類型，相應的運算對象稱為符号對象。如：符号常量，符号變量，以及它們參與的數學表達式等。在進行符号運算前首先要建立符号對象，然後才可以進行符号對象的運算。

一．符号對象

1． 建立符号變量和符号常量

matlab提供了兩個建立符号對象的命令：sym和syms。其用法不同。

①sym函數

sym函數用來建立單個符号量，格式為：

符号量名=sym(符号字符串)

該函數可以建立一個符号字符串，符号字符串可以是常量、變量、函數或表達式。如：a=sym(' a ')将建立符号變量a,此後用戶可以在表達式中使用變量a進行各種運算.符号變量a和在其他過程中建立的非符号變量a是不同的。一個非符号變量在參與運算前必須賦值，變量的運算實際上該變量所對應值的運算，其運算結果是一個和變量類型對應的值，而符号變量參與運算前無須賦值，其結果是一個由參與運算的變量名組成的表達式。

a=sym('a');％定義符号變量

b=sym('b');

c=sym('c');

x=5;％定義數值變量

y=-8;

z=11;

w=a*a+b*b+c*c％符号運算

w =

a^2+b^2+c^2

w=x*x+y*y+z*z％數值運算

w=

210

使用sym函數還可以定義符号常量，使用符号常量進行代數運算時和數值常量進行的運算不同。下面比較差别：

pi1=sym('pi');k1=sym('8');k2=sym('3'); ％定義符号變量

pi2=pi;r1=8;r2=3;％定義數值變量

sin(pi1/3)％符号計算

ans =

1/2*3^(1/2)

sin(pi2/3)％數值計算

ans =

0.8660

sqrt(k1+sqrt(k2))％符号計算

ans =

(8+3^(1/2))^(1/2)

sqrt(r1+sqrt(r2))％數值計算

ans =

3.1196

從命令執行結果來看，用符号常量進行計算像在進行數學演算，所得到的結果是精确的數學表達式，而數值計算将結果近似為一個有限的小數。

②syms命令

sym函數一次隻能定義一個符号變量，使用不方便。Matlab提供了syms命令，一次可以定義多個符号變量。其格式：

syms arg1 arg2 … argn

arg1 arg2 … argn是定義的符号變量名，注意變量間用空格而不要用逗号分隔。

2．建立符号表達式

含有符号對象的表達式稱為符号表達式，建立符号表達式有以下3種方法：

①利用單引号生成符号表達式

y='1/sqrt(2*x)'％生成一般的符号表達式

y =

1/sqrt(2*x)

f='cos(x^2)-sin(2*x)=0'％生成符号方程

f =

cos(x^2)-sin(2*x)=0

②用sym函數建立符号表達式

U=sym('3*x^2-5*y+2*x*y+6')％生成符号表達式，不需定義變量

U =

3*x^2-5*y+2*x*y+6

M=sym('[a,b;c,d]')％生成符号矩陣

M =

[ a, b]

[ c, d]

③使用已定義的符号變量組成符号表達式

syms x y;

V=3*x^2-5*y+2*x*y+6

V =

3*x^2-5*y+2*x*y+6

二．基本的符号運算

1．四則運算：

符号表達式的加減乘除可以分别利用函數symadd、symsub、symmul、symdiv來實現，幂運算可以由sympow來實現。

例：

f=‘2*x^2+3*x-5’%定義符号表達式

g=‘x^2-x+7’

U=symadd(f,g)%求f+g

V=symsub(f,g) %求f-g

W=symmul(f,g)%求f*g

X=symdiv(f,g)%求f/g

Y=sympow(f,’3*x’)%求f^(3x)

另外，與數值運算一樣，也可以用+ - * / ^運算符來實現符号運算。如：

syms x y z;

f=2*x+x^2*x-5*x+x^3+exp(2)

f=2*x/(5*x)

f=(x*x-y*y)/(x-y)

f =

(x^2-y^2)/(x-y)%有時，Matlab并未将結果化為最簡形式。

2．符号表達式的提取分子和分母運算

如果符号表達式是一個有理分式或可以展開為有理分式，可利用numden函數來提取符号表達式中的分子、分母。調用格式

[n,d]=numden(s)

該函數提取符号表達式s的分子和分母，分别存入n和d中

如a=sym(0.3333)

a =

3333/10000

[n,d]=numden(a)

n =

3333

d =

10000

再如：f=sym('a*x^2/(b+x)')

f =

a*x^2/(b+x)

[n,d]=numden(f)

n =

a*x^2

d =

b+x

3．符号表達式的因式分解和展開

Matlab提供了符号表達式的因式分解和展開的函數，

factor(s)，對符号表達式s進行分解因式

expand(s)，對s進行展開

collect(s)，對s進行合并同類項

collect(s,v)，對s按變量v合并同類項。

例

syms a b x y;

A=a^3-b^3;

factor(A)

ans =

(a-b)*(a^2+a*b+b^2)

s=(-7*x^2-8*y^2)*(-x^2+3*y^2)

expand(s)％對s展開

ans =

7*x^4-13*x^2*y^2-24*y^4

collect(s,x)％對s按變量x合并同類項

ans =

7*x^4-13*x^2*y^2-24*y^4

4．符号表達式的化簡

Matlab提供的對符号表達式化簡的函數有：

simplify(s) ，應用函數規則對s進行化簡

simple(s)，調用Matlab的其他函數對表達式進行綜合化簡，并顯示化簡過程。

syms x y a

s=log(2*x/y);

simplify(s)

ans =

log(2)+log(x/y)

s=(-a^2+1)/(1-a)

simplify(s)

ans =

a+1

函數simple試用幾種不同的化簡工具，然後選擇在結果中含有最少字符的那種形式。如下例：

syms x y;

s=(x^2+y^2)^2+(x^2-y^2)^2

simple(s)％下面是自動調用多種函數對s進行化簡，并顯示每步結果，以下是化簡過程

simplify:

2*x^4+2*y^4

radsimp:

2*x^4+2*y^4

combine(trig):

2*x^4+2*y^4

factor:

2*x^4+2*y^4

expand:

2*x^4+2*y^4

combine:

(x^2+y^2)^2+(x^2-y^2)^2

convert(exp):

(x^2+y^2)^2+(x^2-y^2)^2

convert(sincos):

(x^2+y^2)^2+(x^2-y^2)^2

convert(tan):

(x^2+y^2)^2+(x^2-y^2)^2

collect(x):

2*x^4+2*y^4

ans =

2*x^4+2*y^4

5．符号表達式與數值表達式之間的轉換

利用函數sym可以将數值表達式變換為它的符号表達式：

sym(1.5)

ans =

3/2

函數numeric或eval可以将符号表達式變換成數值表達式：

phi='(1+sqrt(5))/2'

phi =

(1+sqrt(5))/2

numeric(phi)

ans =

1.6180

三．符号表達式中變量的确定

Matlab中的符号可以表示符号變量和符号常量，findsym可以幫助用戶查找一個符号表達式中的符号變量。該函數的調用格式為：

findsym(s,n)

函數返回符号表達式s中的n個符号變量，若沒有指定n則返回s中的全部符号變量。

syms x a y z b;%定義5個符号變量

s1=3*x+y;s2=a*y+b%定義兩個符号表達式

findsym(s1)

findsym(s2,2)

syms x y;

s=2*x+3*y;

findsym(s)

ans =

x, y

syms a b x y;％定義符号變量

c=sym('3');％定義符号常量c

findsym(a*x+b*y+c)

ans =％c不在結果中出現

a, b, x, y

注：Matlab按離字母x最近原則确定默認變量。

四．符号矩陣

二．符号函數及其應用

一．符号函數的極限

Matlab中求極限的函數是limit，可以用來求函數在指定點的極限值和左右極限。對于沒有定義的極限，Matlab給出的結果為NaN，極限值為無窮大時，Matlab給出的結果為inf。limit的調用格式：

①limit(f,x,a)求符号函數f(x)的極限。當x趨向于a時，f(x)的極限值。

②limit(f,a) 求符号函數f(x)的極限。由于沒有指定自變量，則使用該格式時，符号函數f(x)的變量為函數findsym(f)确定的默認自變量，即變量x趨向于a。

③limit(f) 求符号函數f(x)的極限。沒有指定變量的目标值，系統默認變量趨近于0時的情況。

④limit(f,x,a,’right’)，求極限，’right’表示變量x從右邊趨近于a。

⑤limit(f,x,a,’left’)，求極限，’left’表示變量x從左邊趨近于a。

例：求下列極限

syms a m x;

f=(x^(1/m)-a^(1/m))/(x-a);

limit(f,x,a)

f=(sin(a+x)-sin(a-x))/x;

limit(f)

f=x*(sqrt(x^2+1)-x);

limit(f,x,inf,’left’)

f=(sqrt(x)-sqrt(a)-sqrt(x-a))/sqrt(x*x-a*a);

limit(f,x,a,’right’)

二．符号函數求導及其應用

diff函數用于對符号表達式求導數，其調用格式為：

①diff(f)沒有指定變量和導數階數，則系統按findsym(f)函數指示的默認變量對符号表達式求一階導數

②diff(f,v)以v為自變量對符号表達式f求一階導數

③diff(f,n)對f求n階導數，n為正整數

④diff(f,v,n) 以v為自變量，對f求n階導數。

例：求下列函數的導數

1．

2． 求

3． 求

4． 求

5． 求

syms x

f=sqrt(1+exp(x));

diff(f)%對默認變量求一階導數

f=x*cos(x)

diff(f,x,2) %對x求二階導數

diff(f,x,3) %對x求三階導數

三．符号積分

一．符号函數的不定積分

在Matlab中，int函數用于求符号函數的不定積分，有兩種格式：

int(f)，沒有指定積分變量和積分階數時，系統按findsym函數指示的默認變量對被積函數求不定積分

int(f,v)，以v為自變量對被積函數求不定積分

例，求下列不定積分

1．

2，

3．

4．

命令如下：

x=sym('x')

f=(3-x^2)^3;

int(f)

f=sin(x)^2

int(f)

syms alpha t;

f=exp(alpha*t);

int(f)

f=5*x*t/(1+x^2)

int(f,t)

二．符号函數的定積分

在Matlab中，求符号函數的定積分也是使用int函數，調用格式為：

int(f,v,a,b)

其中，a,b分别表示定積分的上下限。該函數求被積函數f在區間[a,b]上的定積分。a,b可以是兩個具體的數，也可以是一個符号表達式，還可以是無窮(inf)。當函數f關于變量x在閉區間[a,b]上可積時，返回一個定積分結果。當a,b有一個是inf時，返回一個廣義積分。當ab有一個是符号表達式時，函數返回一個符号函數。

例：求下列定積分

1．

2．

3．

4．

syms x t;

int(abs(1-x),1,2)

f=1/(1+x^2);

int(f,-inf,inf)

f=x^3/(x-1)^10;

I=int(f,2,3)

double(I)% 把結果轉換為數值

int(4*x/t,t,2,sin(x))

例 求橢球的體積

用平面Z=z0，去截上述橢球，其相交線是一個橢圓，該橢圓在xy平面投影的面積是：

，橢球的體積為：

例：

三．積分變換

積分變換就是通過積分運算把一個函數f（原函數）變成另一個函數F（像函數），變換過程是：

其中K(x,t)稱為變換的核，變換的核決定了變換的不同名稱。一一對應，相互轉化。

積分變換的一項基本應用是求解微分方程，求解過程是基于這樣這樣一種想法：。。。

常見的積分變換有：傅立葉變換、拉普拉斯變換和Z變換

1．傅立葉（Fourier）變換

當積分變換的核 （i為虛數單位），稱積分變換

為傅立葉變換，其逆變換為：

在Matlab中，進行傅立葉變換的函數是

fourier(f,x,t)；求f的傅立葉像函數F(t)

ifourier(F,t,x)；求傅立葉像函數F(t)的原函數f(x)

例：求函數y=|x|的傅立葉變換及其逆變換

命令如下：

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的拉普拉斯變換的像函數F(t)

ilaplace(Fw,t,x) ：求拉普拉斯變換的像函數F(t)的原函數f（x）

例：求 的拉普拉斯變換及其逆變換。

命令如下：

x=sym(‘x’);

y=x^2;

Ft=laplace(y,x,t)%對y進行拉普拉斯變換

Fx=ilaplace(Ft,t,x) %進行逆變換

3．Z變換

當函數f(x)呈現為一個離散的數列f(n)時，稱變換為：

為Z變換，其逆變換為： （i為虛數單位）

對數列f(n)進行Z變換的Matlab函數是：

ztrans(fn,n,z) ：求fn的Z變換像函數F(z)

iztrans(Fz,z,n) ：求Fz的Z變換原函數f(n)

例： 求數列 的Z變換和逆變換。

syms n z;

fn=exp(-n);

Fz=ztrans(fn,n,z)% 求Z變換

f=iztrans(Fz,z,n)%求逆變換

四．級數

一．級數的符号求和

（sum處理的級數是有窮級數）

對無窮級數求和，sum是無能為力的，需要使用符号表達式求和函數symsum，調用格式為：

symsum(a,v,m,n)

其中a表示一個級數的通項，是一個符号表達式。v是求和變量，m，n是求和的開始項和末項。（無窮inf）

例：求下列級數之和

1．

2．

3．

4．

命令如下：

n=sym(‘n’)

s1=symsum(1/n^2,n,1,inf)

s2=symsum((-1)^(n+1)/n,1,inf)

…………

二．函數的泰勒級數

泰勒（Taylor）級數将一個任意函數表示為一個幂級數，并且，在許多情況下，隻需要取幂級數的前有限項來表示該函數，這對于大多數工程應用來說，精度已經足夠。Matlab提供了taylor函數将函數展開為幂級數，調用格式為

taylor(f,v,n,a)

該函數将f按變量v展開為泰勒級數，展開到第n項（即變量v的n-1次幂）為止。n的默認值為6，v的默認值與diff函數相同，參數a指定将函數f在自變量v=a處展開，a的默認值為0。

例：求函數的泰勒級數展開

1．求 的5階泰勒級數展開

2．将 在x=1處按5次多項式展開（n=6）

命令如下：

x=sym(‘x’);

f1=sqrt(1-2*x+x^3)-(1-3x+x^2)^(1/3)

f2=(1+x+x^2)/(1-x+x^2)

taylor(f1,x,5)%求(1)

taylor(f2,6,1)%求(2)，展開到x-1的5次幂時應選擇n=6

五．符号方程求解

一．符号代數方程求解

代數方程是指未涉及微積分運算的方程，相對比較簡單。在Matlab中，求解用符号表達式表示的代數方程可由函數solve實現，調用格式為：

①solve(eq) 求解符号表達式表示的代數方程eq，求解變量為默認變量，當方程右端為0時，方程eq中可以不包含右端項和等号，而僅列出方程左邊的表達式。

②solve(eq,v) 求解符号表達式表示的代數方程eq，求解變量為v

③solve(eq1,eq2,…,v1,v2,…) 求解符号表達式eq1，eq2…組成的代數方程組，求解變量為v1,v2,v3,….

例：解下列方程

1．

2．

3．

4．

命令如下：

x=solve(‘1/(x+2)+4*x/(x^2-4)=1+2/(x-2)’,’x’)

f=sym(‘x-(x^3-4*x-7)^(1/3)=1’)

x=solve(f)

x=solve(‘2*sin(3*x-pi/4)=1’)

x=solve(‘x-x*exp(x)-10’,’x’)

例：求下列方程組的解

1．

2．

3．

4．

[x y]=solve(‘1/x^3+1/y^3=28’,’1/x+1/y=4’,’x,y’)

[x y]=solve(‘x+y-98’,’x^(1/3)+y^(1/3)-2’,’x,y’)

二．符号常微分方程求解

在Matlab中，用大寫字母D表示導數，如Dy表示y的導數，D2y表示y的2階導數等。符号常微分方程求解可以通過函數dsolve來實現，調用格式：

dsolve(eq,c,v)

該函數求解常微分方程eq在初值條件c下的特解。v是方程中的自變量，可以省略。若沒有給出初值條件c則求方程的通解。

dsolve在求常微分方程組時的調用格式：

dsolve(eq1，eq2，eq3，…，c1，c2，c3，…，v1，v2，v3，…)

該函數求解常微分方程組eq1，eq2，eq3，…在初值條件c1，c2，c3，…下的特解，若不給出初值條件則求方程的通解。

例：

1．求 的通解

2．求 的通解

3．求 的特解，y(2)=1

4．求 的通解

5．求 的通解。

命令如下：

y=dsolve(‘Dy-(x^2+y^2)/x^2/2’,’x’)

y=dsolve('Dy*x^2+2*x*y-exp(x) ', 'x')

y=dsolve('Dy-x^2/(1+y^2) ', 'y(2)=1', 'x')

[x,y]=dsolve('Dx=4*x-2*y','Dy=2*x-y','t')

[x,y]=dsolve('D2x-y', 'D2y+x', 't')