第5章MATLAB符号计算



 

 


符号计算又称计算机代数，通俗地说就是用计算机推导数学公式，如对表达式进行因式分解、化简、微分、积分、解代数方程、求解常微分方程等。与数值计算不同，符号计算对操作对象不进行离散化和近似化处理，因此其计算结果是完全准确而没有误差的。
MATLAB是通过调用第三方软件进行符号计算的，在MATLAB R2008b(MATLAB 7．7)之前，MATLAB符号计算引擎是基于Maple内核的，从MATLAB R2008b开始，符号计算引擎开始采用MuPAD内核。对于普通的MATLAB用户，符号计算引擎的改变并没有带来什么不同，用户只需要调用MATLAB函数进行符号运算即可。本章结合具体案例介绍MATLAB符号计算。
5．1符号对象和符号表达式
符号对象和符号表达式是MATLAB进行符号计算的基本元素，要进行符号计算首先要创建符号对象，对符号对象进行数学上的运算操作便得到符号表达式。
5．1．1创建符号对象
MATLAB中提供了sym和syms函数，用来创建符号对象，其用法参见下例。

【例5．11】调用sym和syms函数创建符号对象。



>> a = sym('6.01');% 定义符号常数

>> b = sym('b'，'real'); % 定义实数域上的符号变量

>> A = ［1， 2; 3， 4］; % 定义数值矩阵

>> B = sym(A);% 把数值矩阵转为符号矩阵

>> C = sym('c%d%d'，［3，4］)% 定义3行4列的符号矩阵

C =

［ c11， c12， c13， c14］

［ c21， c22， c23， c24］

［ c31， c32， c33， c34］

>> symsxy % 同时定义多个复数域上的符号变量

>> symszpositive % 定义正实数域上的符号变量

>> symsf(x，y)% 定义符号函数

>> f(x，y) = x + y^2; % 指定符号函数表达式






>> c = f(1， 2)% 计算符号函数在(1，2)点处的函数值

c =

5

>> zv = solve(z^2 == 1， z)% 求方程z^2 = 1的解(只有大于零的解)

zv =

1

>> symsz% 撤销对符号变量取值域的限定，将其恢复为复数域上的符号变量

>> zv = solve(z^2 == 1， z)% 求方程z^2 = 1的解(有两个解)

zv =

1

 -1



由例5．11可知，sym函数可用来创建单个符号对象，而syms函数可用来创建多个符号对象。关于sym和syms函数的详细用法，读者可以参考其帮助文档。
5．1．2符号变量取值域的限定
从例5．11可以看出，MATLAB中定义的符号变量有取值域的限定，不同的取值域会对符号计算的结果造成不同影响，默认情况下均为复数域。在调用sym或syms函数定义符号变量时，可通过附加real或positive改变符号变量的取值域，除此之外，还可调用assume和assumeAlso函数对符号变量取值域实行更为精细的控制。
【例5．12】解不等式x2＞12，其中x为整数变量，并满足0＜x＜5。



>> syms x% 定义符号变量

>> assume(x>0 & x<5);% 对符号变量的取值域进行限定，0<x<5

>> assumeAlso(x，'integer');% 对符号变量的取值域增加别的限定，x取整数

>> assumptions(x)% 查看符号变量取值域的限定

ans =

［ in(x， 'integer')， 0 < x， x < 5］

>> result = solve(x^2>12) % 求解不等式

result =

4

>> symsx% 重新定义符号变量，以撤销对符号变量取值域的限定



若需撤销对符号变量取值域的限定，将其恢复为复数域上的符号变量，只需用命令“syms 变量名”重新定义符号变量即可。
【说明】在MATLAB R2018a及其以前的版本中，若需撤销对符号变量取值域的限定，将其恢复为复数域上的符号变量，可用命令“syms 变量名 clear”实现，单纯用命令“clear 变量名”只是从MATLAB工作区中删除了变量，并没有改变符号计算引擎MuPAD中变量取值域的限定。
5．1．3创建符号表达式
对符号对象进行各种运算(算术运算、关系运算、逻辑运算)，即可创建符号表达式。
1． 符号计算中的算术运算与转置
由于MATLAB采用了重载(overload)技术，使得用来构成符号表达式的运算符，无论在拼写还是在使用方法上，都与数值计算中的运算符完全相同。常用的算术运算符与转置符如表5．11所示。


表5．11符号计算中的算术运算符与转置符



运算符
说明
运算符
说明

+
加
．*
点乘
-
减
．/
右点除
*
乘
．＼
左点除
/
右除
．^
点乘方
＼
左除
'
共轭转置
^
乘方
．'
非共轭转置


【例5．13】在符号对象的基础上通过算术运算创建符号表达式。



>> syms a b c x y z% 定义多个符号变量

>> f1 = a*x^2+b*x-c;% 创建符号表达式f1

>> f2 = sin(x)*cos(y);% 创建符号表达式f2

>> f3 = (x+y)/z;% 创建符号表达式f3

>> f4 = ［x+1， x^2; x^3， x^4］ % 创建符号表达式矩阵f4

f4 =

［ x + 1， x^2］

［ x^3， x^4］

>> f5 = f4'% 符号表达式矩阵的共轭转置(')

f5 =

［ conj(x) + 1， conj(x)^3］

［ conj(x)^2， conj(x)^4］

>> f6 = f4.'% 符号表达式矩阵的转置(.')

f6 =

［ x + 1， x^3］

［x^2， x^4］



2． 符号计算中的关系运算
符号计算中的关系运算符如表5．12所示。


表5．12符号计算中的关系运算符



