第二章-Matlab軟件的數(shù)值計算方法2

二、函數(shù)求極值1.fmin的命令格式X=fmin(‘Fun’,x1,x2)舉例352.fmins的命令格式X=fmins(‘Fun’,x0)舉例36三、函數(shù)求零點命令格式：

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

ub])繪出上下限之間的函數(shù)fmin(‘fname’,[lb

ub])尋找上下限內(nèi)的標(biāo)量最小值fmins(‘fname’,xo)尋找xo附近的向量最小值fminbnd由一有范圍限制的變量找出函數(shù)的最小值fminsearch由幾個變量找出函數(shù)的最小值fzero(‘fname’,xo)尋找xo附近的標(biāo)量函數(shù)的零點(解)frapz(x,y)給定數(shù)據(jù)點x和y,計算y=f(x)下的梯形面積積分diff(x)數(shù)組元素間的差分[t,y]=ode23(‘fname’,to,tf,yo)用2階/3階龍格-庫塔算法解微分方程組[t,y]=ode45(‘fname’,to,tf,yo)補(bǔ)充：符號計算！1、符號對象和符號表達(dá)式 在MATLAB中，數(shù)值和數(shù)值變量用于數(shù)值的存儲和各種數(shù)值計算。而符號常數(shù)、符號變量、符號函數(shù)、符號操作等則是用來形成符號表達(dá)式，嚴(yán)格按照代數(shù)、微積分等課程中的規(guī)則公式進(jìn)行運(yùn)算，并盡可能給出解析表達(dá)結(jié)果。1.1符號對象的生成和使用 數(shù)值計算中，數(shù)值表達(dá)式所用的變量必須事先被賦過值，否則無法計算。SymbolicMathToolbox2.1版沿用這種模式，規(guī)定：在進(jìn)行符號運(yùn)算時， 首先要定義基本的符號對象（可以是常數(shù)、變量、表達(dá)式），然后利用這些基本符號對象去構(gòu)成新的表達(dá)式，進(jìn)而從事所需的符號運(yùn)算。 運(yùn)算中，凡是由包含符號對象的表達(dá)式所生成的衍生對象也都是符號對象。定義基本符號對象的指令有兩個：sym,syms。f=sym(arg)f=sym(argn,flagn)argv=sym(‘a(chǎn)rgv’,flagv)syms(‘a(chǎn)rgv1’,’argv2’,’argvk’)symsargv1argv2argvk【例1】符號常數(shù)形成中的差異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】把字符表達(dá)式轉(zhuǎn)換為符號變量y=sym('2*sin(x)*cos(x)')y=simple(y) y=2*sin(x)*cos(x)y=sin(2*x)說明：符號運(yùn)算中，如事先沒對表達(dá)式中的獨立符號變量進(jìn)行定義，軟件會自動檢查哪些字符是軟件的符號函數(shù)，哪些是變量符號，且總把在英文字母表中離x最近的字母認(rèn)為獨立符號變量。【例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符號計算中的算符和基本函數(shù) 由于matlab中采用了重載技術(shù)，使得用來構(gòu)成符號計算表達(dá)式的算符和基本函數(shù)，無論在形狀、名稱上，還是在使用方法上，都與數(shù)值計算中的算符和基本函數(shù)幾乎完全相同。

1）基本運(yùn)算符

+-*\/矩陣

.*./.\元素對元素

‘.’矩陣

2）關(guān)系運(yùn)算符 符號對象中沒有大于、大于等于、小于、小于等于的概念，而只有是否等于的概念。 算符＝＝、～＝分別對算符兩邊的對象進(jìn)行相等、不等的比較。當(dāng)事實為真時，比較結(jié)果用1表示，假用0表示。

3）三角函數(shù)、雙曲函數(shù)及它們的反函數(shù) 除atan2只用于數(shù)值計算外，其余的三角函數(shù)（如sin）、雙曲函數(shù)（如cosh）以及它們的反函數(shù)（如asin,acosh）在數(shù)值計算和符號計算中使用方法都相同。

4）指數(shù)、對數(shù)函數(shù) 數(shù)值、符號計算中，函數(shù)sqrt,exp,expm的使用方法完全相同，至于對數(shù)函數(shù)，符號計算中只有自然對數(shù)log，而沒有數(shù)值計算中的log2,log10。

5）復(fù)數(shù)函數(shù)

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

6）矩陣代數(shù)指令

diag,triu,tril,inv,det,rank,rref,null,colspace,poly,expm,eig基本與數(shù)值計算相同，只有svd稍微不同。

1.3符號表達(dá)式中自由變量的確定

findsym可實現(xiàn)對表達(dá)式中所有自由符號變量或指定數(shù)目的獨立自變量的自動認(rèn)定。

findsym(expr)

findsym(expr,n)【例7】對獨立自由符號變量的自動辨認(rèn)。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確定自由變量是對整個矩陣進(jìn)行的。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符號對象的操作和轉(zhuǎn)換

2.1符號表達(dá)式的操作

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置換操作

一、子表達(dá)式置換操作

[RS,ssub]=subexpr(S,ssub)【例10】把復(fù)雜表達(dá)式中所含的多個相同子表達(dá)式用一個符號代替，使表達(dá)簡潔。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的置換規(guī)則。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符號數(shù)值精度控制和任意精度計算 數(shù)值計算受字長的限制，每次數(shù)值操作都帶截斷誤差，任何數(shù)值計算不管什么算法都將產(chǎn)生積累誤差。Matlab中每個算術(shù)操作結(jié)果的相對精度約為16位數(shù)字。但符號計算的結(jié)果是絕對準(zhǔn)確的，不包含任何計算誤差。

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符號對象與其他數(shù)據(jù)對象間的轉(zhuǎn)換 數(shù)值、符號、字符是matlab中的三種不同的數(shù)據(jù)類型。Matlab為每種數(shù)據(jù)類型提供了各自特定的生成指令和操作指令。 轉(zhuǎn)換：vpadoublevpasymdoublestr2num!3符號微積分3.1符號序列的求和對數(shù)學(xué)上的通式求和：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】求定積分，注意內(nèi)積分上下限都是函數(shù)。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符號代數(shù)方程的求解4.1線性方程組的符號解

矩陣計算是求解線性方程組最簡便有效的方法。不管數(shù)據(jù)對象是數(shù)值還是符號，實現(xiàn)矩陣運(yùn)算的指令形式幾乎完全相同?！纠?7】求線性方程組的解。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一般代數(shù)方程組的解 這里所講的一般代數(shù)方程包括線性、非線性等，求解指令是solve。當(dāng)方程組不存在符號解又無其他自由參數(shù)時，solve將給出數(shù)值解。

s=solve(‘eq1’,’eq2’,…,’eqn’,’v1’,’v2’,…,’vn’) s=solve(exp1,exp2,…,expn,v1,v2,…,vn)【例18】求方程組，關(guān)于的解。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指令求構(gòu)成的“欠定”方程組解。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(chǎn)_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]四、數(shù)值積分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’)六、求函數(shù)極值命令格式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級數(shù)分析方法命令格式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('傳感器疊加有隨機(jī)噪聲信號的輸出信號');Y=fft(y,512);f=1000*(0:256)/512;subplot(2,1,2)plot(f,Y(1:257));title('傳感器疊加有隨機(jī)噪聲信號的輸出信號的頻域特性');>>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('時間函數(shù)信號');X=fft(x,512);f=1000*(0:256)/512;subplot(2,1,2);plot(f,X(1:257));title('時間函數(shù)信號的頻域特性');將函數(shù)用傅立葉級數(shù)展開

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,