matlab ch2例题doc文档格式.docx

1、x,y,theta 【例2.1-4】findsym确定自由变量是对整个矩阵进行的。syms a b t u v x yA=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【例2.1-5】数据对象及其识别指令的使用。cleara=1;b=2;c=3;d=4; Mn=a,b;c,d Mc=a,b;c,dMs=sym（Mc） Mn = 1 2 3 4Mc =c,dMs = a, b c, d SizeMn=size（Mn）SizeMc=size（Mc）SizeMs=size（M

2、s） SizeMn = 2 2SizeMc = 1 9SizeMs = 2 2 CMn=class（Mn）CMc=class（Mc）CMs=class（Ms） CMn =CMc =charCMs =sym isa（Mn,double）isa（Mc,charisa（Ms,sym） 1 1 （5）whos Mn Mc Ms Name Size Bytes Class Attributes Mc 1x9 18 char Mn 2x2 32 double Ms 2x2 312 sym 【例2.2-1】digits, vpa指令的使用。digits p0=sym（1+sqrt（5）/2） pr=sym（

3、1+sqrt（5）/2） %pd=sym（1+sqrt（5）/2,de32r=vpa（abs（p0-pr） ） e16=vpa（abs（p0-pd）,16） e32d=vpa（abs（p0-pd） Digits = 32p0 =（1+sqrt（5）/2pr =7286977268806824*2（-52）pd =1.6180339887498949025257388711907e32r =.543211520368251e-16e16 =0.e32d =.543211520368251e-16 【例2.2-2】简化。syms xf=（1/x3+6/x2+12/x+8）（1/3）;g1=simp

4、le（f）g2=simple（g1） g1 =（2*x+1）/xg2 =2+1/x 【例2.2-3】对符号矩阵进行特征向量分解。clear all syms a b c d WV,D=eig（a b;c d） RVD,W=subexpr（V;D,W） V = -（1/2*d-1/2*a-1/2*（d2-2*a*d+a2+4*b*c）（1/2）/c, -（1/2*d-1/2*a+1/2*（d2-2*a*d+a2+4*b*c）（1/2）/c 1, 1D = 1/2*d+1/2*a+1/2*（d2-2*a*d+a2+4*b*c）（1/2）, 0 0, 1/2*d+1/2*a-1/2*（d2-2*a*

5、d+a2+4*b*c）（1/2）RVD = -（1/2*d-1/2*a-1/2*W）/c, -（1/2*d-1/2*a+1/2*W）/c 1/2*d+1/2*a+1/2*W, 0 0, 1/2*d+1/2*a-1/2*WW =（d2-2*a*d+a2+4*b*c）（1/2） 数学背景：设A为n阶方阵,若存在数和非零的 n维列向量x,使得Ax=x或（A-I）x=0 ，则称数为矩阵 A的特征值,称 x为矩阵A对应于特征值的特征向量.【例2.2-4】用简单算例演示subs的置换规则。（1）产生符号函数syms a x;f=a*sin（x）+5 f =a*sin（x）+5 （2）符号表达式置换f1=s

6、ubs（f,sin（x）,sym（y） %class（f1） f1 =a*y+5（3）被双精度数字置换，被符号数字置换f2=subs（f,a,x,2,sym（pi/3） %class（f2） f2 =3（1/2）+5sym （4），均被双精度数值取代f3=subs（f,a,x,2,pi/3） %class（f3） f3 = 6.7321double （5）被数值数组取代f4=subs（subs（f,a,2）,x,0:pi/6:pi） %class（f4） f4 = 5.0000 6.0000 6.7321 7.0000 6.7321 6.0000 5.0000（6），被同样大小的数值数组取代，

7、这时执行的是数组乘f5=subs（f,a,x,0:6,0:pi） %class（f5） f5 = 5.0000 5.5000 6.7321 8.0000 8.4641 7.5000 5.0000表达式能和符号变量同时包含在old中被置换，如subs（f,a,sin（x）,2,sym（） 2*y+5 【例2.3-1】试求syms x kLim_f=limit（1-1/x）（k*x）,x,inf） Lim_f =exp（-k） 【例2.3-2】求求，syms a t x;f=a,t3;t*cos（x）, log（x）;df=diff（f） %求矩阵f对x的导数dfdt2=diff（f,t,2） %

8、求矩阵f对t的二阶导数dfdxdt=diff（diff（f,x）,t） %求二阶混合导数df = 0, 0 -t*sin（x）, 1/xdfdt2 = 0, 6*tdfdxdt = -sin（x）, 0 结论：求导运算是对矩阵元素逐个进行的【例2.3-3】求的Jacobian矩阵syms x1 x2;f=x1*exp（x2）;x2;cos（x1）*sin（x2）;v=x1 x2;fjac=jacobian（f,v） fjac = exp（x2）, x1*exp（x2） 0, 1 -sin（x1）*sin（x2）, cos（x1）*cos（x2） 【例2.3-4】，求需要分区间求导（1）对于x0

