数值计算和符号计算.ppt

《数值计算和符号计算.ppt》由会员分享，可在线阅读，更多相关《数值计算和符号计算.ppt（25页珍藏版）》请在三一办公上搜索。

1、1,五、数值计算和符号计算,两种计算的特点数值计算符号对象和符号表达式符号计算符号函数的可视化Maple函数的使用,2,5.1 两种计算的特点,数值计算特点：1）以数值数组作为运算对象，给出数值解；2）计算过程中产生误差累积问题，影响计算结果的精确性；3）计算速度快，占用资源少。,符号计算特点：1）以符号对象和符号表达式作为运算对象，给出解析解；2）运算不受计算误差累积问题的影响；3）计算指令简单；4）占用资源多，计算耗时长。,5.2 数值计算,MATLAB具有强大的数值计算功能，可完成矩阵分析、线性代数、多元函数分析、数值微积分、方程求解、边值问题求解、数理统计等常见的数值计算。,数值计算的

2、常用运算单元是数值数组。MATLAB给出了大量的数值计算函数，基本上与理论数学、数值数学的数学描述式表达方式相同，便于编程和掌握。,本节主要以例题的形式给出一些常用的数值计算问题的MATLAB解算过程，以便熟悉MATLAB的计算指令。相对于具体的应用环境，需要根据实际情况查阅MATLAB函数列表，选择合适的函数和参数进行处理。,3,【例51】矩阵常见运算,%exm05_01.mA=1 2 3;4 7 2;7 4 3;%A为33矩阵b=2;4;5;%b为13矩阵%矩阵的分解L U P=lu(A)%矩阵的LU分解，分成下三角和上三角阵,LU=PAQ,R=qr(A)%矩阵的QR分解,分成正交方阵和上

3、三角阵,A=QR,%矩阵的特征参数Adet=det(A)%求矩阵的行列式Arank=rank(A)%求矩阵的秩Anorm=norm(A)%求矩阵的范数，通过带不同的参数可以求不同的范数P=poly(A)%求矩阵特征多项式Aroots=roots(P)%求特征根Aroots2=eig(A)%特征根的又一种求法,%线性方程组求解x=Ab%求方程组AXb的解,x=inv(A)*b,4,【例52】求函数的零点,%exm05_02.my=inline(cos(t)*exp(-0.1*t)-0.1,t)%构造内联函数 ye-0.1t cost-0.1t=-10:0.01:10;%对自变量采样，采样步长不宜

4、太大。yt=feval(vectorize(y),t);%计算得到相应的 y 值ysign=sign(yt);%利用sign函数判断正负变化n=0;%找出ysign中发生正负变化的下标，即yt中穿越y=0水平线前一时刻的下标for i=2:length(t)if ysign(i)=ysign(i-1)n=n+1;yzero(n)=i-1;%与前一函数值符号相反，则表示有一零点 end%yzero(n)存放第n个零点对应的下标end,n,yzeron=6yzero=220 523 852 1146 1488 1764,%对应yzero中每一个下标找出穿越零水平线前后的两个时间%利用fzero寻找

