Matlab符号计算(含作业).docx

第 2 章 符号计算符号计算：解算数学表达式、方程不是在离散化的数值点上进行，而是凭借一系列恒等式，数学定理，通过推理和演绎，获得解析结果。符号计算建立在数值完全准确表达和推演严格解析的基础之上，所得结果完全准确。特点：一相对于MATLAB的数值计算“引擎”和“函数库”而言，符号计算的“引擎”和“函数库”是独立的。二在相当一些场合，符号计算解算问题的命令和过程，显得比数值计算更自然、更简明。三大多数理工科的本科学生在学过高等数学和其他专业基础课以后，比较习惯符号计算的解题理念和模式。 2.1 符号对象和符号表达式 MATLAB依靠基本符号对象（包括数字、参数、变量）、运算符及一些预定义函数来构造和衍生符号表达式和符号方程。2.1.1 基本符号对象和运算算符1. 生成符号对象的基本规则l 任何基本符号对象（数字、参数、变量、表达式、函数）都必须借助专门的符号命令sym、syms、symfun定义。l 任何包含符号对象的表达式或方程，将继承符号对象的属性。2. 精准符号数字和符号常数符号（类）数字的定义：sym(Num)采用精准数值类数创建精准的符号数字（推荐格式!）sc=sym(Num)采用精准数值类数创建精准的符号常数sc（推荐格式!）说明：若输入量Num是精准的浮点数（如0.321、10/3等），能生成精准的符号数字； 若输入量Num是诸如sin(0.3)的数值表达式，那么就只能生成由数字表达式获得的16位精度的近似符号数字。sym(Num) 采用有理分数字符串创建精准的符号数字sc=sym(Num) 采用有理分数字符串创建精准的符号常数sc说明： Num必须处于（英文状态下的）单引号内，构成字符串（关于字符串参见附录A）；只有当字符串数字Num采用诸如321/1000、10/3等整数构成的有理分数形式表达时，sym(Num) 才能生成精准的符号数字；若字符串数字用诸如0.321、3.21e-1等“普通小数或科学记述数”表达，那么只能产生“近似符号数字”。在默认情况下，该近似符号数字为32位精度。【例2.1-1】（1）创建完全精准的符号数字或数字表达式clear allR1=sin(sym(0.3)% 输入量为普通小数R2=sin(sym(3e-1)% 输入量为科学记述数R3=sin(sym(3/10)% 输入量为有理分数R4=sin(sym(3/10)% 输入量为“整数构成的有理分数”字符串数字disp(R1属于什么类别？ 答：,class(R1)disp(R1与R4是否相等？（是为1，否为0） 答：,int2str(logical(R1=R4)R1 =sin(3/10)R2 =sin(3/10)R3 =sin(3/10)R4 =sin(3/10)R1属于什么类别？ 答：symR1与R4是否相等？（是为1，否为0） 答：1 （2）产生具有32位精度的“近似”符号数字（杜绝使用！）S1=sin(sym(0.3)% sym的输入量是字符串小数，生成32位精度下的% 近似符号数，进而在sin作用下给出近似符号数。S2=sin(sym(3e-1)% syms的输入量是字符串科学记述数。eRS=vpa(abs(R1-S1),64);disp(S1属于什么类别？ 答：,class(S1)disp(S1与R1是否相同？ 答: ,int2str(logical(R1=S1)disp(S1与R1的误差为)disp(double(eRS)S1 =0.29552020666133957510532074568503S2 =0.29552020666133957510532074568503S1属于什么类别？ 答：symS1与R1是否相同？ 答: 0S1与R1的误差为 6.3494e-41（3）产生具有16位精度的“近似”符号数字（杜绝使用！）F1=sym(sin(3/10)% sym的输入量为双精度表达式sin(3/10)，% 就只能创建出仅16位精度的近似符号数。F2=sym(sin(0.3)% 同上eFS=vpa(abs(F1-S1),32);disp(F1属于什么类别？ 答：,class(F1)disp(S1与F1是否相同？ 答: ,int2str(logical(F1=S1)disp(F1与S1的误差为)disp(double(eFS) F1 =5323618770401843/18014398509481984F2 =5323618770401843/18014398509481984F1属于什么类别？ 答：symS1与F1是否相同？ 答: 0F1与S1的误差为 2.8922e-17 3. 基本符号变量经典教科书里，表达式e-axsinbx中的a,b称为参数，x为变量。在MATLAB的符号计算中，a、b、x统称为基本符号变量，其中，x总被默认为“待解（自由）符号变量”，其他被作为“符号参数”处理。定义基本符号变量的命令格式：para=sym( para) 定义单个复数域符号变量parapara=sym( para, Flag) 定义单个Flag指定域符号变量parasyms para 定义单个复数域符号变量para的另一种方式syms para Flag 定义单个Flag指定域符号变量para的另一种方式syms para1 para2 paraN 定义多个复数域符号变量para1 para2 paraNsyms para1 para2 paraN Flag 定义多个Flag指定域符号变量para1 para2 paraNl 符号参数名不要用处于“字母表中小写字母x及其两侧的英文字母”开头。l Flag表示数域的限定性假设，可具体取以下词条：positive 正实数域； real 实数域。l 默认值是复数域符号变量l sym命令只能对单变量作用，syms不能用于对数值、常数相关的定义。4. 符号计算中的各种算符l 与数值计算中的算符在形状、名称和使用方法上相同。2.1.2 符号计算中的函数命令表2.1-2 MATLAB中可调用的符号计算函数指令类别情 况 描 述与数值计算对应关系符号工具包函数三角函数、双曲函数及反函数；除atan2外相同指数、对数函数（如exp, expm）相同复数函数（如abs, angle）相同矩阵分解函数（如eig, svd）相同方程求解函数solve不相同微积分函数（如diff, int）不完全相同积分变换和反变换函数（如laplace, ilaplace）只有离散Fourier变换绘图函数（如ezplot, ezsurf）数值绘图命令更丰富特殊函数单位脉冲和阶跃函数（如dirac, heaviside） 不完全对应函数（如beta, gamma）椭圆积分（如ellipke）贝塞尔函数（如besseli, besselj）MuPAD库函数借助evalin和feval指令可调用比mfunlist所列范围更广泛的MuPAD库函数；需要具备MuPAD语言知识。无对应函数说明2.1.3 符号表达式和符号函数1. 符号表达式和符号函数（1） 为表达某种数学算式、实现某种计算目的，采用基本符号对象（数字、常数、变量）、运算符、MATLAB函数命令等基本要素编写而成的M码。（2）为表达变量间抽象（或具体）约束关系而编写的M码。在符号函数中，构成函数关系的变量名必须明确指定。即，在定义符号函数时，不仅要指定函数名，而且要指定变量名。比如syms f（x,y）就定义了一个以x、y为变量的抽象符号函数f。2 自由符号变量 解题通常是围绕自由符号变量进行。 确定自由符号变量的规则：l 在专门指定变量名的符号运算中，解题一定围绕指定变量名进行。l x是首选自由符号变量，其后的次序规则是：与x的ASCII码值之差的绝对值小的字母优先；差绝对值相同时，ASCII码值大的字母优先。l 自动识别符号变量时，字母的优先次序为x,y,w,z,v等，大写后排。自动识别表达式中自由、独立的符号变量的命令： symvar(expression) 列出表达式中的所有基本符号变量 symvar(expression, n) 列出表达式中n个认定的自由变量【例2.1-2】1）各种符号对象的创建clearsyms a b c x y u v% 定义基本符号对象syms F(X,Y,Z)% 定义“抽象”符号函数k=sym(3)% 定义符号常数G=sym(p*sqrt(q)+r*sin(t)% 创建符号表达式EXPR=a*G*u+(b*x2+k)*v% 创建“衍生”符号表达式f(x,y)=a*x2+b*y2-c% 创建“具体”符号函数disp(F)% 显示抽象符号函数 k =3G =p*q(1/2) + r*sin(t)EXPR =v*(b*x2 + 3) + a*u*(p*q(1/2) + r*sin(t)f(x, y) =a*x2 + b*y2 - cF(X, Y, Z)symbolic function inputs: X, Y, Z 2）symvar对EXPR符号表达式的作用symvar(EXPR)%ans = a, b, p, q, r, t, u, v, x symvar(EXPR,20) %ans = x, v, u, t, r, q, p, b, a symvar(EXPR,1)% ans =x 3）symvar对符号函数的作用disp(symvar(f)% a, b, c, x, y disp(symvar(f,2)% x, y 【例2.1-3】用符号法求方程的解。1）产生符号表达式和符号函数clearsyms u v w z% Eq=u*w2+z*w-v% 表达式g(z)=u*w2+z*w=v% 函数Eq =u*w2 + z*w - vg(z) =u*w2 + z*w = v 2）symvar认定的自由变量symvar(Eq,1)%ans =w symvar(g(z),1)% ans =w 3）solve对默认自由变量解方程R1=solve(Eq)%关于w解方程u*w2+z*w-v=0 R2=solve(g)%关于w解g（z）所表达的方程 R1 = -(z + (z2 + 4*u*v)(1/2)/(2*u) -(z - (z2 + 4*u*v)(1/2)/(2*u)R2 = -(z + (z2 + 4*u*v)(1/2)/(2*u) -(z - (z2 + 4*u*v)(1/2)/(2*u) 4）对变量z求解S1=solve(Eq,z)%S2=solve(g(z),z)% S1 =(- u*w2 + v)/wS2 =(- u*w2 + v)/w 5）检验求解结果的正确性disp(simplify(subs(Eq,z,S1)% S1代替z，观察Eq是否为00 disp(simplify(g(S2)% S2代替z，观察g(S2)方程式是否成立TRUE 说明不要把g(z)理解为以z为自由变量的符号函数。【例2.1-4】symvar确定自由变量是对整个矩阵进行的。syms a b t u v x yA=a+b*x,sin(t)+u;x*exp(-t),log(y)+v%symvar(A,1)% A = a + b*x, u + sin(t) x*exp(-t), v + log(y)ans =x 2.1.4 符号对象的识别为了函数命令与数据对象的适配，MATLAB提供了用于识别数据对象属性的命令：class(var) 给出变量var的数据类别（如double，sym等）isa(var, Obj) 若变量var是Obj代表的类型，给出1，表示“真”whos 给出所有MATLAB内存变量的属性【例2.1-5】数据对象及其识别命令的使用。（1）cleara=1;b=2;c=3;d=4; % 产生4个数值变量Mn=a,b;c,d % 利用已赋值变量构成数值矩阵Mc=a,b;c,d % 字符串中的a,b,c,d与前面输入的数值变量无关Ms=sym(Mc) % Ms是一个符号矩阵，它与前面各变量无关Mn = 1 2 3 4Mc =a,b;c,dMs = a, b c, d （2）SizeMn=size(Mn)SizeMc=size(Mc)SizeMs=size(Ms) SizeMn = 2 2SizeMc = 1 9SizeMs = 2 2 （3）CMn=class(Mn)CMc=class(Mc)CMs=class(Ms) CMn =doubleCMc =charCMs =sym （4）isa(Mn,double)isa(Mc,char)isa(Ms,sym) ans = 1ans = 1ans = 1 （5）whos Mn Mc Ms Name Size Bytes Class Attributes Mc 1x9 18 char Mn 2x2 32 double Ms 2x2 60 sym 2.1.5 符号运算机理和变量假设1. 符号运算的工作机理符号计算是由MuPAD引擎在其专有的内存工作空间中执行，而只是把计算结果送回到MATLAB的内存空间。每当借助sym或syms命令定义一个带限定性假设的符号变量时，就发生以下过程：l 启动MuPAD引擎，并开启一个专供MuPAD使用的内存空间l 被定义变量保存至Matlab内存空间l 对变量的限定性假设被保存在MuPAD内存空间中，并对此后的MuPAD的工作方式进行约束。2. 对符号变量的限定性假设assume（assumption） 清空后由assumption定义的新假设assume（expr,set） 清空后进行“符号表达式expr属于集合set”的新假设设置assumeAlso(assumption) 继续追加由assumption定义的新假设assumeAlso(expr,set) 继续追加由“符号表达式expr属于集合set”的新假设a=sym(a,res) 创建带res限定性假设的符号变量asyms a res 同上说明：l 在不对符号变量进行专门设置的情况下，MuPAD符号计算总把变量默认为“复数变量”。l assumption是那些可以用符号表达式、符号方程、符号关系、符号逻辑描述的假设。l set可取字符串integer、rational、real中的任何一个，分别表示整数集、有理数集、实数集。l res只能取字符串positive或real中的任何一个，表示创建变量是实数或正数。3. 清除变量和撤销假设由于符号变量和其假设存放在不同的内存空间，因此删除符号变量和撤销关于变量的假设是两件需要分别处理的事。具体执行命令如下：clear x % 清除MATLAB内存中的x变量syms x clear % 撤销MuPAD内存中对变量x的任何假设，而恢复为“复数”变量syms（x,clear） % 功能同上assumptions(x) % 显示符号变量x的限定性假设assumptions % 显示MuPUD内存中已带限定性假设的全部符号变量reset(symengine) % 重启MuPUD引擎，清空MuPUD内存中所有内容【例2.1-6】符号变量的默认数域是复数域。1）在默认的复数域求根clear all% 清空MATLAB内存，清除MuPAD中的所有假设syms x%f=x3+475*x/100+5/2;%r=solve(f,x)% 求使f=0的全部根r = -1/2 (79(1/2)*i)/4 + 1/4 1/4 - (79(1/2)*i)/4 assumptions(x)% 获悉MuPAD内存中关于x的假设ans =Empty sym: 1-by-0 2）求实数根assume(x,real)% 限定x为实数的方法一r21=solve(f,x)% r21 =-1/2 syms x clear%assume(imag(x)=0)% 限定x为实数的方法二r22=solve(f,x) r22 =-1/2 disp(assumptions(x)% imag(x) = 0 3）第一、四象限根sym(x,clear)%assume(real(x)0)%r3=solve(f,x) ans =xr3 = (79(1/2)*i)/4 + 1/4 1/4 - (79(1/2)*i)/4 disp(assumptions(x)% 0 0)%r4=solve(f,x) r4 =(79(1/2)*i)/4 + 1/4 disp(assumptions(x) 0 imag(x), 0 real(x) 2.2 符号数字及表达式的操作2.2.1 符号数字转换成双精度数字Nd=double(Num_sym) 把符号数字Num_sym转换为双精度数字Nd说明：l 一般情况下，Nd是符号数字Num_sym双精度近似。l 2.1.1-2节中，已经介绍了sym(Num)能把数字类数字转换成符号数字l 注意：命令double(Num)是把字符串数字Num转换为各数字字符的ASCII码值数组。2.2.2 符号数字的任意精度表达形式为了兼顾计算精度和速度，或使某些无法用“封闭解析式”表达的计算结果以简洁的“任意精度符号数”表达。MATLAB提供了控制符号数字或表达式数字精度的命令：digits 显示当前环境下十进制符号数字的有效位数digits(n) 把十进制符号数字有效位数设定为nxs=vpa(x) 据表达式x得到digits指定精度下的符号数字xsxs=vpa(x,n) 据表达式x得到n位有效数字的符号数字xs说明： MATLAB对digits命令的默认精度设置是32位。 vpa(x,n)只在运行的当时起作用。【例2.2-1】digits, vpa命令的使用。（1）重新启动符号计算引擎，产生准确符号数字reset(symengine) % 重新启动符号计算引擎 sa=sym(1/3+sqrt(2) % 定义准确符号数字表达式sa =2(1/2) + 1/3 （2）变精度算法的计算结果，有效位数的含义digits% 观察当前有效位数 Digits = 32 format longa=1/3+sqrt(2)% 定义双精度数sa_Plus_a=vpa(sa+a,20)% 给出20位有效数字结果sa_Minus_a=vpa(sa-a,20)% a = 1.747546895706428sa_Plus_a =3.4950937914128567869sa_Minus_a =-0.000000000000000022658064826339973669 （3）digits设置和vpa指定对“数位”的不同影响sa32=vpa(sa)% 采用默认设置的32位有效数字 digits(48)% 设置48位有效数字sa5=vpa(sa,5)% 仅影响sa5数位，对其后无影响。sa48=vpa(sa)% 仍为48位有效数字 sa32 =1.747546895706428382135022057543sa5 =1.7475sa48 =1.74754689570642838213502205754303141190300520871 2.2.3 符号表达式的基本操作 collect（合并同类项） factor（进行因式或因子分解） numden（提取公因式）等最常用：simplify(EXPR) 对EXPR（符号表达式或矩阵）运用多种方法进行一轮简化simplify(EXPR,Steps,value,IgnoreAnalyticConstraints,true) 多轮纯粹表达形式简化【例2.2-2】简化。syms xf=(1/x3+6/x2+12/x+8)(1/3)g1=simplify(f)f =(12/x + 6/x2 + 1/x3 + 8)(1/3)g1 =(2*x + 1)3/x3)(1/3) g2=simplify(f,Steps,10,IgnoreAnalyticConstraints, true)% g2 =1/x + 2 2.2.4 表达式中的置换操作1. 公因子法简化表达RS=subexpr(S) 从S中自动提取公因子sigma，并把采用sigma重写的S赋给RSRS=subexpr(S,w) 从S中自动提取公因子，记为w，并把采用w重写的S赋给RSRS,w=subexpr(S,w) 该调用格式的效果与RS=subexpr(S, w)相同说明：S：符号表达式、符号表达式矩阵【例2.2-3】对符号矩阵进行特征向量分解。（）clear% 清空所有内存变量A=sym(a b;c d)% 经字符串直接定义符号矩阵V,D=eig(A)% 符号矩阵的特征值、特征向量分解A = a, b c, dV =(a/2+d/2-a2-2*a*d+d2+4*b*c)(1/2)/2)/c-d/c, (a/2+d/2+(a2-2*a*d+ d2+4*b*c)(1/2)/2)/c-d/c 1, 1D =a/2 + d/2 - (a2 - 2*a*d + d2 + 4*b*c)(1/2)/2, 0 0, a/2 + d/2 + (a2 - 2*a*d + d2 + 4*b*c)(1/2)/2 （2）subexpr(V;D)% 自动提取公因式who % 列出工作内存中的变量sigma = (a2 - 2*a*d + d2 + 4*b*c)(1/2)ans = (a/2 + d/2 - sigma/2)/c - d/c, (a/2 + d/2 + sigma/2)/c - d/c 1, 1 a/2 + d/2 - sigma/2, 0 0, a/2 + d/2 + sigma/2Your variables are:A D V ans sigma （3）Dw=subexpr(D,w) % 把自动提取的公因式记为w,Dw是用w重记D后的表达w = (a2 - 2*a*d + d2 + 4*b*c)(1/2)Dw = a/2 + d/2 - w/2, 0 0, a/2 + d/2 + w/2 （4）RVD,w=subexpr(V;D,w)% % 给出合成矩阵V;D的公因式表达方式RVD = (a/2 + d/2 - w/2)/c - d/c, (a/2 + d/2 + w/2)/c - d/c 1, 1 a/2 + d/2 - w/2, 0 0, a/2 + d/2 + w/2w =(a2 - 2*a*d + d2 + 4*b*c)(1/2) 2. 通用置换命令RES=subs(ES,old,new) 用new置换ES中的old后产生符号结果RESRES=subs(ES,new) 用new置换ES中的自由变量后产生符号结果RES【例2.2-4】用简单算例演示subs的置换规则。1）产生符号函数clearsyms a b x;f1=a*sin(x)+b f1 =b + a*sin(x) 2）被字符串置换f2=subs(f1,sin(x),log(y)%class(f2)% f2 =b + a*log(y)ans =sym 3）单个符号参数被置换f3=subs(f1,a,sym(3.11)%class(f3)% f3 =b + (311*sin(x)/100ans =sym 4）单个符号变量被数组置换f4=subs(f1,x,0,pi/2,pi)%class(f4) f4 = b, a + b, bans =sym 5）所有变量被置换format% 恢复双精度数字显示的默认设置format compact% 在Notebook中紧凑显示t=0:pi/10:2*pi;% （1*21）双精度数组f5=subs(f1,a,b,x,2,3,t);% 置换得到符号数字数组class(f5)plot(t,f5,r:,LineWidth,5)% ans =sym图2.2-1 利用符号表达式变量置换产生的单根曲线6）两次置换k=0.6;0.8;1; %（3*1）数组f6=subs(subs(f1,a,b,k,2),x,t);% class(f6)plot(t,f6) % ans =sym图2.2-2 利用两次subs置换产生的多根曲线2.3 符号微积分2.3.1 极限和导数的符号计算大学本科高等数学中的大多数微积分问题，都能用符号计算解决，手工笔算演绎的烦劳都可以由计算机完成。limit(f,v,a) 求极限 limit(f,v,a,right) 求右极限 limit(f,v,a,left) 求左极限 【例2.3-1】两种重要极限和。syms t x ks=sin(k*t)/(k*t);f=(1-1/x)(k*x);Lsk=limit(s,0)% t趋于0Ls1=subs(Lsk,k,1)%Lf=limit(f,x,inf)%Lf1=vpa(subs(Lf,k,sym(-1),48)%给出48位精度的自然数 Lsk =1Ls1 =1Lf =exp(-k)Lf1 =2.7182818284590452353602874713526624977572470937 diff(f,v,n) 求 (n缺省时，默认n=1)【例2.3-2】已知，求、 、。syms a t x;f=a,t3;t*cos(x), log(x);df=diff(f)%dfdt2=diff(f,t,2)%dfdxdt=diff(diff(f,x),t)% df = 0, 0 -t*sin(x), 1/xdfdt2 = 0, 6*t 0, 0dfdxdt = 0, 0 -sin(x), 0 2.3.2 序列/