第二章-Matlab软件的数值计算方法2

二、函数求极值1.fmin的命令格式X=fmin(‘Fun’,x1,x2)举例352.fmins的命令格式X=fmins(‘Fun’,x0)举例36三、函数求零点命令格式：

X=fzero(‘Fun’,x0)举例37举例39表2－6数值分析的部分MATLAB命令函数数值分析的M|ATLAB命令函数命令函数的格式描述fplot(‘fname’,[lb

ub])绘出上下限之间的函数fmin(‘fname’,[lb

ub])寻找上下限内的标量最小值fmins(‘fname’,xo)寻找xo附近的向量最小值fminbnd由一有范围限制的变量找出函数的最小值fminsearch由几个变量找出函数的最小值fzero(‘fname’,xo)寻找xo附近的标量函数的零点(解)frapz(x,y)给定数据点x和y,计算y=f(x)下的梯形面积积分diff(x)数组元素间的差分[t,y]=ode23(‘fname’,to,tf,yo)用2阶/3阶龙格-库塔算法解微分方程组[t,y]=ode45(‘fname’,to,tf,yo)补充：符号计算！1、符号对象和符号表达式 在MATLAB中，数值和数值变量用于数值的存储和各种数值计算。而符号常数、符号变量、符号函数、符号操作等则是用来形成符号表达式，严格按照代数、微积分等课程中的规则公式进行运算，并尽可能给出解析表达结果。1.1符号对象的生成和使用 数值计算中，数值表达式所用的变量必须事先被赋过值，否则无法计算。SymbolicMathToolbox2.1版沿用这种模式，规定：在进行符号运算时， 首先要定义基本的符号对象（可以是常数、变量、表达式），然后利用这些基本符号对象去构成新的表达式，进而从事所需的符号运算。 运算中，凡是由包含符号对象的表达式所生成的衍生对象也都是符号对象。定义基本符号对象的指令有两个：sym,syms。f=sym(arg)f=sym(argn,flagn)argv=sym(‘argv’,flagv)syms(‘argv1’,’argv2’,’argvk’)symsargv1argv2argvk【例1】符号常数形成中的差异a1=[1/3,pi/7,sqrt(5),pi+sqrt(5)] a2=sym([1/3,pi/7,sqrt(5),pi+sqrt(5)]) a3=sym('[1/3,pi/7,sqrt(5),pi+sqrt(5)]')a23=a2-a3a1=0.33330.44882.23615.3777a2=[1/3,pi/7,sqrt(5),6054707603575008*2^(-50)]a3=[1/3,pi/7,sqrt(5),pi+sqrt(5)]a23=[0,0,0,189209612611719/35184372088832-pi-5^(1/2)]【例2】把字符表达式转换为符号变量y=sym('2*sin(x)*cos(x)')y=simple(y) y=2*sin(x)*cos(x)y=sin(2*x)说明：符号运算中，如事先没对表达式中的独立符号变量进行定义，软件会自动检查哪些字符是软件的符号函数，哪些是变量符号，且总把在英文字母表中离x最近的字母认为独立符号变量。【例3】用符号计算验证三角等式。symsfai1fai2y=simple(sin(fai1)*cos(fai2)-cos(fai1)*sin(fai2))y=sin(fai1-fai2)【例4】求矩阵的行列式值、逆和特征根symsa11a12a21a22A=[a11,a12;a21,a22]DA=det(A),IA=inv(A),EA=eig(A)A=[a11,a12][a21,a22]DA=a11*a22-a12*a21IA=[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)]1.2符号计算中的算符和基本函数 由于matlab中采用了重载技术，使得用来构成符号计算表达式的算符和基本函数，无论在形状、名称上，还是在使用方法上，都与数值计算中的算符和基本函数几乎完全相同。

1）基本运算符

+-*\/矩阵

.*./.\元素对元素

‘.’矩阵

2）关系运算符 符号对象中没有大于、大于等于、小于、小于等于的概念，而只有是否等于的概念。 算符＝＝、～＝分别对算符两边的对象进行相等、不等的比较。当事实为真时，比较结果用1表示，假用0表示。