5、该时间范围内的精确零点的横坐标for n=1:length(yzero)index=yzero(n);z(n)=fzero(y,t(index),t(index+1);%在区间t(index),t(index+1)内寻找零点end,内联函数是MATLAB可实现函数功能的一个对象.其输入变量不能是数组,但可用命令vectorize使其适用于数组运算.可直接用feval命令执行内联函数.,5,z%显示零点横坐标plot(t,yt)hold onplot(t,zeros(1,length(t),k)%画一黑横线(y=0)xx=ginput(6)%用鼠标从图形上获取n个点的坐标(t,y),z=-7.8

6、082-4.7745-1.4845 1.4549 4.8760 7.6377,用鼠标取得的坐标：xx=-7.8111 0.0044-4.7696 0.0044-1.4977 0.0044 1.4977 0.0044 4.9078 0.0044 7.6728 0.0044,6,clear%被积函数y=inline(exp(-0.8*t.*abs(sin(t),t);n=1;%将 x 取值离散化，分别求出 x 等于不同的值时的积分结果x=0:0.1:10;for xi=x yint(n)=quad(y,0,xi);n=n+1;end,【例53】数值积分:,%以折线拟和离散的积分结果，得到 y(x)

7、的近似曲线plot(x,yint),quad(exp(-0.8*x.*abs(sin(x),0,10)ans=2.6597,7,%exm05_04.mclear%得到掺杂了均值为0，方差为0.2的高斯白噪声的信号y%原始 y 信号有三个频率分量：10，100，180 Hzrandn(state,0);Ts=0.002;t=0:Ts:0.4;y=sin(2*pi*10*t)+cos(2*pi*100*t)+2*sin(2*pi*180*t)+0.2*randn(size(t);%利用FFT求 y 的离散 Fourier 系数 YY=fft(y);ws=2*pi/Ts;%采样角频率wn=ws/2;

8、%Nyquest频率%得到10阶ButterWorth带通滤波器，截止频率是90110HzB,A=butter(10,90/wn,110/wn);y100=filter(B,A,y);%滤波，得到100Hz的信号subplot(2,1,1)plot(t,y,b-,t,y100,r-)%作图显示噪声信号和滤波后信号axis(0.2,0.3,-4,4)subplot(2,1,2)w=linspace(0,wn,length(t)/2);%半采样频率中相应的刻度Ya=Ts*abs(Y(1:length(t)/2);plot(w,Ya)%做出半边频谱,【例54】Fourier分析,8,5.3 符号对象

9、和符号表达式,【例55】符号对象入门,syms a11 a12 a21 a22;%声明基本符号对象A=a11,a12;a21,a22A=a11,a12 a21,a22,参加数值计算的数值数组必须先进行赋值才能参加计算，参加符号运算的基本符号也必须先经过定义，才能构成符号表达式，进而进行相应的符号计算。,class(A)ans=sym,isa(A,sym)ans=1,DA=det(A)%行列式的值DA=a11*a22-a12*a21,IA=inv(A)%求A的逆IA=a22/(a11*a22-a12*a21),-a12/(a11*a22-a12*a21)-a21/(a11*a22-a12*a21

10、),a11/(a11*a22-a12*a21),9,5.3.1 符号对象的生成和使用,定义基本符号对象的命令是：sym和syms。,1.创建符号常量,符号常量是不含变量的符号表达式，用sym命令来创建符号常量。,sym(常量)创建符号常量 sym(常量,参数)把常量按“参数”指定的格式转换为符号常量,2.创建符号变量和表达式,sym(变量,参数)把变量定义为符号对象，参数用来设置限定符号变量 的数学特性，可以选择为positive、real和unreal。sym(表达式)创建符号表达式,syms(arg1,arg2,参数)把字符变量定义为符号变量syms arg1 arg2,参数 把字符变量定

11、义为符号变量的简洁形式 syms用来创建多个符号变量，这两种方式创建的符号对象是相同的。参数设置和前面的sym命令相同，省略时符号表达式直接由各符号变量组成。,f1=sym(a*x2+b*x+c)f1=a*x2+b*x+c,syms a b c x f2=a*x2+b*x+c f2=a*x2+b*x+c,syms(a,b,c,x)f3=a*x2+b*x+cf3=a*x2+b*x+c,符号计算步骤：定义符号对象；构造符号表达式；进行符号计算。,10,5.3.2 符号表达式及符号函数的操作和转换,1.符号表达式中自由变量的确定,自由变量的确定原则：小写字母i和j不能作为自由变量。符号表达式中如果有

12、多个字符变量，则按照以下顺序选择自由变量：首先选择x作为自由变量；如果没有x，则选择在字母顺序中最接近x的字符变量；如果与x相同距离，则在x后面的优先。大写字母比所有的小写字母都靠后。,如果不确定符号表达式中的自由符号变量，可以用findsym函数来自动确定。findsym(EXPR,n)确定自由符号变量 EXPR可以是符号表达式或符号矩阵；n为按顺序得出符号变量的个数，当n省略时，则不按顺序得出EXPR中所有的符号变量。,K=3;f=sym(a*x2+b*x+c+2*K);findsym(f)%得出所有的符号变量 ans=K,a,b,c,x,g=sym(sin(z)+cos(v)g=sin(

13、z)+cos(v),findsym(g,1)%得出第一个符号变量 ans=z,符号变量z和v距离x相同，以在x后面的z为自由符号变量。,syms a b c x;k=3;f=a*x2+b*x+c+2*k;findsym(f)ans=a,b,c,x,11,2.符号表达式的因式分解、简化和展开,同一个数学函数的符号表达式可以表示成三种形式，如 多项式形式的表达方式：f(x)=x3-6x2+11x-6 因式形式的表达方式：f(x)=(x-1)(x-2)(x-3)嵌套形式的表达方式：f(x)=x(x(x-6)+11)-6,collect(EXPR,v)合并v的幂类项，如省略v则默认xfactor(EX

14、PR)因式分解simple（EXPR）化简表达式EXPR,f1=sym(x3+2*x2*y+4*x*y+6)f1=x3+2*x2*y+4*x*y+6collect(f1,y)%按y来合并同类项 ans=(2*x2+4*x)*y+x3+6,f=sym(x3-6*x2+11*x-6);factor(f)ans=(x-1)*(x-2)*(x-3),f2=sym(x2-a2)/(x-a);simple(f2)ans=x+a,simple 函数给出多种简化形式，给出除了pretty、collect、expand、factor、simplify简化形式之外的radsimp、combine、combine(

15、trig)、convert形式，并寻求包含最少数目字符的表达式简化形式。,12,3.求反函数和复合函数,finverse(f,v)对指定自变量v的函数f(v)求反函数 当v省略，则对默认的自由符号变量求反函数。,fg=compose(f,g)把g作为f的自变量得到复合函数fg,f=sym(t*ex);%原函数 tex g=finverse(f)%对默认自由变量求反函数 g=log(x/t)/log(e),g=finverse(f,t)%对t求反函数,若先定义t为符号变量，则单引号可去掉 g=t/(ex),f=sym(t*ex);%创建符号表达式g=sym(a*y2+b*y+c);%创建符号表达

16、式h1=compose(f,g)%计算f(g(x)h1=t*e(a*y2+b*y+c),h2=compose(g,f)%计算g(f(x)h2=a*t2*(ex)2+b*t*ex+c,h3=compose(f,g,z)%计算f(g(z)h3=t*e(a*z2+b*z+c),13,4.符号表达式的替换,RS,ssub=subexpr(S,ssub)用符号变量ssub来置换S中的子表达式,得到RS subexpr函数对子表达式是自动寻找的，只有比较长的子表达式才被置换，比较短的子表达式，即使重复出现多次，也不被置换。,syms a b c d D IA=inv(a b;c d)IA=d/(a*d-b

17、*c),-b/(a*d-b*c)-c/(a*d-b*c),a/(a*d-b*c),syms a11 a12 a13 a21 a22 a23 a31 a32 a33 D;A=a11,a12,a13;a21,a22,a23;a31,a32,a33;IA=inv(A);IAS,D=subexpr(IA,D)IAS=(a22*a33-a23*a32)/D,-(a12*a33-a13*a32)/D,(a12*a23-a13*a22)/D-(a21*a33-a23*a31)/D,(a11*a33-a13*a31)/D,-(a11*a23-a13*a21)/D(a21*a32-a22*a31)/D,-(a1

18、1*a32-a12*a31)/D,(a11*a22-a12*a21)/DD=a11*a22*a33-a11*a23*a32-a21*a12*a33+a21*a13*a32+a31*a12*a23-a31*a13*a22,IAS,D=subexpr(IA,D)IAS=d/(a*d-b*c),-b/(a*d-b*c)-c/(a*d-b*c),a/(a*d-b*c)D=empty sym,14,subs函数可用来进行对符号表达式中符号变量的替换。RES=subs(ES)用给定值替换符号表达式ES中的所有变量后产生RESRES=subs(ES,new)用new替换符号表达式ES中的自由变量后产生RES

19、RES=subs(ES,old,new)用new替换符号表达式ES中的old变量后产生RES,【例56】符号置换,%exm05_06.mclearsyms a x y;f=a*sin(x)+5;f1=subs(f,sin(x),y);%用y替换f中的sin(x)f2=subs(f,a,x,2,sym(2*pi/3);%用2和符号变量2*pi/3分别替换f中的a和xf1,f2f1=a*y+5f2=3(1/2)+5,class(f2),double(f2)ans=symans=6.7321,f3=subs(f,a,x,0:3,0:pi/3:pi)%用数值数组替f3=%换符号变量 5.0000 5.

20、8660 6.7321 5.0000class(f3)ans=double,15,5.3.3 符号对象和其他数据对象之间的转换,符号、数值、字符串间的转换：,符号多项式、系数向量、字符串多项式之间的转换：,syms x;F=x3+2*x+5;,F=x3+2*x+5,F1=1 0 2 5,F3=x3+2*x+5,3 x+2 x+5,F=sym(x3+2*x+5)F2=char(F),F1=sym2poly(F)；,F2=poly2str(F1,x);,F3=poly2sym(F1);,pretty(F),16,5.4 符号计算,5.4.1 符号微分,diff(f)求f对自由变量的一阶微分 dif

21、f(f,v)求f对符号变量v的一阶微分 diff(f,n)求f对自由变量的n阶微分 diff(f,v,n)求f对符号变量v的n阶微分,5.4.2 符号积分,int(f,v)求符号变量v的不定积分,当v省略时则为默认自由变量 int(f,v,a,b)求符号变量v的积分,a和b为数值,a,b为积分区间 int(f,v,m,n)求符号变量v的积分,m和n为符号对象，m,n为积分区间,【例57】符号微积分,syms a b t x;f=a t3;t*cos(x)log(x);%22符号矩阵df=diff(f)%求f对x的微分df=0,0-t*sin(x),1/x,dfdxdt=diff(diff(f,

22、x),t)dfdxdt=0,0-sin(x),0,dfdt2=diff(f,t,2)%求f对t的二阶微分dfdt2=0,6*t 0,0,F1=int(f)%对f求不定积分F1=a*x,t3*x t*sin(x),x*log(x)-x,F2=int(f,x,0,b)%在0,b对f求定积分F2=a*b,t3*b t*sin(b),b*log(b)-b,17,【例58】符号卷积:,syms T t tao;%定义符号变量ut=exp(-t);%生成u(t)ht=exp(-t/T)/T;%生成h(t)uh=subs(ut,t,tao)*subs(ht,t,t-tao)%产生被积函数uh=u(tao)*

23、h(t-tao)uh=exp(-tao)*exp(-(t-tao)/T)/T,卷积定义式：,5.4.3 符号积分变换,1.fourier变换及其反变换,Fwfourier(ft,t,w)求时域函数ft的fourier变换Fw ft=ifourier(Fw,w,t)求频域函数F的fourier反变换f(t)时域中的默认变量是x，频域中的默认变量是w,yt=int(uh,tao,0,t)%在区间0,t对uh求定积分yt=-(exp(-t)-exp(-t/T)/(T-1),yt=subs(yt,T,0.5)%令y(t)中T=0.5yt=2*exp(-t)-2*exp(-2*t),ezplot(yt,

24、0,2*pi)%画出函数y(t)在区间0,2*pi上的波形,18,【例59】Fouier变换,%宽度为T，高度为1的矩形脉冲的Fourier变换syms t T w;ft1=sym(A*(Heaviside(t+T/2)-Heaviside(t-T/2);%用阶跃函数表示矩形脉冲Fw1=fourier(ft1,t,w)%求ft1的Fourier变换Fw1=A*(exp(1/2*i*T*w)*(pi*Dirac(w)-i/w)-exp(-1/2*i*T*w)*(pi*Dirac(w)-i/w),%由Fourier变换定义式,直接用积分命令求syms T A w;Fw=int(A*exp(-i*w

25、*t),t,-T/2,T/2);Fw=simple(Fw)Fw=2*A*sin(1/2*T*w)/w,Fw1=maple(convert,Fw1,piecewise,w)%进入maple工作空间计算，结果返回Fw1=-i*A*exp(1/2*i*T*w)/w+i*A*exp(-1/2*i*T*w)/w,Fw1=simple(Fw1)%化简表达式Fw1Fw1=2*A*sin(1/2*T*w)/w,Fw1=subs(Fw1,T,A,1,1)Fw1=2*sin(1/2*w)/w,Heaviside(t)和Dirac(t)在MATLAB符号资源中代表了单位阶跃函数和单位冲激函数。但由于MATLAB仅借

26、用了这两个函数的数学定义，它们可以在符号计算中完成正确的计算过程，却无法与其它指令相结合使用，如绘图指令。,19,2.Laplace变换及其反变换,Fs=laplace(ft,t,s)求时域函数ft的Laplace变换Fs ftilaplace(Fs,s,t)求Fs的Laplace反变换ft,3.Z变换及其反变换,Fzztrans(fn,n,z)求时域序列fn的Z变换Fzfniztrans(Fz,z,n)求Fz的z反变换fn,【例510】Laplace变换、Z变换,syms a t sFs1=laplace(sin(a*t),t,s)%求sin(at)的Laplace变换 Fs1=a/(s2+

27、a2),f1=ilaplace(1/(s+a),s,t)%求1/s+a的Laplace反变换 f1=exp(-a*t),syms z n;fn=sym(2*charfcn0(n)+6*(1-(1/2)n);Fz=simple(ztrans(fn,n,z)%求f(n)的Z变换Fz=(4*z2+2)/(2*z2-3*z+1)Fz_n=iztrans(Fz,z,n)%求Fz的反Z变换Fz_n=2*charfcn0(n)+6-6*(1/2)n,charfcnA(x)是Maple定义的指征函数，定义式为利用charfcn可以表示冲激序列：charfcn0(n),charfcn3(n)等。,20,5.4.

28、4 符号方程的求解,solve(eq,v)求方程关于指定变量的解solve(eq1,eq2,v1,v2,)求方程组关于指定变量的解 eq可以是含等号的符号表达式的方程，也可以是不含等号的符号表达式，但所指的仍是令eq=0的方程；当参数v省略时，默认为方程中的自由变量；其输出结果为结构数组类型。,dsolve(eq,con,v)求解微分方程dsolve(eq1,eq2,con1,con2,v1,v2)求解微分方程组 eq为微分方程；con是微分初始条件，可省略；v为指定自由变量，省略时则默认为x或t为自由变量；输出结果为结构数组类型。,【例511】符号代数方程求解,%数值解法A=1,0.5,-1

29、;1 0 1;-2-1 1;b=0;2;1;%代数方程Ax=bx=Ab,%符号解法的两种调用格式S=solve(x+1/2*y-z=0,x+z=2,-2*x-y+z=1)%返回结构数组S.x=3 S.y=-8S.z=-1,x=3-8-1,21,x,y,z=solve(x+1/2*y-z=0,x+z=2,-2*x-y+z=1)%返回符号数值x,y,zx=3y=-8z=-1,【例512】符号微分方程求解,%求解微分方程，含有待定系数S=dsolve(Dx=y,Dy=-x);%默认独立变量t,dx/dt=y,dy/dt=-xdisp(S.x=,char(S.x),disp(S.y=,char(S.y

30、)S.x=cos(t)*C1+sin(t)*C2S.y=-sin(t)*C1+cos(t)*C2,%给定初始条件S=dsolve(Dx=y,Dy=-x,x(0)=-0.5,y(0)=0.5,s);%dx/ds=y,dy/ds=-x,x(0)=y(0)=0.5disp(S.x=,char(S.x),disp(S.y=,char(S.y)S.x=-1/2*cos(s)+1/2*sin(s)S.y=1/2*sin(s)+1/2*cos(s),%二阶微分方程求解S=dsolve(x*D2y-3*Dy=x2,y(1)=0,y(5)=0,x)%xd2y/dx2-3dy/dx=x2S=-1/3*x3+125

31、/468+31/468*x4,1)输入宗量必须是字符形式，其内容可以包含三个内容：微分方程、初始条件、独立变量。其中只有微分方程是必须的；Dny表示y的n阶微分，Dy表示y的一阶微分；2)初始条件以字符串的形式给出，如:y(0)=1。初始条件不足时，解中包含有任意常数C1,C2,;3)最后一个输入宗量指出独立变量，默认是t。,22,5.5 符号函数的可视化,5.5.1 符号函数的绘图命令,1.ezplot和 ezplot3命令,ezplot命令是绘制符号表达式的自变量和对应各函数值的二维曲线，ezplot3命令用于绘制三维曲线。,ezplot(F,xmin,xmax,fig)画符号表达式的图形

32、 F是将要画的符号函数；xmin,xmax是绘图的自变量范围，省略时默认值为2，2；fig是指定的图形窗口，省略时默认为当前图形窗口。,x=sym(sin(t);z=sym(t);y=sym(cos(t);ezplot3(x,y,z,0,10*pi,animate)绘制t在0,10*pi范围的三维曲线,2.其它常用绘图命令,23,5.5.2 图形化的符号函数计算器,Symbolic Math Toolbox还提供了另一种符号计算方式即图形化的符号函数计算器，由funtool.m文件生成。在MATLAB命令窗口输入命令“funtool”，就会出现该图形化函数计算器。,Figure No.1:f(

33、x)波形Figure No.2:g(x)波形Figure No.3:由鼠标控制界面操作，可对f(x)和g(x)进行处理并将结果在图形窗口1和窗口2中图示出来。,24,5.6 Maple函数的使用,MATLAB的符号运算实际上是由Maple软件所支持的符号数学工具箱完成的，Maple软件有2000多条符号运算命令。MATLAB提供了五个与MAPLE之间进行交互的命令。,5.6.1 访问Maple函数,1.maple函数,maple函数用于进行符号运算，并将计算结果返回到MATLAB的工作空间。maple(MapleStatement)运行Maple格式的语句MapleStatement mapl

34、e(fun,arg1,arg2)运行以arg1,arg2为参数的Maple的fun函数,maple(discrim(a*x2+b*x+c,x)%利用Maple的discrim函数，计算多项式的判别式 ans=-4*a*c+b2,2.mfun函数,mfun函数用于对Maple中的经典函数进行数值运算。mfun(fun,p1,p2,)fun为函数名；p1,p2,为函数的参数,mfun(gcd,20,30)%利用gcd函数计算最大公约数 ans=10,syms a b c x%discrim函数的又一用法 maple(discrim,a*x2+b*x+c,x)ans=-4*a*c+b2,25,5.6

35、.2 获得Maple的帮助,1.mfunlist命令,mfunlist命令用来列出能被“mfun”命令计算的经典特殊Maple函数。,2.mhelp命令,mhelp命令用来寻求关于Maple库函数及其调用方法的帮助：(1)使用“mhelp index”可以查看Maple的索引目录。(2)使用“mhelp index 分类名”可以进一步深入查看Maple的某个具体类别。当输入“mhelp index”命令会出现如下表所示的类别名。,mhelp gcd 查看求最大公约数的函数gcd,gcd-greatest common divisor of polynomialslcm-least common multiple of polynomialsCalling Sequence:gcd(a,b,cofa,cofb)lcm(a,b,.)Parameters:a,b-multivariate polynomials over an algebraic number field or analgebraic function field.cofa,cofb-(optional)unevaluated names Description:,