数学软件 Matlab —— Matlab 符号运算

主要内容 Matlab 符号运算介绍 符号对象与基本符号运算 symvar、 subs 和 vpa 常见的符号计算（重点内容）

符号运算 符号运算的特点 计算以推理方式进行，不受计算误差累积所带来的困扰 符号计算可以给出完全正确的封闭解，或任意精度的数值解（封闭解不存在时 ） 符号计算指令的调用比较简单，与教科书上的公式相近 符号计算所需的运行时间相对较长

Matlab 符号运算 Matlab 符号运算 Matlab 符号运算是通过符号数学工具箱（Symbolic Math Toolbox）来实现的。 Matlab 的符号数学工具箱可以完成几乎所有得符号运算功能，如：符号表达式的运算，符号矩阵的运算，符号微积分，符号作图，符号代数方程求解，符号微分方程求解等。 此外，该工具箱还支持可变精度运算，即支持以指定的精度返回结果。

符号运算举例 solve('a*x^2+b*x+c=0') x=sym('x'); diff(cos(x)^2) syms a b x; 求的根 f (x) = (cos x)2 的一次导数 x=sym('x'); diff(cos(x)^2) 计算 f (x) = x2 在区间 [a, b] 上的定积分 syms a b x; int(x^2,a,b)

内容提要 Matlab 符号运算介绍 符号对象与基本符号运算 symvar、 subs 和 vpa 常见的符号计算

符号对象 Matlab 符号对象 在进行符号运算时，必须先定义基本的符号对象，可以是 符号变量、符号表达式等 符号对象是一种数据结构 在进行符号运算时，必须先定义基本的符号对象，可以是 符号变量、符号表达式等 符号对象是一种数据结构 符号表达式：含有符号对象的表达式称 符号矩阵/数组：元素为符号表达式的矩阵/数组

符号对象的建立 符号对象的定义/声明：sym、syms a 是符号变量 b 是符号常量 C 是符号矩阵 例： 符号变量 = sym(x) 参数 x 可以是一个数或数值矩阵，也可以是字符串 例： a=sym('a') a 是符号变量 b=sym('1/3') b 是符号常量 C=sym('[1 ab; c d]') C 是符号矩阵

符号对象的建立 符号对象的定义/声明：sym、syms syms 符号变量1 符号变量2 ... 符号变量n 例： a=sym('a'); b=sym('b'); c=sym('c'); 例： syms a b c;

符号表达式 符号表达式：含符号对象的表达式 例： 建立符号表达式通常有以下 2 种方法： (1) 用 sym 函数直接建立符号表达式 (2) 使用已经定义的符号变量组成符号表达式 例： y=sym('sin(x)+cos(x)') x=sym('x'); y=sin(x)+cos(x) syms x; y=sin(x)+cos(x)

Matlab 符号运算采用的运算符和基本函数，在形状、名称和使用上，都与数值计算中的运算符和基本函数完全相同 基本符号运算 Matlab 符号运算采用的运算符和基本函数，在形状、名称和使用上，都与数值计算中的运算符和基本函数完全相同 基本运算 普通运算： + - * \ / ^ 数组运算： .* .\ ./ .^ 矩阵转置： ' .' 基本数学函数 三角函数与反三角函数、指数函数、对数函数等 sin,cos,asin,acos,exp,log,abs,diag,tril,triu, ...

符号矩阵 符号矩阵的生成 符号矩阵中元素的引用和修改 使用 sym 函数直接生成 A=sym('[1+x, sin(x); 5, exp(x)]') 将数值矩阵转化成符号矩阵 B=[2/3, sqrt(2); 5.2, log(3)]; C=sym(B) 符号矩阵中元素的引用和修改 A=sym('[1+x, sin(x); 5, exp(x)]'); A(1,2) % 引用 A(2,2)=sym('cos(x)') % 重新赋值

内容提要 Matlab 符号运算介绍 符号对象与基本符号运算 symvar、 subs 和 vpa 常见的符号计算

