500 likes | 707 Views
第 2 章 MATLAB 符号计算 2.1 符号对象和符号表达式 2.2 符号数字及表达式的操作 2.3 符号微积分 2.4 符号微分方程的求解 2.5 符号变换和符号卷积 2.6 符号代数方程的求解. 什么是符号计算?. 符号计算是指: 解算数学表达式、方程 不是在 离散化 的数值点上进行 ,而是凭借一系列恒等式,数学定理,通过推理和演绎,获得解析结果。. 符号计算的特点: 1 )符号表达式计算前必须 先定义符号变量 。 2 )符号计算是 精确 计算。 3 )符号计算的计算 速度较慢 。
E N D
第2章 MATLAB符号计算 2.1 符号对象和符号表达式 2.2 符号数字及表达式的操作 2.3 符号微积分 2.4 符号微分方程的求解 2.5 符号变换和符号卷积 2.6 符号代数方程的求解
什么是符号计算? • 符号计算是指: 解算数学表达式、方程不是在离散化的数值点上进行,而是凭借一系列恒等式,数学定理,通过推理和演绎,获得解析结果。
符号计算的特点: 1)符号表达式计算前必须先定义符号变量。 2)符号计算是精确计算。 3)符号计算的计算速度较慢。 4)符号计算的运算符和函数与数值计算中的运算符和函数几乎完全相同
2.1 符号对象和符号表达式 • 符号对象的创建 (1)sym(‘变量’,参数) 功能:把变量定义为符号对象。 参数有三种选择: ’positive’ 表示为“正”符号变量, ’real’ 表示为“实”符号变量, ’unreal’ 表示为“非实”符号变量。 例如:S=sym('88') sa=sym('pi+sqrt(5)')
(2)syms函数 格式:syms arg1 arg2 …参数 功能:创建多个符号变量。 例如: syms r_l real syms a b c d positive
一、符号数字与数值数字 【例2.1-1】符号(类)数字与数值(类)数字之间的差异。 a=pi+sqrt(5) sa=sym('pi+sqrt(5)') Ca=class(a) Csa=class(sa) vpa(sa-a) 结果: a =5.3777 sa =pi+sqrt(5) Ca =double Csa =sym ans = .138223758410852e-16
二、符号变量 1、符号表达式允许使用自由变量。确定自 由变量的原则: 1)小写字母i和j不能作为自由变量。 2)符号表达式中如果有多个字符变量,则按照以下顺序选择自由变量: 首先选择x作为自由变量;如果没有x,则选择在字母顺序中最接近x的字符变量;如果与x相同距离,则在x后面的优先。
【例2.1-2】 用符号计算研究方程 的解。 (1)不指定变量(默认) syms u v w z Eq=u*z^2+v*z+w; result_1=solve(Eq) findsym(Eq,1) 2)指定变量为Z result_2=solve(Eq,z) result_1 = -u*z^2-v*z ans = w result_2 = 1/2/u*(-v+(v^2-4*u*w)^(1/2)) 1/2/u*(-v-(v^2-4*u*w)^(1/2))
2、确定自由变量的指令 findsym (11a版为symvar)的格式为: findsym(EXPR,n) 功能:确定EXPR中的自由变量。 其中EXPR可以是符号表达式或符号矩 阵;n为按顺序得出符号变量的个数, 当n省略时,则不按顺序给出EXPR中所 有的符号变量。 注意:(11版指令换成:symvar)
【例】findsym确定自由变量是对整个矩阵进行的。【例】findsym确定自由变量是对整个矩阵进行的。 syms a b t u v x y A=[a+b*x,sin(t)+u;x*exp(-t),log(y)+v] findsym(A,1) A = [ a+b*x, sin(t)+u] [ x*exp(-t), log(y)+v] ans = x
2.1.1符号计算中的算符 与数字计算的符号相同 2.1.2符号计算中的函数指令 见表2.1-1,P45 注意与数值计算的函数和指令的异同
2.1.3 符号对象的识别 例2.1-5 数据对象及其识别指令的使用。 (1) clear a=1;b=2;c=3;d=4; Mn=[a,b;c,d] ; Mc='[a,b;c,d]'; Ms=sym(Mc) ; (2) SizeMn=size(Mn) SizeMc=size(Mc) SizeMs=size(Ms) (3) CMn=class(Mn) CMc=class(Mc) CMs=class(Ms) SizeMn = 2 2 SizeMc =1 9 SizeMs =2 2 CMn =double CMc =char CMs =sym
2.2符号数字及表达式的操作 • 2.2.1符号数字的任意精度计算 • 1、digits(n) • 功能:设定有效位数函数。其中n为所 • 期望的有效位数,默认值为32位。 • 2、vpa(s,n) • 功能:将s表示为n位有效位数的符号 • 对象。
例:应用digits和vpa函数设置运算精度。 a=sym('2*sqrt(5)+pi') %创建符号对象 digits %显示默认的有效位数 vpa(a) %用默认的位数计算并显示 vpa(a,20) %按指定的精度计算并显示 digits(15) %改变默认的有效位数 vpa(a) %按digits指定的精度计算并显示 结果:a =2*sqrt(5)+pi Digits = 32 ans =7.6137286085893726312809907207421 ans =7.6137286085893726313 ans =7.61372860858937
2.2.2符号表达式的基本操作 指令:simple(EXPR) 【例2.2-2】简化。 syms x f=(1/x^3+6/x^2+12/x+8)^(1/3); g1=simple(f) g2=simple(g1) g1 = (2*x+1)/x g2 = 2+1/x
2.2.3表达式中的置换操作 1、通用置换指令 subs (S) 功能:用MATLAB工作空间中的变量替换S符号 表达式中的所有变量。 subs (S,NEW) 功能:用变量NEW替换符号表达式S中的自由 变量。 subs (S,OLD,NEW) 功能:用变量NEW替换符号表达式S中的变量 OLD。
【例2.2-4】演示subs的置换规则。 • (1)syms a x • f=a*sin(x)+5 • (2)f1=subs(f,'sin(x)',sym('y')) • class(f1) • (3)f2=subs(f,{a,x},{2,sym('pi/3')}) • class(f2) • (4)f3=subs(f,{a,x},{2,pi/3}) • class(f3) • (5)数组置换 • f4=subs(subs(f,a,2),x,0:pi/6:pi) • class(f4) • (6)f5=subs(f,{a,x},{0:6,0:pi/6:pi}) • class(f5) f1 =a*y+5 ans =sym f2=3^(1/2)+5 ans =sym f3 =6.7321 ans =double f4 = 5.0000 6.0000 6.7321 7.0000 6.7321 6.0000 5.0000 ans =double
2.3符号微积分 2.3.1极限和导数的符号计算 1、求函数极限的函数:limit limit(f,x,a) 功能:求符号函数f(x)的极限值。 即计算当自变量x趋近于常数a时,f(x)函数的极 限值。 limit(f,x,a,'right') 功能:求符号函数f的极限值。 'right'表示变量x从右边趋近于a。 limit(f,x,a,'left') 功能:求符号函数f的极限值。 'left'表示变量x从左边趋近于a。
【例2.3-1】试求 syms x k Lim_f=limit((1-1/x)^(k*x),x,inf) Lim_f =exp(-k)
2、求符号表达式的微分的函数:diff diff(f) 功能:求f对自由变量的一阶微分 diff(f,x) 功能:求f对符号变量x的一阶微分 diff(f,n) 功能:求f对自由变量的n阶微分 diff(f,x,n) 功能:求f对符号变量t的n阶微分
syms a t x; f=[a,t^3;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/x] dfdt2 = [ 0, 6*t] [ 0, 0] dfdxdt = [ 0, 0] [ -sin(x), 0]
clear syms x g=sym('cos(x+sin(y(x)))=sin(y(x))') dgdx=diff(g,x) %结果为:dgdx = -sin(x+sin(y(x)))*(1+cos(y(x))*diff(y(x),x)) = cos(y(x))*diff(y(x),x) dgdx1=subs(dgdx,'diff(y(x),x)','dydx') % diff(y(x),x))不能作为变量名 dydx=solve(dgdx1,'dydx') %关于dydx求解,使dy/dx通过x,y表示 dydx =-sin(x+sin(y(x)))/cos(y(x))/(sin(x+sin(y(x)))+1)
3、r=taylor(f,n,x,a) 功能:把f(x)在x=a处展开成n次幂级数,以上x缺少时由findsys自动确认,n缺少时默认为1 syms x r=taylor(x*exp(x),9,x,0) %忽略9阶及9阶以上小量的展开 求f(x)=sin(x)在x=0处的泰勒级数展开?
2.3.1序列/级数的符号求和 级数求和运算函数:symsum symsum(f,v,a,b) 功能:计算符号表达式f的级数和。 其中f为通项式,v自变量,v省略则默认为对自由变量求和;[a,b]为参数v的取值范围,缺少时求和区间默认为[0~v-1]
syms k t; f1=[t k^3]; f2=[1/(2*k-1)^2,(-1)^k/k]; s1=simple(symsum(f1)) s2=simple(symsum(f2,1,inf)) s1 = [ 1/2*t*(t-1), k^3*t] s2 = [ 1/8*pi^2, -log(2)]
2.3.2符号积分 符号积分函数:int int(f,v,a,b) 功能:以v为自变量,以a、b分别表示定积分的下 限和上限,对被积函数的符号表达式f求定积分。 int(f,v) 功能:以v为自变量,对被积函数的符号表达式f求 不定积分。 注意:没有指定积分变量v时,按findsym函数确 定的默认变量对被积函数的符号表达式f求积分. int函数的嵌套使用可实现二重积分的计算。
clear syms x f=sqrt((1+x)/x)/x s=int(f,x) s=simple(simple(s)) f =((1+x)/x)^(1/2)/x s = log(1/2+x+((1+x)*x)^(1/2))-2*((1+x)*x)^(1/2)/x
syms a b x; f=[a*x,b*x^2;1/x,sin(x)]; disp('The integral of f is'); pretty(int(f)) %一行显示变两行显示 直接执行int(f)时显示: [ 1/2*a*x^2, 1/3*b*x^3] [ log(x), -cos(x)] The integral of f is [ 2 3] [1/2 a x 1/3 b x ] [ ] [ log(x) -cos(x) ]
syms x y z F2=int(int(int(x^2+y^2+z^2,z,sqrt(x*y),x^2*y),y,sqrt(x),x^2),x,1,2) VF2=vpa(F2) F2 = 14912/4641*2^(1/4)-6072064/348075*2^(1/2)+64/225*2^(3/4)+1610027357/6563700 VF2 = 224.92153573331143159790710032805
2.4 微分方程的符号解法 2.4.2 求微分方程符号解的一般指令 • 符号微分方程的求解函数:dsolve dsolve('eqn1','eqn2',...,'eqn',’condition1’,‘ condition2’, …,‘conditionN’,'v’) 功能:求微分方程组eqn1,eqn2, ...,eqn的通解。初值条件为condition1, condition2,…, conditionN下,变量为v (缺少时默认为t)。
dsolve函数也可求解二阶微分方程。 clear S=dsolve('Dx=y,Dy=-x') disp(['解',blanks(12),'x',blanks(21),'y']) disp([S.x,S.y]) S = x: [1x1 sym] y: [1x1 sym] 解 x y [ -C1*cos(t)+C2*sin(t), C1*sin(t)+C2*cos(t)]
(1)指令 y=dsolve('x*D2y-3*Dy=x^2','y(1)=0,y(5)=0','x') y =31/468*x^4-1/3*x^3+125/468 (2)图型化 ezplot(y,[-1,6]) %画y函数的图像,范围[-1~6] hold on plot([1,5],[0,0],'.r','MarkerSize',20) %在(1,0)和(5,0)位置画出红点 text(1,1,'y(1)=0') text(4,1,'y(5)=0') title(['x*D2y-3*Dy=x^2',', y(1)=0,y(5)=0']) hold off
2.5符号变换和符号卷积 2.5.1Fourier变换及其反变换 正变换fourier F=fourier(f,t,w) 功能:返回函数f(t)的fourier变换F。其中返回结 果F是符号变量w的函数,当参数w省略,默认返 回结果为w的函数;当参数t省略,默认自由变量为x。 傅里叶反变换函数:ifourier f=ifourier(F,w,t) 功能:返回函数F(w)的fourier反变换f(t)。参数 含义同fourier函数。
(1)傅里叶变换 syms t w ut=heaviside(t); %阶跃函数 UT=fourier(ut,t,w) (2)反变换验算 Ut=ifourier(UT,w,t) UT = pi*dirac(w)-i/w Ut = heaviside(t)
(1)求傅里叶变换 syms A t w syms tao positive yt=heaviside(t+tao/2)-heaviside(t-tao/2); Yw=fourier(A*yt,t,w) (2)反变换验算 Yt=ifourier(Yw,w,t)
(3)图形展示 设tao=3,A=1 yt3=subs(yt,tao,3) Yw3=subs(Yw,[A,tao],[1,3]) subplot(2,1,1) Ht=ezplot(yt3,[-3,3]);%图形句柄Ht set(Ht,'Color','r','LineWidth',3) subplot(2,1,2),ezplot(Yw3)
2.5.2 Laplace变换及其反变换 Laplace和ilaplace F=laplace(f,t,s) 功能: 求函数f(t)的Laplace变换F(s)。 其中返回结果F为s的函数 f=ilaplace(F,s,t) 功能: 求函数F(s)的iLaplace变换f(t)。参数的含义同上laplace函数。
syms t s ; syms a b positive Dt=dirac(t-a); Ut=heaviside(t-b); Mt=[Dt,Ut;exp(-a*t)*sin(b*t),t^2*exp(-t)]; MS=laplace(Mt,t,s) MS = [ exp(-s*a), exp(-s*b)/s] [ 1/b/((s+a)^2/b^2+1), 2/(s+1)^3]
2.5.3 Z变换及其反变换 ztrans和iztrans FZ=ztrans(fn,n, z) 功能:求时域序列fn的Z变换FZ。 FZ以符号变量z为自变量 fn=iztrans(FZ,z,n) 功能:求Z域函数FZ的逆变换函数fn fn为时域离散序列,变量为n
【例2.5-5】一组Z变换、反变换算例 (1)求6(1-0.5^n)序列的Z变换 clear syms n z clear gn=6*(1-(1/2)^n) G=simple(ztrans(gn,n,z)); pretty(G) (2) clear syms n z clear delta=sym(‘kroneckerDelta(n, 0)’); %定义单位脉冲序列 KD=ztrans(delta,n,z) inv_KD=iztrans(KD) gn = 6 - 6*(1/2)^n 6 z -------------- 2 2 z - 3 z + 1 KD = 1 inv_KD = kroneckerDelta(n, 0) 冲激序列的Z变换及反变换
2.5.4 符号卷积 syms T t tao ut=exp(-t); ht=exp(-t/T)/T; uh_tao=subs(ut,t,tao)*subs(ht,t,t-tao); yt=simple(simple(int(uh_tao,tao,0,t))) yt = -(1/exp(t) - 1/exp(t/T))/(T - 1)
【例2.5-8】采用Laplace变换和反变换求上例的输出响应。【例2.5-8】采用Laplace变换和反变换求上例的输出响应。 syms s yt=ilaplace(laplace(ut,t,s)*laplace(ht,t,s),s,t); yt=simple(yt) syms tao t=sym('t','positive'); ut=heaviside(t)-heaviside(t-1); ht=t*exp(-t); yt=int(subs(ut,t,tao)*subs(ht,t,t-tao),tao,0,t)
2.6 符号矩阵分析和代数方程解 • 2.6.1符号矩阵分析 • 常用指令: • Det(A): 求矩阵A的行列式 • Diag(A): 取对角元素构成向量 • [V,D]=eig(A): 特征值分解,使AV=VD • expm(A): 矩阵指数eA • inv(A): A的逆 • ploy(A): 矩阵的特征多项式 • rank(A): 矩阵的秩
syms a11 a12 a21 a22 A=[a11,a12;a21,a22] DA=det(A) IA=inv(A) EA=eig(A) A = [ a11, a12] [ a21, a22] DA = a11*a22-a12*a21 IA = [ a22/(a11*a22-a12*a21), -a12/(a11*a22-a12*a21)] [ -a21/(a11*a22-a12*a21), a11/(a11*a22-a12*a21)] EA = 1/2*a11+1/2*a22+1/2*(a11^2-2*a11*a22+a22^2+4*a12*a21)^(1/2) 1/2*a11+1/2*a22-1/2*(a11^2-2*a11*a22+a22^2+4*a12*a21)^(1/2)
2.6 线性方程组的符号解 • 一般代数方程组的解 • Solve指令 • S=solve(‘eq1’,’eq2’ …’eqn’,’v1’,’v2’ …’vn’) syms d n p q eq1=d+n/2+p/2-q; eq2=n+d+q-p-10; eq3=q+d-n/4-p; S=solve(eq1,eq2,eq3,d,n,p,q); disp([' S.d',' S.n',' S.p',' S.q']) disp([S.d,S.n,S.p,S.q]) S.d S.n S.p S.q [ d, 8, 4*d+4, 3*d+6]
S=solve('u*y^2+v*z+w=0','y+z+w=0','y','z') disp('S.y'),disp(S.y) disp('S.z'),disp(S.z) S = y: [2x1 sym] z: [2x1 sym] S.y -1/2/u*(-2*u*w-v+(4*u*w*v+v^2-4*u*w)^(1/2))-w -1/2/u*(-2*u*w-v-(4*u*w*v+v^2-4*u*w)^(1/2))-w S.z 1/2/u*(-2*u*w-v+(4*u*w*v+v^2-4*u*w)^(1/2)) 1/2/u*(-2*u*w-v-(4*u*w*v+v^2-4*u*w)^(1/2))
clear all,syms x; s=solve('(x+2)^x=2','x') s = .69829942170241042826920133106081 本方程无解析解,只有尽似解
2.8 符号计算结果的可视化 注意:此部分仅仅作为数值化绘图指令的补充, 故不作详细介绍 • 2.8.1直接可视化符号表达式 • 见表2.8-1(P95) • 1.单独立变量符号函数的可视化 • Ezplot(fx,[xmin,xmax,ymin,ymax]): • 二维曲线 • Ezplot3(xt,yt,zt,[tmin,tmax]): • 三维曲线
syms t tao y=2/3*exp(-t/2)*cos(sqrt(3)/2*t); s=subs(int(y,t,0,tao),tao,t); %s=int(y,t,0,t) 错误 subplot(2,1,1) ezplot(y,[0,4*pi]),ylim([-0.2,0.7]) %确定y轴范围 grid on subplot(2,1,2) ezplot(s,[0,4*pi]) grid on title('s = \inty(t)dt')
2 双独立变量符号函数的可视化 ezsurf(x,y,z,[坐标范围]):曲面 【例2.8-2】使用球坐标参量画部分球壳。 x='cos(s)*cos(t)'; y='cos(s)*sin(t)'; z='sin(s)'; ezsurf(x,y,z,[0,pi/2,0,3*pi/2]) view(17,40) shading interp colormap(spring) light('position',[0,0,-9],'style','local') light('position',[-1,-0.5,2],'style','local') material([0.5,0.5,0.5,10,0.3])