运算符
说明
运算符
说明

==
等于
＜=
小于或等于
＞=
大于或等于
＜
小于
＞
大于
~=
不等于

3． 符号计算中的逻辑运算
符号计算中的逻辑运算符如表5．13所示。


表5．13符号计算中的逻辑运算符



运算符
说明

|
逻辑或
&
逻辑与
~
逻辑非
xor
逻辑异或

【例5．14】在符号对象的基础上通过算术运算、关系运算和逻辑运算创建符号表达式。



>> syms x y% 定义符号变量

>> f1 = abs(x) >= 0 % 创建符号表达式

f1 =

0 <= abs(x)

>> f2 = x^2 + y^2 == 1 % 创建符号表达式

f2 =

x^2 + y^2 == 1

>> f3 = ~(y - sqrt(x) > 0) % 创建符号表达式

f3 =

~0 < y - x^(1/2)

>> f4 = x > 0 | y < -1 % 创建符号表达式

f4 =

0 < x | y < -1

>> f5 = x > 0 & y < -1 % 创建符号表达式

f5 =

0 < x & y < -1



通过算术运算、关系运算和逻辑运算创建的符号表达式是由逻辑运算符连接的一系列等式或不等式。利用isAlways函数或logical函数可判断等式或不等式是否成立，利用isequaln函数可判断两个符号表达式是否相等，它们均返回逻辑型结果。
【例5．15】isAlways、logical和isequaln函数的用法示例。



>> syms x % 定义符号变量

>> f = abs(x) >= 0;% 创建符号表达式

>> result1 = isAlways(f)% 判断不等式|x|>=0是否成立

result1 =

1

>> result2 = isequaln(abs(x)， x)% 判断|x|是否等于x

result2 =

0

>> assume(x>0);% 限定x>0

>> result3 = isequaln(abs(x)， x)% 重新判断|x|是否等于x

result3 =

1

>> syms x % 撤销对符号变量取值域的限定



5．1．4符号表达式的常用运算
这里要介绍的符号表达式的常用运算包括： 因式(或因子)分解、合并同类项、对指定项展开、提取符号多项式系数、提取分式的分子和分母、对符号表达式进行化简、对符号分式进行约分、复合函数、嵌套多项式等，用到的MATLAB函数如表5．14所示。


表5．14符号表达式的常用运算函数



函数名
说明
函数名
说明

factor
因式(或因子)分解
horner
嵌套多项式
collect
合并同类项
simplify
对符号表达式进行化简
expand
对指定项展开
simplifyFraction
对符号分式进行约分续表


函数名
说明
函数名
说明

combine
把相同的代数结构结合在一起
compose
复合函数
coeffs
提取符号多项式系数
finverse
反函数
numden
提取分式的分子和分母
pretty
美化显示符号表达式

【例5．16】因式(或因子)分解。
(1) 对f=x3－y3进行因式分解； (2)对符号数12345678901234567890 进行因子分解。



>> syms x y % 定义符号变量

>> f = factor(x^3-y^3)% 对符号表达式进行因式分解

f =

［ x - y， x^2 + x*y + y^2］

>> fa = factor(sym('12345678901234567890'))% 对符号数进行因子分解

fa =

［ 2， 3， 3， 5， 101， 3541， 3607， 3803， 27961］




【例5．17】把f=(x+y)(x2+y2+1)按变量y合并同类项。



>> syms x y % 定义符号变量

>> f = (x+y)*(x^2+y^2+1); % 定义符号表达式

>> collect(f，y)% 对符号表达式按变量y合并同类项

ans =

y^3 + x*y^2 + (x^2 + 1)*y + x*(x^2 + 1)



【例5．18】把符号表达式展开。
(1) f1=cos(x+y)； (2)f2=(a+b)e(a－b)2。



>> syms x y a b% 定义符号变量

>> f = ［cos(x+y); (a+b)*exp((a-b)^2)］;% 定义符号表达式向量

>> expand(f)% 把符号表达式展开

ans =

cos(x)*cos(y) - sin(x)*sin(y)

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


【例5．19】对符号表达式进行化简。
(1) f1=4x2+4x+1；(2)f2=cos(3arccos(x))。



>> syms x % 定义符号变量

>> f1 = sqrt(4/x^2+4/x+1); % 定义符号表达式

>> g1 = simplify(f1)% 按默认设置进行化简

g1 =

((x + 2)^2/x^2)^(1/2)

>> g2 = simplify(f1，'IgnoreAnalyticConstraints'，1)% 忽略分析约束进行化简

g2 =

(x + 2)/x

>> pretty(g2)% 把符号表达式显示为数学公式形式

x + 2

-----





x

>> f2 = cos(3*acos(x)); % 定义符号表达式

>> g3 = simplify(f2， 'Steps'， 4)% 进行4步化简

g3 =

4*x^3 - 3*x



【说明】
pretty函数用来美化符号表达式在屏幕上的显示样式，把符号表达式显示为数学公式的形式。在MATLAB R2006b及其以后的版本中，用户可利用实时编辑器窗口进行符号计算，结果展示更为直观，如图5．11所示。



图5．11MATLAB实时编辑器窗口



【例5．110】已知f(x)=ex，g(x)=sinx，求复合函数y=f(g(x))，并计算 f(g(π))。



>> syms f(x) g(x)% 定义符号函数

>> f(x) = exp(x);% 定义符号函数表达式

>> g(x) = sin(x);% 定义符号函数表达式



>> y1(x) = f(g(x))% 求复合函数(方法1)

y1(x) =

exp(sin(x))



>> y2(x) = compose(f，g) % 求复合函数(方法2)