symvar 列出符号表达式中的符号变量 symvar(s) symvar(s, N) 例： 按字母顺序列出符号表达式 s 中的所有符号变量 列出符号表达式 s 中离 x 最近的 N 个符号变量 若有两个符号变量与 x 的距离相等，则ASCII 码大者优先 常量 pi, i, j 不作为符号变量 例： f=sym('2*v-3*y+z^2+5*a') symvar(f) symvar(f,2)

subs 符号替换 subs(s,x,a) 例： 用给定的数据替换符号表达式中的指定的符号变量 用 a 替换符号表达式 s 中的符号变量 x Matlab 演示 syms x y u v; f1=2*x+y-1; f2=subs(f1,x,u) f3=subs(f1,y,2+3) f3=subs(f1,{x,y},{u,v})

subs 举例 例：指出下面各条语句的输出结果 f=sym('2*u'); f1=subs(f,'u',2) f2=subs(f,'u','u+2') f3=subs(f,'u',[1,2]) a=3; f4=subs(f2,'u',a+2) f5=subs(f2,'u','a+2') syms x y; f6=subs(f,'u',x+y) f7=subs(f6,{x,y},{1,2}) f8=subs(f6,{x,y},{x+y,x+y}) f=2*u f1=4 f2=2*u+4 f3=[2,4] f4=14 f5=2*a+8 f6=2*x+2*y f7=6 f8=4*x+4*y Matlab 演示 符号变量若没有声明，则需要加上单引号！

vpa 可变精度 vpa(s,n) 例： 计算表达式 s 的值，保留 n 位有效数字 返回值是符号对象 x1=vpa(sin(pi/2),10) x2=vpa(pi^3,3) x3=vpa(pi,100) Matlab 演示

内容提要 Matlab 符号运算介绍 符号对象与基本符号运算 symvar、 subs 和 vpa 常见的符号计算（六类运算） 因式分解、展开、合并、简化、通分和反函数等 计算极限 计算导数 计算积分 符号级数求和 代数方程和微分方程的求解（重点与难点）

因式分解 因式分解 factor(f) 例： factor 也可用于正整数的分解 例： syms x; f=x^6+1; factor(f) s=factor(100) 对大整数进行因式分解时可以先将其转化成符号常量 factor(12345678901234567890) % ERROR factor(sym('12345678901234567890'))

函数展开 函数展开 expand(f) 例： 例： 多项式展开 syms x; f=(x+1)^6; expand(f) 三角函数展开 syms x y; f=sin(x+y); expand(f)

合并同类项 合并同类项 collect(f,v) % 按指定变量 v 进行合并 collect(f) % 按默认变量进行合并 例： 默认变量：symvar(f) 的返回结果 例： syms x y; f=x^2*y+y*x+y^2+2*x ; collect(f) collect(f,y) syms u v; g=u^2*v+u*v^3-u^2+v; collect(g)

函数简化 函数简化 y=simplify(f) 例： 对符号表达式 f 进行简化 syms x; f=sin(x)^2 + cos(x)^2; y=simplify(f)

函数简化 函数简化 y=simple(f) 例：化简 对 f 尝试多种不同的方法（包括 simplify）进行简化， 以寻求其最简短形式 syms x; f=(cos(x)^2-sin(x)^2)*sin(2*x)*(exp(2*x) ... -2*exp(x)+1)/(exp(2*x)-1); y1=simplify(f) y2=simple(f)

函数简化 通分 [N,D]=numden(f) 例： N 为通分后的分子，D 为通分后的分母 syms x y; f=x/y+y/x; [n,d]=numden(sym(112/1024))

horner 多项式 horner 多项式：嵌套形式的多项式 例： syms x; f=x^4+2*x^3+4*x^2+x+1; g=horner(f)

求反函数 反函数 finverse(f,v) % 求 f 关于指定变量 v 的反函数 finverse(f) % 求 f 关于默认变量的反函数 例：计算函数 的反函数 syms x t; f=x^2+2*t; g1=finverse(f,x) g2=finverse(f,t)