9、,据导数定义求导syms d positivef_p=sin（x）; %x0时，f（x）改写而成df_p=limit（subs（f_p,x,x+d）-f_p）/d,d,0） %求x0区间的 导数 df_p0=limit（subs（f_p,x,d）-subs（f_p,x,0）/d,d,0） %x=0+的右导数 df_p =cos（x）df_p0 =1 （2）对于x0,据导数定义求导f_n=sin（-x）; % x0时，f（x）改写而成 df_n=limit（f_n-subs（f_n,x,x-d）/d,d,0） %求x0区间的导数 df_n0=limit（subs（f_n,x,0）-subs（f_

10、n,x,-d）/d,d,0） %x=0-的左导数df_n =-cos（x）df_n0 =-1 （3）直接利用diff求导f=sin（abs（x）;dfdx=diff（f,x） %只能得到时的导函数 dfdx =cos（abs（x）*abs（1,x） abs（1,x）表示abs（x）的一阶导数，只在时有值，等价于sign（x）.用abs（1,sym（x）可观察（4）图形观察xn=-3/2*pi:pi/50:0;xp=0:3/2*pi;xnp=xn,xp（2:end）;hold onplot（xnp,subs（f,x,xnp）,k,LineWidth,2.5） % plot（xn,subs（df_

11、n,x,xn）,-r,2.5） plot（xp,subs（df_p,x,xp）,:r,2.5）legend（char（f）,char（df_n）,char（df_p）,LocationNorthEast）% grid onxlabel（hold off 图 2.3-1 函数及其导函数【例2.3-6】求在处展开的8阶Maclaurin级数。r=taylor（x*exp（x）,9,x,0） %忽略9阶及9阶以上小量的展开 r =x+x2+1/2*x3+1/6*x4+1/24*x5+1/120*x6+1/720*x7+1/5040*x8 【例2.3-8】求syms k t;f1=t k3;f2=1/

12、（2*k-1）2,（-1）k/k; findsym（f1,1）findsym（f2）s1=simple（symsum（f1） %f1的自变量被确认为ts2=simple（symsum（f2,1,inf） %f2的自变量被确认为ktks1 = 1/2*t*（t-1）, k3*ts2 = 1/8*pi2, -log（2） 说明 通式中的自变量只取整数值 求和指令中的f可以是符号矩阵，此时求和操作将对矩阵中的“元素通式”逐个进行，但通式矩阵的自变量及其取值区间对各“元素通式”是相同的。 由于findsym对f符号矩阵的自变量进行确认时，不是对该矩阵元素分别进行的，而是把矩阵f作为一个统一的对象处理的

13、。所以当f矩阵中含有一个以上非常数符号变量时要特别注意，求和对哪个符号变量进行。【例2.3-9】求f=sqrt（1+x）/x）/xs=int（f,x） %求不定积分s=simple（simple（s） %简化计算结果（1+x）/x）（1/2）/xs =（1+x）/x）（1/2）/x*（-2*（x2+x）（3/2）+2*（x2+x）（1/2）*x2+log（1/2+x+（x2+x）（1/2）*x2）/（1+x）*x）（1/2）log（1/2+x+（1+x）*x）（1/2）-2*（1+x）*x）（1/2）/x 【例2.3-10】求syms a b x;f=a*x,b*x2;1/x,sin（x）;d

14、isp（The integral of f ispretty（int（f） %求符号函数矩阵The integral of f is 2 3 1/2 a x 1/3 b x log（x） -cos（x） pretty指令把通常在一个物理行里显示的表达式扩展成用两个物理行显示的表达式。目的是读起来更方便。【例2.3-11】求积分syms x y zF2=int（int（int（x2+y2+z2,z,sqrt（x*y）,x2*y）,y,sqrt（x）,x2）,x,1,2）VF2=vpa（F2） %积分结果用32位数字表示F2 =1610027357/6563700+64/225*2（3/4）-60

15、72064/348075*2（1/2）+14912/4641*2（1/4）VF2 =224.92153573331143159790710032805 对于内积分上下限为函数的多重积分，采用数值积分时编程较复杂。【例2.3-12】求阿基米德（Archimedes）螺线到间的曲线长度函数，并求出时的曲线长度。syms a r theta phi positive %限定符号参数和变量取正值x=r*cos（theta）;x=subs（x,r,a*theta）;y=r*sin（theta）;y=subs（y,r,a*theta）;dLdth=sqrt（diff（x,theta）2+diff（y,th

16、eta）2）;L=simple（int（dLdth,theta,0,phi） L =1/2*a*（phi*（phi2+1）（1/2）-log（-phi+（phi2+1）（1/2） L_2pi=subs（L,a,phi,sym（1,2*pi） %获得完全准确值L_2pi_vpa=vpa（L_2pi） %计算32精度近似值L_2pi =pi*（1+4*pi2）（1/2）+1/2*asinh（2*pi）L_2pi_vpa =21.256294148209098800702512272566 L1=subs（L,a,sym（1）;ezplot（L1*cos（phi）,L1*sin（phi）,0,2*p