y2(x) =

exp(sin(x))



>> y = y1(pi) % 计算复合函数在x = pi处的函数值

y =

1



【例5．111】求符号函数f(x)=esinx的反函数。



>> syms f(x)% 定义符号函数

>> f(x) = exp(sin(x));% 定义符号函数表达式

>> g(x) = finverse(f(x))% 求符号函数的反函数

g(x) =

asin(log(x))




5．1．5符号运算中的转换操作
1． 符号数与数值型数(或字符)的转换

当利用符号计算得到计算结果时，有时需要将其转化为数值型结果，以便在后续数值计算中加以利用。符号数与数值型数(或字符)的转换函数如表5．15所示。


表5．15符号数与数值型数(或字符)的转换函数



函数名
说明
函数名
说明

sym
创建符号对象
int8， int16， 

int32， int64
把符号矩阵转为有符号整型矩阵
double
把符号矩阵转为双精度矩阵
uint8， uint16， 

uint32，uint64
把符号矩阵转为无符号整型矩阵
eval
执行MATLAB运算
poly2sym
根据系数向量得到符号多项式
single
把符号矩阵转为单精度矩阵
sym2poly
根据符号多项式得到系数向量
vpa
按指定的有效数字位数来显示符号数值对象
char
把符号对象转为字符串


【例5．112】计算符号函数f(x)=ln(5．2)ex在x=3处的函数值。



>> syms f(x) % 定义符号函数

>> f(x) = log(sym(5.2))*exp(x);% 指定符号函数表达式

>> y = f(3)% 计算符号函数在x = 3处的函数值

y =

exp(3)*log(26/5)

>> y1 = double(y)% 把符号数转为双精度数

y1 =

33.1142

>> y2 = vpa(y，10)% 以10位有效数字形式显示符号数

y2 =

33.1141937

>> x = 3;% 指定x的值

>> y3 = eval(f) % 执行MATLAB运算，得到函数值

y3 =

33.1142



2． 符号表达式中的变量替换
MATLAB中提供了subs函数，用来对符号表达式中的变量(或符号项)进行替换。
【例5．113】把符号表达式f=asin(x)+b中的a和b分别换为2和5，并把sin(x)换为ln(y)。



>> syms a b x y% 定义符号变量

>> f = a*sin(x)+b; % 定义符号表达式

>> f1 = subs(f，sin(x)，log(y))% 符号项替换

f1 =

b + a*log(y)

% 变量替换方式一





>> f2 = subs(f1，［a，b］，［2，5］) % 同时替换变量a和b的值

f2 =

2*log(y) + 5

% 变量替换方式二

>> f3 = subs(f1，{a，b}，{2，5}) % 同时替换变量a和b的值

f3 =

2*log(y) + 5



【例5．114】把符号表达式f=asin(x)+b中的a和b分别换为2和5，并求f在x=1，2，3三点处的值。



>> syms a b x % 定义符号变量

>> f = a*sin(x)+b;% 定义符号表达式

>> y = subs(f， {a，b，x}， {2， 5， 1:3})% 同时替换多个符号变量的值

y =

［ 2*sin(1) + 5， 2*sin(2) + 5， 2*sin(3) + 5］

>> y = double(y) % 将计算结果转为双精度值

y =

6.68296.81865.2822



3． 将符号表达式转为函数
当通过符号计算得到一个表达式时，希望把它转化成关于其中某个变量的函数，这里的函数可以是符号函数，也可以是匿名函数或M文件函数。MATLAB中的symfun函数用来将符号表达式转为符号函数，matlabFunction函数用来将符号表达式转为匿名函数或M文件函数。
【例5．115】把符号表达式asin(x)+b转为关于x的符号函数f(x)，并求其在x=1，2，3三点处的函数值。



>> syms a b x

>> f(x) = symfun(a*sin(x)+b， x); % 把符号表达式转为符号函数

>> y = f(1:3)

y =

［ b + a*sin(1)， b + a*sin(2)， b + a*sin(3)］



【例5．116】创建符号表达式f=a(x+b)c+d，并将a，b，c，d分别替换为2，－1，1/2，3，并将替换后表达式转为x的函数f(x)，然后求函数在x=10处的函数值。



>> syms a b c d x% 定义符号变量

>> f = a*(x+b)^c+d; % 定义符号表达式

>> g = subs(f，{a，b，c，d}，{2，-1，sym(1/2)，3});% 同时替换多个变量

>> FunFromSym1 = matlabFunction(g)% 将符号表达式转为匿名函数

FunFromSym1 = 

@(x)sqrt(x-1.0).*2.0+3.0

>> y = FunFromSym1(10)% 调用匿名函数计算函数值

y =

9

% 将符号表达式转为M文件函数FunFromSym2.m

>> matlabFunction(g，'file'，［pwd，'＼FunFromSym2.m'］，...

'vars'，{'x'}，'outputs'，{'y'});

>> y = FunFromSym2(10)% 调用M文件函数计算函数值

y =

9



本例中由符号表达式转化成的M文件函数FunFromSym2．m的代码如下： 



function y = FunFromSym2(x)

%FUNFROMSYM2

%Y = FUNFROMSYM2(X)



%This function was generated by the Symbolic Math Toolbox version 8.5.

%12-Apr-2021 09:58:54



y = sqrt(x-1.0).*2.0+3.0;



5．1．6符号函数绘图
根据符号表达式或符号函数进行绘图的MATLAB函数如表5．16所示。


表5．16根据符号表达式或符号函数进行绘图的MATLAB函数



函数名
说明
函数名
说明