计算极限 计算极限 limit(f,x,a) % 计算 limit(f,a) % 当默认变量趋向于 a 时的极限 limit(f,x,a,'right') % 计算右极限 limit(f,x,a,'left') % 计算左极限 例：计算 ， syms x h n; L=limit((log(x+h)-log(x))/h,h,0) M=limit((1-x/n)^n,n,inf)

计算导数 计算导数 g=diff(f,v) % 求符号表达式 f 关于变量 v 的导数 g=diff(f) % 计算关于默认变量的导数 g=diff(f,v,n) % 求 f 关于 v 的 n 阶导数 例： syms x; f=sin(x)+3*x^2; g1=diff(f,x) g2=diff(f,x,3)

计算积分 计算积分 int(f,v,a,b) % 计算定积分 int(f,a,b) % 计算关于默认变量的定积分 例：计算 和 syms x; f=(x^2+1)/(x^2-2*x+2)^2; I=int(f,x) K=int(exp(-x^2),x,0,inf)

符号级数求和 符号级数求和 symsum(f,v,a,b) % 级数求和 symsum(f,a,b) % 关于默认变量求和 syms n; f=1/n^2; S=symsum(f,n,1,inf) S100=symsum(f,n,1,100) 例：计算函数级数 syms n x; f=x/n^2; S=symsum(f,n,1,inf)

代数方程求解 代数方程求解 solve(f,v) % 求方程关于指定自变量的解 例： 这里 f 可以用字符串表示或符号表达式 得不到解析解时，给出数值解 例： solve('2*x-3') % 或 solve('2*x-3=0') syms x; solve(2*x-3) % 不能写成 solve(2*x-3=0) syms x; solve(2*x-sin(x)+1)

微分方程求解 ode23s、ode23tb 用 Maltab自带函数 解初值问题 求微分方程解析解：dsolve 求微分方程数值解*（自学，选学）： ode45、ode23、 ode113、ode23t、ode15s、 ode23s、ode23tb

dsolve 求解析解 dsolve 的使用 例 1：求微分方程 的通解，并验证。 y=dsolve('eq1','eq2', ... ,'cond1','cond2', ... ,'v') 其中 y 为输出， eq1、eq2、...为微分方程，cond1、cond2、...为初值条件，v 为自变量。 例 1：求微分方程 的通解，并验证。 y=dsolve('Dy+2*x*y=x*exp(-x^2)','x') syms x; diff(y)+2*x*y - x*exp(-x^2)

dsolve 的使用 几点说明 微分方程中用 D 表示对 自变量 的导数，如： Dy y'； D2y y''； D3y y''' 如果省略初值条件，则表示求通解； 如果省略自变量，则默认自变量为 t dsolve('Dy=2*x','x'); ％ dy/dx = 2x dsolve('Dy=2*x'); ％ dy/dt = 2x 若找不到解析解，则返回其积分形式。

dsolve 的使用 例 2：求微分方程 在初值条件 下的特解，并画出解函数的图形。 例 2：求微分方程 在初值条件 下的特解，并画出解函数的图形。 y=dsolve('x*Dy+y-exp(x)=0','y(1)=2*exp(1)','x') ezplot(y);

dsolve 的使用 例3：求微分方程组 在初值条件 下的特解，并画出解函数的图形。 例3：求微分方程组 在初值条件 下的特解，并画出解函数的图形。 [x,y]=dsolve('Dx+5*x+y=exp(t)','Dy-x-3*y=0', ... 'x(0)=1', 'y(0)=0', 't') ezplot(x,y,[0,1.3]); 注：解微分方程组时，如果所给的输出个数与方程个数相同，则方程组的解按词典顺序输出；如果只给一个输出，则输出的是一个包含解的结构（structure）类型的数据。