3）三角函数、双曲函数及它们的反函数 除atan2只用于数值计算外，其余的三角函数（如sin）、双曲函数（如cosh）以及它们的反函数（如asin,acosh）在数值计算和符号计算中使用方法都相同。

4）指数、对数函数 数值、符号计算中，函数sqrt,exp,expm的使用方法完全相同，至于对数函数，符号计算中只有自然对数log，而没有数值计算中的log2,log10。

5）复数函数

conj,real,imag,abs相同，但无求相角的指令。

6）矩阵代数指令

diag,triu,tril,inv,det,rank,rref,null,colspace,poly,expm,eig基本与数值计算相同，只有svd稍微不同。

1.3符号表达式中自由变量的确定

findsym可实现对表达式中所有自由符号变量或指定数目的独立自变量的自动认定。

findsym(expr)

findsym(expr,n)【例7】对独立自由符号变量的自动辨认。symsabxXY;k=sym('3');z=sym('c*sqrt(delta)+y*sin(theta)');EXPR=a*z*X+(b*x^2+k)*Y;findsym(EXPR) ans=X,Y,a,b,c,delta,theta,x,yfindsym(EXPR,1)ans=xfindsym(EXPR,2),findsym(EXPR,3)ans=x,yans=x,y,theta

【例8】findsym确定自由变量是对整个矩阵进行的。symsabtuvxy;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符号对象的操作和转换

2.1符号表达式的操作

collect expand factor

horner

numden simplify pretty

最常用的式simple(expr)【例9】简化symsxf=(1/x^3+6/x^2+12/x+8)^(1/3);g1=simple(f)g2=simple(g1)g1=(2*x+1)/xg2=2+1/x

2.2置换操作

一、子表达式置换操作

[RS,ssub]=subexpr(S,ssub)【例10】把复杂表达式中所含的多个相同子表达式用一个符号代替，使表达简洁。clearallsymsabcdW[V,D]=eig([a

b;cd]);[RVD,W]=subexpr([V;D],W)RVD=[-(1/2*d-1/2*a-1/2*W)/c,-(1/2*d-1/2*a+1/2*W)/c][1,1][1/2*d+1/2*a+1/2*W,0][0,1/2*d+1/2*a-1/2*W]W=(d^2-2*a*d+a^2+4*b*c)^(1/2)二、通用置换指令

RES=subs(ES,old,new) RES=subs(ES,new)【例11】用简单算例演示subs的置换规则。symsax;f=a*sin(x)+5;f1=subs(f,'sin(x)',sym('y'))f1=a*y+5f2=subs(f,{a,x},{2,sym(pi/3)}) f2=3^(1/2)+5f3=subs(f,{a,x},{2,pi/3}) f3=6.7321f4=subs(subs(f,a,2),x,0:pi/6:pi) f4=5.00006.00006.73217.00006.73216.00005.0000f5=subs(f,{a,x},{0:6,0:pi/6:pi}) f5=5.00005.50006.73218.00008.46417.50005.0000

2.3符号数值精度控制和任意精度计算 数值计算受字长的限制，每次数值操作都带截断误差，任何数值计算不管什么算法都将产生积累误差。Matlab中每个算术操作结果的相对精度约为16位数字。但符号计算的结果是绝对准确的，不包含任何计算误差。

double(x)

digits(n)

xs=vpa(x)

xs=vpa(x,n)【例12】指令使用演示。p0=sym('(1+sqrt(5))/2');p1=sym((1+sqrt(5))/2) pd=double(p0) e01=vpa(abs(p0-p1)) e0d=vpa(abs(p0-pd)) p1=7286977268806824*2^(-52)pd=1.6180e01=.543211520368250e-16e0d=.543211520368250e-16p2=vpa(p0) e02=vpa(abs(p0-p2),64) p2=1.6180339887498948482045868343657

e02=.61882279690820194237137864551377e-31digits Digits=32