fplot
绘制平面曲线
ezpolar
绘制极坐标函数图

fplot3
绘制三维曲线
ezplot
绘制平面曲线

fmesh
绘制三维网格图
ezplot3
绘制三维曲线

fsurf
绘制三维面图
ezmesh
绘制三维网格图

fcontour
绘制等高线图
ezsurf
绘制三维面图

fimplicit
绘制平面隐函数图
ezcontour
绘制等高线图

fimplicit3
绘制三维隐函数图
ezcontourf
绘制填充式等高线图

ezmeshc
绘制带有等高线的三维网格图


ezsurfc
绘制带有等高线的三维面图


注： 以f开头的函数是MATLAB较新版本里提供的绘图函数，以ez开头的函数是MATLAB较老版本中提供的绘图函数。



【例5．117】绘制函数f(x)=1ln|x|在x∈［－6，6］上的图形。



>> syms f(x) % 定义符号函数

>> f(x) = 1/log(abs(x));% 指定符号函数表达式

>> fplot(f，［-6，6］); % 绘制函数图形

>> xlabel('x'); % 添加x轴标签

>> ylabel('$$ f(x) = ＼frac{1}{ln|x|} $$'，'Interpreter'，'Latex');% 添加y轴标签



以上命令做出的函数图形如图5．12所示。


【例5．118】绘制函数f(x，y)=xyx2+y2，x2+y2≠0
0，x2+y2=0在－1≤x，y≤1上的图形。



>> syms f(x，y) % 定义符号函数

>> f(x，y) = x*y/(x^2+y^2);% 指定符号函数表达式

>> fsurf(f，［-1，1，-1，1］)% 绘制函数图形

>> xlabel('x');ylabel('y');zlabel('z');% 添加坐标轴标签

>> hold on % 开启图形保持

>> plot3(0，0，0，'r.'，'MarkerSize'，20)% 绘制原点

>> title('$$f(x，y)=＼frac{xy}{x^2+y^2}$$'，'Interpreter'，'latex') % 添加标题

>> view(127，34)% 设置视点位置



以上命令做出的函数图形如图5．13所示。



图5．12一元符号函数图形




图5．13二元符号函数图形



5．2符号微积分


MATLAB符号计算功能强大，可以解决高等数学中大多数微积分问题，而且求解命令简单，符合人们求解问题的思路。

5．2．1极限、导数和级数的符号计算
MATLAB中完成极限、导数和级数符号计算的函数如表5．21所示。


表5．21求极限、导数和级数的MATLAB函数



函数名使 用 示 例说明

limit


limit(f，x，a)求极限limx→afx
limit(f，x，a，'left')求左极限limx→a－fx
limit(f，x，a，'right')求右极限limx→a+fx


diff
diff(f，x，n)求n阶导数f(n)(x)
jacobianjacobian(f，v)求多元向量函数fv 的Jacobian矩阵
taylortaylor(f， x， x0， 'Order'， n)把fx在x=x0处作n－1阶泰勒展开
symsum
symsum(f，k，a，b)求∑bk=afk


【例5．21】求下列极限。
(1)limn→∞(－1)n(n+1)2； (2)limx→0－sinaxax； (3)limx→∞1－2xkx； 
(4)limx→b
y→ca1+x2+y2。



>> syms n a b c k x y% 定义符号变量

>> xn = (-1)^n/(n+1)^2;% 数列通项

>> L1 = limit(xn，n，inf)% 求数列极限





L1 =

0



>> f1 = sin(a*x)/(a*x);% 定义符号表达式

>> L2 = limit(f1，x，0，'left')% 求函数极限

L2 =

1



>> f2 = (1-2/x)^(k*x); % 定义符号表达式

>> L3 = limit(f2，x，inf)% 求函数极限

L3 =

exp(-2*k)



>> f3 = a/(1+x^2+y^2); % 定义符号表达式

>> L4 = limit(limit(f3，x，b)，y，c)% 求函数极限

L4 =

a/(b^2 + c^2 + 1)



【例5．22】求显函数和隐函数的导数。
(1) 已知f(x)=sin2x，求f′(1)和f″(x)； (2)已知cos(x+siny)=siny，求dydx。



>> syms x y% 定义符号变量

>> f(x) = sin(x)^2; % 定义一元符号函数

>> df = diff(f，x)% 求一阶导函数

df(x) =

2*cos(x)*sin(x)



>> df_1 = df(1)% 求一阶导数值

df_1 = 

2*cos(1)*sin(1)

 

>> ddf = diff(f，x，2)% 求二阶导函数

ddf(x) = 

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

 

>> F(x，y) = cos(x+sin(y))-sin(y); % 定义二元符号函数

>> dy = -diff(F，x)/diff(F，y)% 隐函数求导

dy(x， y) = 

-sin(x + sin(y))/(cos(y) + sin(x + sin(y))*cos(y))



【例5．23】求fx1，x2=x1+x2
x2lnx1的Jacobian矩阵： f1x1f1x2
f2x1f2x2。



>> syms x1 x2% 定义符号变量

>> f = ［x1+x2;x2*log(x1)］;% 定义符号表达式向量

>> v = ［x1;x2］;% 定义自变量向量

>> jac = jacobian(f，v)% 求Jacobian矩阵

jac = 

［ 1，1］

［ x2/x1， log(x1)］



【例5．24】求函数f(x)=ex在x=0处展开的5阶Maclaurin公式（函数在x=0处的泰勒公式称为函数的Maclaurin公式）。



>> syms x % 定义符号变量

>> f = exp(x);% 定义符号表达式