dsolve 的使用 例： [x,y]=dsolve('Dx+5*x=0','Dy-3*y=0', ... 'x(0)=1', 'y(0)=1','t') r = dsolve('Dx+5*x=0','Dy-3*y=0', ... 'x(0)=1', 'y(0)=1','t') 这里返回的 r 是一个 结构类型 的数据 r.x %查看解函数 x(t) r.y %查看解函数 y(t) dsolve的输出个数只能为一个 或 与方程个数相等 只有很少一部分微分方程（组）能求出解析解。 大部分微分方程（组）只能利用数值方法求数值解。

上机作业 1. 教材P27的习题4，写入文件m04_1.m； 2. 求解下列微分方程（组），写入文件m04_2.m，并且对于初值问题还要求画出解函数的图形；(5)(6)选做； Matlab 演示 3.（选做）教材P28的习题6~10，写入文件m04_3.m。

上机要求 上机要求强调 将完成每题所用的命令写入一个规定文件名的文件中 然后将这些文件作为附件，通过 foxmail 以邮件形式发给 admin@system.mail 邮件主题为：机号-学号-姓名 其中机号为 两位数 三个字段之间用英文状态下的减号链接

书面作业 教材P35-38有4段关于Matlab用法的英文材料，试将其中的(Topic 3) 段（有关 Matlab 图形）翻译为中文. 【具体要求】将译文的手写版写在作业本上，欢迎给出进行中文、英文对照的作业！ 40

求微分方程数值解*（自学，选学） [T,Y] = solver(odefun,tspan,y0) 其中 y0 为初值条件，tspan为求解区间；Matlab在数值求解时自动对求解区间进行分割，T (列向量) 中返回的是分割点的值(自变量)，Y (数组) 中返回的是这些分割点上的近似解，其列数等于因变量的个数。 solver 为Matlab的ODE求解器（可以是 ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb） 没有一种算法可以有效地解决所有的 ODE 问题，因此MATLAB 提供了多种ODE求解器，对于不同的ODE，可以调用不同的求解器。

Matlab的ODE求解器*（自学，选学） 特点 说明 ode45 非刚性 单步法；4，5 阶 R-K 方法；累计截断误差为 (△x)3 大部分场合的首选方法 ode23 单步法；2，3 阶 R-K 方法；累计截断误差为 (△x)3 使用于精度较低的情形 ode113 多步法；Adams算法；高低精度均可到 10-3～10-6 计算时间比 ode45 短 ode23t 适度刚性 采用梯形算法 适度刚性情形 ode15s 刚性 多步法；Gear’s 反向数值微分；精度中等 若 ode45 失效时，可尝试使用 ode23s 单步法；2 阶Rosebrock 算法；低精度 当精度较低时，计算时间比 ode15s 短 ode23tb 梯形算法；低精度 当精度较低时，计算时间比ode15s短

参数说明*（自学，选学） [T,Y] = solver(odefun,tspan,y0) 例 ： odefun 为显式常微分方程，可以用命令 inline 定义，或在函数文件中定义，然后通过函数句柄调用。 求初值问题 的数值解，求解范围为 [0,0.5] 例 ： fun=inline('-2*y+2*x^2+2*x','x','y'); [x,y]=ode23(fun,[0,0.5],1); 注：也可以在 tspan 中指定对求解区间的分割，如： [x,y]=ode23(fun,[0:0.1:0.5],1); %此时 x=[0:0.1:0.5]

数值求解举例*（自学，选学） 如果需求解的问题是高阶常微分方程，则需将其化为一阶常微分方程组，此时必须用函数文件来定义该常微分方程组。 求解 Ver der Pol 初值问题 例 ： 令 ，则原方程可化为

数值求解举例*（自学，选学） 先编写函数文件 verderpol.m 再编写脚本文件 vdpl.m，在命令窗口直接运行该文件。 function xprime=verderpol(t,x) global mu; xprime=[x(2); mu*(1-x(1)^2)*x(2) - x(1)]; 再编写脚本文件 vdpl.m，在命令窗口直接运行该文件。 clear; global mu; mu=7; y0=[1;0]; [t,x]=ode45('verderpol',[0,40],y0); plot(t,x(:,1),'r-', t,x(:,2),'b-');