2.4符号对象与其他数据对象间的转换 数值、符号、字符是matlab中的三种不同的数据类型。Matlab为每种数据类型提供了各自特定的生成指令和操作指令。 转换：vpadoublevpasymdoublestr2num!3符号微积分3.1符号序列的求和对数学上的通式求和：S=symsum(f,v,a,b)【例13】求symskt;f1=[tk^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)]3.2符号微分dfdvn=diff(f,v,n)【例14】求、和symsatx;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]3.3符号积分 积分有不定积分、定积分和重积分等。一般来说，积分比微分更难求取。

intf=int(f,v)

intf=int(f,v,a,b)【例15】求。symsabx;f=[a*x,b*x^2;1/x,sin(x)];disp('Theintegraloffis');pretty(int(f))Theintegraloffis[23][1/2ax1/3bx][][log(x)-cos(x)] 【例16】求定积分，注意内积分上下限都是函数。symsxyzF2=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=1610027357/6563700-6072064/348075*2^(1/2)+14912/4641*2^(1/4)+64/225*2^(3/4)VF2=224.92153573331143159790710032805！4符号代数方程的求解4.1线性方程组的符号解

矩阵计算是求解线性方程组最简便有效的方法。不管数据对象是数值还是符号，实现矩阵运算的指令形式几乎完全相同。【例17】求线性方程组的解。A=sym([11/21/2-1;11-11;1-1/4-11;-8-111]);b=sym([0;10;0;1]);X1=A\bX1=[1][8][8][9]4.2一般代数方程组的解 这里所讲的一般代数方程包括线性、非线性等，求解指令是solve。当方程组不存在符号解又无其他自由参数时，solve将给出数值解。

s=solve(‘eq1’,’eq2’,…,’eqn’,’v1’,’v2’,…,’vn’) s=solve(exp1,exp2,…,expn,v1,v2,…,vn)【例18】求方程组，关于的解。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:[2x1sym]z:[2x1sym]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))]【例19】solve指令求构成的“欠定”方程组解。symsdnpqeq1=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);S.d,S.n,S.p,S.q

ans=dans=8ans=4*d+4ans=3*d+6【例20】求的解。clearall,syms

x;s=solve('(x+2)^x=2','x')s=.69829942170241042826920133106081!5符号微分方程的求解S=dsolve(‘a_1’,’a_2’,…,’a_n’)【例21】求的解。S=dsolve('Dx=y,Dy=-x');disp([blanks(12),'x',blanks(21),'y']),disp([S.x,S.y])xy[cos(t)*C1+sin(t)*C2,-sin(t)*C1+cos(t)*C2]四、数值积分1、符号解析法2、利用quad和quadl求积分值 命令格式1:quad(‘fun’,A,B)

命令格式2：quad(‘fun’,A,B,TOL)

命令格式3：quadl(‘fun’,A,B)

命令格式4：quadl(‘fun’,A,B,TOL)五、方程（组）求解1、solve的命令格式 命令格式1：g=solve(‘eqn’)

命令格式2：g=solve(‘eqn’,var)2、其他命令格式 命令格式1：g=solve(‘eqn1’,’eqn2’,…,’eqnn’)

命令格式2：g=solve(‘eqn1’,’eqn2’,…,’eqnn’,’var1,var2,…,varn’)六、求函数极值命令格式1：limit(F,x,a)命令格式2：limit(F,a)命令格式3：limit(F)命令格式4：limit(F,x,a,’right’)命令格式5：limit(F,x,a,’left’)七、Taylor级数分析方法命令格式1：r=taylor(f,n,v)命令格式2：r=taylor(f)命令格式3：r=taylor(f,n,v,a)八、傅立叶及其他重要变换1、fft的命令格式命令格式1：fft(X)命令格式2：fft(X,N)命令格式3：fft(X,[],DIM)或者fft(X,N,DIM)举例：