>> g = taylor(f， x， 0， 'Order'， 6)% 泰勒展开

g = 

x^5/120 + x^4/24 + x^3/6 + x^2/2 + x + 1



【例5．25】求下列无穷级数。
(1) ∑+∞k=3k－22k； (2)∑+∞k=11(2k+1)2，－1k3k。



>> syms k % 定义符号变量

>> f1 = (k-2)/2^k;% 级数的一般项

>> s1 = symsum(f1，k，3，inf) % 级数求和

s1 = 

1/2



>> f2 = ［1/(2*k+1)^2，(-1)^k/3^k］;% 定义向量形式的一般项

>> s2 = symsum(f2，k，1，inf) % 级数求和

s2 = 

［ pi^2/8 - 1， -1/4］



5．2．2符号积分计算
MATLAB中提供了int函数，用来求符号函数的积分，其调用格式如下： 



intf = int(f，x)	% 求以x为自变量的函数f的不定积分

intf = int(f，x，a，b)	% 求以x为自变量的函数f从a到b的定积分




【例5．26】求下列积分。
(1) ∫xln(ax)dx； (2)∫1－11－x2dx； (3)∫+∞－∞e－x22dx； (4)∫21∫2xx∫2xyxyx+yzdzdydx。



>> syms x y z a% 定义符号变量

>> F = int(x*log(a*x)，x)% 求不定积分

F = 

(x^2*(log(a*x) - 1/2))/2

 

>> f1 = sqrt(1-x^2);% 定义符号表达式

>> s1 = int(f1，x，-1，1)% 求定积分

s1 = 

pi/2

 

>> f2 = exp(-x^2/2);% 定义符号表达式





>> s2 = int(f2，x，-inf，inf)% 求反常积分

s2 = 

2^(1/2)*pi^(1/2)

 

>> f3 = (x+y)/z;% 定义符号表达式

>> s3 = int(int(int(f3，z，x*y，2*x*y)，y，x，2*x)，x，1，2) % 求多重积分

s3 = 

(35*log(2))/6

 

>> s4 = double(s3) % 符号数转为双精度数

s4 =

4.0434



5．3符号方程求解

5．3．1符号代数方程求解
MATLAB中提供了solve函数，用来求解代数方程的符号解。solve函数不仅可以求解单个线性/非线性方程，而且还可以求解线性/非线性方程组。它的调用格式如下： 



 ［y1，…，yN］ = solve(eqns，vars，Name，Value)



【例5．31】求解下列方程。
(1) x3－2x2+4x=8； (2)sinx+cos2x=1； (3)x+xex－10=0。



>> syms x

>> Result1 = solve(x^3 - 2*x^2 + 4*x == 8， x)% 求解多项式方程

Result1 = 

2

 -2i

2i

 

>> Result2 = solve(sin(x) + cos(2*x) == 1， x)% 求解三角函数方程

Result2 = 

0

pi/6

 (5*pi)/6

 

% 求解三角函数方程，并返回参数及其条件

>> ［Result3，params，conditions］ = solve(sin(x) + cos(2*x) == 1，...

 x， 'ReturnConditions'，true) 



Result3 =% 方程的所有解，其中带有参数k

pi*k

pi/6 + 2*pi*k

 (5*pi)/6 + 2*pi*k 

 

params =% 参数k

k 






 conditions =% 参数k满足的条件(k均为整数)

 in(k， 'integer')

 in(k， 'integer')

 in(k， 'integer')

 

>> Result4 = solve(x + x*exp(x) == 10， x)% 求解超越方程

警告: Unable to solve symbolically. Returning a numeric solution using vpasolve. 

 

Result4 =% 没有找到符号解，返回最接近的数值解

1.6335061701558463841931651789789



【说明】
对于本例，在求解三角函数方程时，方程的解通常具有周期性，默认情况下，solve函数只返回一个周期上的解，可将ReturnConditions参数值设为true，使其返回所有解。在求解第三个方程时，solve函数返回了一个警告信息，这是因为solve函数没有找到符号解，此时返回最接近的数值解。
【例5．32】求解方程组x－3+y－3=28
x－1+y－1=4。



>> syms x y% 定义符号变量

>> ［X，Y］ = solve(［1/x^3 + 1/y^3 == 28， 1/x + 1/y == 4］， ［x，y］)% 求解方程组

X = 

1

1/3 

 

Y = 

1/3

1



5．3．2符号常微分方程求解
MATLAB中求解符号常微分方程的函数是dsolve，其调用格式如下： 



［y1，…，yN］ = dsolve(eqns，conds，Name，Value)



【例5．33】求解常微分方程d2ydx2=x+y。



>> syms y(x)% 定义符号函数y(x)

>> Y = dsolve(diff(y，2) == x+y)% 求解常微分方程

Y = 

C2*exp(x) - x + C1*exp(-x)


上述结果中的C1、C2都是任意常数。
【例5．34】求解初值问题： dydt=1+y2，y0=1。



>> syms y(t)% 定义符号函数y(t)

>> Y = dsolve(diff(y) == 1 + y^2， y(0) == 1)% 求解初值问题





Y = 

tan(t + pi/4)

 

% 考虑分析上的约束，求解更具一般意义的解

>> Y = dsolve(diff(y) == 1 + y^2， y(0) == 1，…

 'IgnoreAnalyticConstraints'， false)

Y = 

piecewise(in(C1， 'integer')， tan(t + pi/4 + pi*C1))