17、i） %ezplot的参变量使用格式grid on %给坐标纸打方格线hold on %保持上述图形窗，在其中继续绘图x1=subs（x,a,sym（ %使参数a数字化 y1=subs（y,a,sym（h1=ezplot（x1,y1,0,2*pi）; %为改变ezplot绘制曲线色彩必须借助“图柄”set（h1,Color,5） %通过图柄改变曲线图形对象属性legend（螺线长度-幅角曲线阿基米德螺线hold off %在上述图形窗中，不再允许画任何图形图 2.3-2 阿基米德螺线（粗红）和螺线长度函数（细蓝）【例2.4-1】求S=dsolve（Dx=y,Dy=-x） 从的输出，只能看到个域

18、.x和S.ydisp（blanks（12）,blanks（21）,）,disp（S.x,S.y） S = x: 1x1 sym y: x y -C1*cos（t）+C2*sin（t）, C1*sin（t）+C2*cos（t） 【例2.4-2】图示微分方程的通解和奇解的关系。y=dsolve（Dy）2-x*Dy+y=0） % 注意书写规则，指定独立变量clf,hold on,ezplot（y（1）,-6,6,-4,8,1） %画奇解cc=get（gca,Children %取奇解曲线的图柄set（cc,5） %把奇解画成粗红线for k=-2:0.5:2;ezplot（subs（y（2）,C1,

19、k）,-6,6,-4,8,1）;end, %画通解，ezplot不能画多条曲线，必须循环解决hold off,title（fontname隶书fontsize16通解和奇解y = 1/4*x2 -C12+x*C1 dsolve可求通解和特解。 若按题写成y=dsolve（y=x*Dy-（Dy）2）出错图 2.4-1 通解和奇解曲线【例 2.5-1】求的Fourier变换。（1）求Fourier变换syms t wut=heaviside（t）; %MATLAB7以前版本采用sym（Heaviside（t-a）UT=fourier（ut） UT =pi*dirac（w）-i/w （2）求Four

20、ier反变换Ut=ifourier（UT,w,t） %结果应为原函数Ut =heaviside（t） 单位阶跃函数heaviside（t-a）=1,当ta;heaviside（t-a）=0,当ta。【例2.5-2】根据Fourier变换定义，用积分指令求方波脉冲syms A t wsyms tao positive %限定是必须的yt=heaviside（t+tao/2）-heaviside（t-tao/2）;Yw=fourier（A*yt,t,w） Yw =2*A/w*sin（1/2*tao*w） Yt=ifourier（Yw,w,t） Yt =A*（heaviside（t+1/2*tao）

21、-heaviside（t-1/2*tao） （3）绘制曲线yt3=subs（yt,tao,3）Yw3=subs（Yw,A,tao,1,3）subplot（2,1,1）Ht=ezplot（yt3,-3,3）;set（Ht,3） subplot（2,1,2）,ezplot（Yw3） yt3 =heaviside（t+3/2）-heaviside（t-3/2）Yw3 =2/w*sin（3/2*w） 图 2.5-1 时域方波及其Fourier变换【例2.5-3】求的Fourier变换，在此是参数，是时间变量。syms t x w;ft=exp（-（t-x）*heaviside（t-x）;F1=simp

22、le（fourier（ft,t,w）%给出以w为频率变量的正确结果F2=simple（fourier（ft） %误把x当作时间变量F3=simple（fourier（ft,t） %误把x当作时间变量，又误把 %t当作频率变量F1 =1/（1+i*w）/exp（i*x*w）i*exp（-i*t*w）/（i+w）F3 =-exp（-t*（2+i*t）/（-1+i*t） 【例2.5-4】求的Laplace变换。本例演示：符号参数的属性限定；diric, heaviside的调用；Laplace变换是对矩阵元素逐个实施的。%positive（1）在不对a,b限定时，运行结果可能不正确。（2）对a，b进

23、行限定后再运行可得正确结果（3）再去除限定，运行仍然正确。但 当利用指令maple restart重新启动Maple内核时，再在没限定情况下运行，仍将出错。maple restartsyms t s;syms a b %positiveDt=dirac（t-a）; Ut=heaviside（t-b）; Mt=Dt,Ut;exp（-a*t）*sin（b*t）,t2*exp（-t）;MS=laplace（Mt,t,s） 说明 若不对符号参数 a, b进行限定，不能确保heaviside（t-b）被正确变换。 dirac（t-a） 函数的数学含义是在7.0版以前，它采用sym（Dirac（t-a）实现。比较而言，dirac（t-a） 更便于使用。【例2.5-6】求序列 的Z变换，并用反变换验算。离散单位函数的构成；ztrans, iztrans的用法。（1）单位函数及性质验证syms nDelta=sym（charfcn0（n） %定义单位函数 D0=subs（Delta,n,0）; %计算D15=subs（Delta,n,15）; %计算D0,D15disp（D0,D15） （2）求序列的Z变换syms zfn=2*Delta+6*（1-（1/2）n） %借助自定义的