e_fft.mclear,randn('state',0)t=linspace(0,10,512);y=3*sin(5*t)-6*cos(9*t)+5*randn(size(t));plot(t,y)Y=fft(y);Ts=t(2)-t(1)Ws=2*pi/Ts;Wn=Ws/2w=linspace(0,Wn,length(t)/2);Ya=abs(Y(1:length(t)/2));subplot(1,2,1);plot(w,Ya)ii=find(w<=20);subplot(1,2,2);plot(w(ii),Ya(ii))grid,xlabel('Frequency

Rad/s')e50.mt=0:0.001:0.6;x=2*sin(2*pi*50*t)+1*sin(2*pi*120*t)+1.5*sin(2*pi*50*t);y=x+0.5*randn(size(t));subplot(2,1,1)plot(y(1:50));title('传感器叠加有随机噪声信号的输出信号');Y=fft(y,512);f=1000*(0:256)/512;subplot(2,1,2)plot(f,Y(1:257));title('传感器叠加有随机噪声信号的输出信号的频域特性');>>randn(3)ans=1.55190.44270.20181.38360.91110.7629-0.7581-1.0741-1.2882

e51.mt=0.5:0.001:1.5;x=t-0.5;subplot(2,1,1);plot(x(1:50));title('时间函数信号');X=fft(x,512);f=1000*(0:256)/512;subplot(2,1,2);plot(f,X(1:257));title('时间函数信号的频域特性');将函数用傅立叶级数展开

e52.mfunction[a0,ak,bk]=fourier_s(f)symsxna0=int(f,-pi,pi)/pi;an=int(f*cos(n*x),-pi,pi)/pi;an=simple(an)bn=int(f*sin(n*x),-pi,pi)/pi;bn=simple(bn)symsxf=x+2;[a0,ak,bk]=fourier_s(f)>>an=4/n*sin(pi*n)/pi

bn=2/n^2/pi*sin(pi*n)-2/n*cos(pi*n)a0=42、Fourier变换 命令格式1：F=fourier(f)

命令格式2：F=fourier(f,v)

命令格式3：F=fourier(f,u,v)

e53.msymsxwuvf=sin(x)*exp(-x^2);F1=fourier(f)g=log(abs(w));F2=fourier(g)h=x*exp(-abs(x));F3=fourier(h,u)symsxrealk=cosh(-x^2*abs(v))*sinh(u)/v;F4=fourier(k,v,u)>>F1=-1/2*i*pi^(1/2)*exp(-1/4*(w-1)^2)+1/2*i*pi^(1/2)*exp(-1/4*(w+1)^2)

F2=fourier(log(abs(w)),w,t)

F3=-4*i/(1+u^2)^2*u

F4=sinh(u)*(1/2*fourier(1/v*exp(x^2*abs(v)),v,u)-i*atan(u/x^2))3、Fourier逆变换 命令格式1：f=ifourier(F)

命令格式2：f=ifourier(F,u)

命令格式3：f=ifourier(F,v,u)

e54.msymswvxtsymsarealf=sqrt(exp(-w^2/(4*a^2)));IF1=ifourier(f)g=exp(-abs(x));IF2=ifourier(g)h=sinh(-abs(w))-1;IF3=simple(ifourier(h,t))symswrealk=exp(-w^2*abs(v))*sin(v)/v;IF4=ifourier(k,v,t)>>IF1=1/2*8^(1/2)*(pi*a^2)^(1/2)*exp(-2*x^2*a^2)/pi

IF2=1/(1+t^2)/pi

IF3=-1/2*(pi*ifourier(exp(abs(w)),w,t)+pi*ifourier(exp(abs(w)),w,t)*t^2-1+2*pi*Dirac(t))/(1+t^2)/pi

IF4=-1/2*(-atan((t+1)/w^2)+atan((t-1)/w^2))/pi4、Laplace变换 命令格式1：L=laplace(F)

命令格式2：laplace(F,t)

命令格式3：laplace(F,w,z)

e55.msymsaxstvf1=sqrt(t);L1=laplace(f1)f2=1/sqrt(s);L2=laplace(f2)f3=exp(-a*t);L3=laplace(f3,x)f4=1-sin(t*v);L4=laplace(f4,v,x)>>L1=1/2*pi^(1/2)/s^(3/2)

L2=(pi/t)^(1/2)

L3=1/(x+a)

L4=1/x-t/(x^2+t^2)5、Laplace逆变换 命令格式1：F=ilaplace(L)

命令格式2：F=ilaplace(L,y)

命令格式3：F=ilaplace(L,