【说明】本例演示了dsolve函数的两种调用方法，其中第二种调用用到了IgnoreAnalyticConstraints参数，其字面意思是“忽略分析上的约束”，这是出于一些求解结果在一般性上的考虑。譬如下面表达式： lnex，通常我们认为其等于x，这是基于默认x为实数情况下得出的结论。如果x=2πi，则lne2πi=ln1=0而不是等于2πi。上述IgnoreAnalyticConstraints有两个参数值可供选择： true和false。默认情况下为true，此时不对所求结果进行一般意义上的推广，所求解在最一般意义条件下可能不成立，但仍满足原始微分方程以及定解条件。而如果选择false，dsolve返回的解(前提是能够求得解析解)在最一般意义下也会成立，但是会增加求不出统一的解析表达式的概率。
【例5．35】求解两点边值问题： xy''－3y'=x2，y1=0，y5=0，并绘制解曲线。



>> syms y(x)% 定义符号函数

% 求解两点边值问题

>> Y = dsolve(x*diff(y，2)-3*diff(y) == x^2， ［y(1) == 0， y(5) == 0］)

Y = 

(31*x^4)/468 - x^3/3 + 125/468



>> h = fplot(Y，［-1，6］);% 绘制解曲线

>> set(h，'color'，'k'，'LineWidth'，2，'LineStyle'，'--'); % 设置解曲线的属性

>> hold on; % 开启图形保持

>> plot(［1 5］，［0，0］，'p'，'color'，'r'，'markersize'，12); % 画微分方程的两个边值点

>> text(1，1，'y(1) = 0');% 图上标注边值条件

>> text(4，1，'y(5) = 0');

>> title('');% 设置标题为空

>> hold off;% 关闭图形保存


上述代码得到的图形如图5．31所示。


图5．31两点边值问题的解曲线



【例5．36】求解两点边值问题： xy″－3y′=x2，y1=0，y′5=1。



>> syms y(x)% 定义符号函数

>> Dy = diff(y);% 求符号函数的一阶导函数

>> eq = x*diff(y，2)-3*Dy == x^2; % 定义符号微分方程

>> Y = dsolve(eq， ［y(1) == 0， Dy(5) == 1］)% 求解常微分方程

Y =

(13*x^4)/250 - x^3/3 + 211/750



【例5．37】求解常微分方程组dxdt=y
dydt=－x。



>> syms x(t)y(t) % 定义符号函数

>> ［X， Y］ = dsolve(diff(x) == y， diff(y) == -x)% 求解常微分方程组

X = 

C1*cos(t) + C2*sin(t) 

 

Y = 

C2*cos(t) - C1*sin(t)



上述结果中的C1、C2都是任意常数。

5．4建模案例选讲——药物中毒的施救方案

5．4．1问题描述

一名儿童一次性误服11片治疗哮喘病的、剂量为100mg/片的氨茶碱片，到达医院就诊时，距服药时刻已过去两个小时，此时孩子已出现呕吐、头晕等不良症状。按照药品使用说明书，氨茶碱的成人用量是每次100~200mg，儿童用量是3~5mg/kg，如果过量服用，可使血药浓度(单位容积血液中的药量)过高，当血药浓度达到100μg/ml时，会出现严重中毒，达到200μg/ml则可致命。接诊医生需要根据专业知识迅速判断孩子的中毒状况，选择合理的施救方案。
5．4．2问题分析
人体服用药物后的吸收代谢过程如图5．41所示，药物首先进入胃肠道，再由胃肠道的外壁进入血液循环系统，然后随血液循环被人体脏器代谢而排出体外。血液系统对药物的吸收率(药物从胃肠道向血液系统转移的转移率)通常与胃肠道中的药量成正比，血液系统对药物的排除率(药物从血液系统向体外转移的转移率)通常与血液中的药量成正比。


图5．41口服药物的吸收代谢过程



接诊医生需要根据药物的吸收代谢规律及时计算出孩子入院时的血药浓度，以此选择接下来的施救方案。如果血药浓度达到危险的水平，可选的临床施救方案包括： 口服活性炭和体外血液透析，前者通过口服活性炭来吸附药物，可使血液中药物的排除率增至原来(人体自身)的2倍，后者通过在体外对血液进行过滤，可使血液中药物的排除率增至原来的6倍。
5．4．3模型假设与符号约定
为了判断孩子的中毒状况，需要建立数学模型，求出胃肠道和血液系统中药量随时间变化的规律。这里把t时刻胃肠道中的药量记为x(t)，血液系统中的药量记为y(t)，以孩子误服药物的时刻为t=0时刻。在此基础上给出如下模型假设。
(1) 由于药物从胃肠道向血液系统转移的转移率与胃肠道中的药量x(t)成正比，这里记比例系数为λ(λ＞0)，并假设总剂量为1100mg的药物在t=0时刻全部进入胃肠道，即x(0)=1100。
(2) 由于血液系统中药物的排除率与血液系统中的药量y(t)成正比，这里记比例系数为μ(μ＞0)，并假设t=0时刻血液系统中药量为0，即y(0)=0。
(3) 假设任意t时刻，血液系统中的药物分布是均匀的。
(4) 在药物总量确定的情况下，血药浓度与血液总量有关。人体血液总量约占体重的7%~8%(L/kg)，即每千克体重约有70~80ml血液，一个体重为50~60kg的成年人的血液总量为4000ml左右。假设误服药物儿童的血液总量为成人的一半，约为2000ml。
(5) 比例系数λ和μ可由氨茶碱吸收和排除的半衰期确定，从药品说明书可知，氨茶碱吸收的半衰期约为5h，即x(5)=x(0)/2，排除的半衰期约为6h，若只考虑血液系统对药物的排除，则有y(t+6)=y(t)/2。
5．4．4模型建立
由假设（1）可知，x(0)=1100，随着药物从胃肠道向血液系统的转移，胃肠道中药量x(t)会逐渐下降，其下降速度dxdt与x(t)成正比(比例系数λ＞0)，从而建立微分方程模型： 


dxdt=－λx，x(0)=1100(5．41)


由假设（2）可知，y(0)=0，从t=0开始，药物从胃肠道向血液系统转移(血液系统吸收药物)，同时血液系统启动对药物的排除，血液系统中的药量y(t)先增后减。y(t)因吸收作用而增长的速度为λx，因排除作用而减少的速度与y(t)成正比(比例系数μ＞0)，所以y(t)满足微分方程： 


dydt=λx－μy，y(0)=0(5．42)


5．4．5模型求解
求解式(5．41)和式(5．42)的MATLAB代码及结果如下： 




>> syms x(t) y(t) lambda mu % 定义符号函数及参数

>> equ = ［diff(x) == -lambda*x， diff(y) == lambda*x-mu*y］;% 定义微分方程

>> cond = ［x(0) == 1100，y(0) == 0］;% 定义初值条件

>> ［x(t)，y(t)］ = dsolve(equ，cond)% 求解微分方程

x(t)=

1100e-λt

y(t)=

1100λe-μtλ-μ-1100λe-λtλ-μ


由以上代码可知： 


x(t)=1100e－λt(5．43)

y(t)=1100λλ－μe－μt－e－λt(5．44)


接下来由氨茶碱吸收和排除的半衰期确定比例系数λ和μ。由于氨茶碱吸收的半衰期约为5h，可知x(5)=x(0)/2，即1100e－5λ=1100/2，解得λ=(ln2)/5≈0．1386。
为了根据氨茶碱排除的半衰期确定μ，不考虑血液系统对药物的吸收，只考虑血液系统对药物的排除，此时dydt=－μy，假设τ时刻血液系统中的药量y(τ)=a，解得y(t)=ae－μ(t－τ)，t≥τ。由于氨茶碱排除的半衰期约为6h，可知y(τ+6)=y(τ)/2，即ae－6μ=a/2，解得μ=(ln2)/6≈0．1155。
求解λ和μ的MATLAB代码如下： 



>> lambda = solve(x(5) == x(0)/2)% 求解血液对药物的吸收率系数

lambda = 

log（2）5

>> syms y2(t) mu tau a% 定义符号函数及参数

>> y2(t) = dsolve(diff(y2) == -mu*y2， y2(tau) == a)% 求解微分方程

y2(t) = 

ae-μteμτ


>> mu = solve(y2(tau+6) == y2(tau)/2，mu)% 求解血液对药物的排除率

mu = 

log（2）6


将λ和μ的值代入式(5．43)和式(5．44)可得： 


x(t)=1100e－(ln2)t/5(5．45)

y(t)=6600e－(ln2)t/6－e－(ln2)t/5(5．46)

5．4．6结果与分析
下面根据模型求解结果绘制x(t)和y(t)随时间变化的曲线(如图5．42所示)，并计算孩子就诊时(t=2)的血药含量，从而判断孩子的中毒状况。


图5．42胃肠道中的药量x(t)和血液系统中的药量y(t)


根据模型假设（4）和问题描述，孩子的血液总量为2000ml，当血药浓度达到100μg/ml，即血液中药物含量达到200mg时，会出现严重中毒，当血药浓度达到200μg/ml，即血液中药物含量达到400mg时则可致命。将t=2代入式(5．46)可得y(2)=236．5588mg，由此可知孩子就诊时(图5．42中A点)已经出现严重中毒。结合图5．42可知，若不及时进行施救，血液中药物含量将在一段时间后超过400mg，这将会危及孩子性命。
进一步地，还可根据模型式(5．46)计算血液中药物含量分别达到200mg和400mg的具体时刻，以及在不进行施救的情况下血液中药物含量达到峰值的时刻和相应的峰值。记血液中药物含量分别达到200mg、400mg和峰值的具体时刻分别为tB，tC和tD，则有


y(tB)=200，y(tC)=400，y′(tD)=0(5．47)


求解式(5．47)并绘图的MATLAB代码如下： 



>> y2h = double(subs(y，2))% t = 2h 时血液中药量

y2h =

236.5588

>> tB = double(solve(y(t) == 200))% 严重中毒时刻

tB =

1.6090

>> tC = double(solve(y(t) == 400))% 致命中毒时刻

tC =

4.8648

>> tD = double(solve(diff(y) == 0)) % 不加施救情形下，血药量达到峰值时刻

tD =

7.8910

>> ymax = double(subs(y，tD)) % 最大血药量

ymax =

442.0653



% 结果可视化

>> figure% 创建figure窗口

>> fplot(x，［0，25］，'b--')% 胃肠道中药量变化曲线

>> hold on% 开启图形保持

>> fplot(y，［0，25］，'r')% 血液中药量变化曲线

>> grid on% 添加主网格线

>> set(gca，'XMinorGrid'，'on'，'YMinorGrid'，'on')% 添加次网格线

>> xlabel('时间 t(h)'); ylabel('药量(mg)');% 添加坐标轴标签

>> plot(［2，tB，tC，tD］，［y2h，200，400，ymax］，'b*')% 绘制各时刻标记点

>> text(［2，tB，tC，tD］+0.5，［y2h，200，400，ymax］，...

{'A'，'B'，'C'，'D'})% 添加标记点的文本标注

>> legend('x(t)'，'y(t)')% 添加图例



由以上求解结果可知血液中药量y(t)在tB=1．6090h达到200mg，在tC=4．8648h(到达医院后2．8648h)达到400mg，在tD=7．8910h(到达医院后5．8910h)达到峰值442．0653mg。以上各状态分别对应图5．42中B、C、D点。
5．4．7施救方案
模型计算结果表明，若不及时进行施救，孩子会有生命危险，因此需要选择合适的施救方案。下面针对口服活性炭和体外血液透析两种施救方案，重新计算血液中药量的变化规律。
1． 口服活性炭
假设孩子到达医院时(t=2)即开始施救，施救中(t≥2)血液里的药量记为z(t)。由于口服活性炭可使血液中药物的排除率增至原来(人体自身)的2倍，则血液中药物的排除率与血液中药量的比例系数将变为原来的2倍，即μ=(ln2)/3，新的模型为


dxdt=－λx，x(0)=1100
dzdt=λx－μz，t≥2，z(2)=236．5588
λ=ln25，μ=ln23(5．48)



求解模型式(5．48)的MATLAB代码如下： 



% ---------------施救方案：口服活性炭-------------------

>> syms x(t) z(t)% 定义符号函数

>> lambda = log(2)/5;% 定义lambda

>> mu = 2*log(2)/6;% 定义mu

>> eq2 = ［diff(x) == -lambda*x，diff(z) == lambda*x-mu*z］; % 定义微分方程

>> cond2 = ［x(0) == 1100，z(2) == 236.5588］; % 定义初值条件

>> ［x(t)，z(t)］ = dsolve(eq2，cond2);% 求解微分方程

>> z = vpa(z，5) % 以5位有效数字形式显示

z(t) = 

1650.0*exp(-0.13863*t) - 1609.5*exp(-0.23105*t)



% 结果可视化

>> figure% 创建figure窗口

>> fplot(x，［0，25］，'b--')% 胃肠道中药量变化曲线

>> hold on% 开启图形保持

>> fplot(y，［0，25］，'r-.')% 血液中药量变化曲线(不施救)

>> fplot(z，［2，25］，'k')% 血液中药量变化曲线(施救)

>> grid on% 添加主网格线

>> set(gca，'XMinorGrid'，'on'，'YMinorGrid'，'on')% 添加次网格线

>> xlabel('时间 t(h)');ylabel('药量(mg)');% 添加坐标轴标签

>> legend('x(t)'，'y(t)'，'z(t)')% 添加图例



>> t4 = double(solve(diff(z) == 0))% 施救情形下，血药量达到峰值时刻

t4 =

5.2582

>> zmax = double(subs(z，t4))% 最大血药量

zmax =

318.3973



由以上代码可知： 



z(t)=1650e－0．13863t－1609．5e－0．23105t，t≥2(5．49)


x(t)，y(t)，z(t)随时间变化的曲线如图5．43所示。由z′(t)=0解得t=5．2582，相应的z(t)=318．3973。很显然，经口服活性炭施救后，血液中药量z(t)在t=5．2582h时达到峰值318．3973mg，远低于y(t)的峰值(442．0653mg)和致命水平(400mg)。


图5．43经口服活性炭施救后血液系统中的药量z(t)



2． 体外血液透析
由于体外血液透析可使血液中药物的排除率增至原来(人体自身)的6倍，则血液中药物的排除率与血液中药量的比例系数将变为原来的6倍，即μ=ln2，新的模型为： 


dxdt=－λx，x(0)=1100
dzdt=λx－μz，t≥2，z(2)=236．5588
λ=ln25，μ=ln2(5．410)


求解模型式(5．410)的MATLAB代码如下： 



% ---------------施救方案：体外血液透析-------------------

>> syms x(t) z(t)% 定义符号函数

>> lambda = log(2)/5;% 定义lambda

>> mu = 6*log(2)/6;% 定义mu

>> eq2 = ［diff(x) == -lambda*x，diff(z) == lambda*x-mu*z］; % 定义微分方程

>> cond2 = ［x(0) == 1100，z(2) == 236.5588］; % 定义初值条件

>> ［x(t)，z(t)］ = dsolve(eq2，cond2);% 求解微分方程

>> z = vpa(z，5) % 以5位有效数字形式显示

z(t)=

275.0*exp(-0.13863*t) + 112.59*exp(-0.69315*t)



% 结果可视化

>> figure% 创建figure窗口

>> fplot(x，［0，25］，'b--')% 胃肠道中药量变化曲线

>> hold on% 开启图形保持

>> fplot(y，［0，25］，'r-.')% 血液中药量变化曲线(不施救)





>> fplot(z，［2，25］，'k')% 血液中药量变化曲线(施救)

>> grid on% 添加主网格线

>> set(gca，'XMinorGrid'，'on'，'YMinorGrid'，'on')% 添加次网格线

>> xlabel('时间 t(h)');ylabel('药量(mg)');% 添加坐标轴标签

>> legend('x(t)'，'y(t)'，'z(t)')% 添加图例



由以上代码可知： 


z(t)=275e－0．13863t－112．59e－0．69315t，t≥2(5．411)


x(t)，y(t)，z(t)随时间变化的曲线如图5．44所示。很显然，经体外血液透析施救后，血液中药量z(t)关于时间t单调递减，说明体外血液透析的救治方案是十分有效的。


图5．44经体外血液透析施救后血液系统中的药量z(t)


由以上计算可知两种施救方案均能有效地将血液中药量降至致命水平以下，虽然体外血液透析可使血液中药量下降更快，但其治疗成本较高，并且存在安全风险，综合考虑，应该选择口服活性炭方案。