MATLAB上機(jī)實(shí)驗(yàn)報(bào)告
平頂山學(xué)院計(jì)算機(jī)語言類課程實(shí)驗(yàn)報(bào)告(4)
課程名稱院系學(xué)號(hào)實(shí)驗(yàn)日期MATLAB語言及應(yīng)用電氣信息工程學(xué)院201*/5/3實(shí)驗(yàn)機(jī)房專業(yè)電氣工程及其自動(dòng)化姓名任課教師王凱實(shí)驗(yàn)學(xué)時(shí)23305班級(jí)機(jī)器號(hào)實(shí)驗(yàn)成績電氣二班A4一.實(shí)驗(yàn)名稱:MATLAB求解線性方程組和矩陣的初等計(jì)算二.實(shí)驗(yàn)?zāi)康暮鸵?、掌握利用MATLAB程序編輯器編寫應(yīng)用程序的方法;2、掌握MATLAB求解線性方程組的方法;3、掌握MATLAB進(jìn)行矩陣的初等計(jì)算的方法三.實(shí)驗(yàn)內(nèi)容1、給出一個(gè)信號(hào)t=0:0.001:3,其正弦信號(hào)頻率響應(yīng)特性u(píng)=sin(300t)+2cos(200t),求其幅頻特性曲線。2、已知某系統(tǒng)的閉環(huán)傳遞函數(shù)為(s)根,并判別其穩(wěn)定性。3、已知某系統(tǒng)的開環(huán)傳遞函數(shù)為G(s)圖,并判斷系統(tǒng)的穩(wěn)定性。4、已知某系統(tǒng)的開環(huán)傳遞函數(shù)為G(s)5,計(jì)算該系統(tǒng)的相位裕度和幅值裕度。s(s4)(s5)1,畫出該系統(tǒng)的開環(huán)Nyquist圖和Bodes(3s1)(4s1)s4試求該系統(tǒng)的特征5432s10s20s30s40s5025、已知某系統(tǒng)的閉環(huán)傳遞函數(shù)為(s)2,設(shè)ζ=0.707,ωn=1:1:5,畫出該2s2wnsnn系統(tǒng)的單位階躍響應(yīng)曲線。
四.實(shí)驗(yàn)設(shè)計(jì)方案(實(shí)驗(yàn)步驟或開發(fā)過程)1、clearallcloseallclct=0:0.001:3;u=sin(300*t)+2*cos(200*t);y=fft(u);yy=abs(y);plot(t,yy)2、d=[11020304154];r=roots(d)3、clearallcloseallnum=1;den=[12710];bode(num,den)figure(2)nyquist(num,den)w=logspace(0,4,50);bode(sys,w);grid;[Gm,Pm,Wg,Wc]=margin(sys)5、i=0;t=0:0.1:10;forwn=1:5i=i+1;zeta=[0.707];num=[wn^2];den=[1,2*zeta*wn,wn^2];sys=tf(num,den);y(:,i)=step(sys,t);endplot(t,y(:,1),t,y(:,2),t,y(:,3),t,y(:,4),t,y(:,5));gridxlabel("time(s)")ylabel("stepresponse")4、num=5;den=conv([1,0],conv([1,4],[1,5]));sys=tf(num,den);五.實(shí)驗(yàn)中存在問題及解決辦法這一次上機(jī)實(shí)驗(yàn)遇到很多問題,首先關(guān)于用MATLAB求自動(dòng)控制相關(guān)問題以前接觸到一點(diǎn),但是現(xiàn)在完全不會(huì),相關(guān)函數(shù)都不清楚,在上網(wǎng)查閱資料,仿照網(wǎng)友給出的程序最終得除了正確結(jié)果,但是對(duì)于怎么用程序函數(shù)判斷系統(tǒng)的穩(wěn)定性還是沒有掌握。最大的問題還是對(duì)于MATLAB軟件不熟悉。
六.實(shí)驗(yàn)結(jié)果1、5、25001.4201*1.21stepresponse15000.810000.60.45000.201*.511.522.5302、r=-7.8739-1.4201+0.9942i-1.4201-0.9942i0.3570+1.4679i0.3570-1.4679i由此看出系統(tǒng)不穩(wěn)定3、NyquistDiagram10080Magnitude(dB)012345time(s)678910BodeDiagram506040ImaginaryAxis0-50200-20-100-90-135-180-225-270-60-80-100-7-6-5-4RealAxis-3-2-10Phase(deg)-4010-210-110Frequency(rad/s)01014、Gm=36(幅值裕度)Pm=83.5835(相角裕度)
七.附錄(源程序清單)1、clearallcloseallclct=0:0.001:3;u=sin(300*t)+2*cos(200*t);y=fft(u);yy=abs(y);plot(t,yy)2、d=[11020304154];r=roots(d)3、clearallcloseallnum=1;den=[12710];bode(num,den)figure(2)nyquist(num,den)4、num=5;den=conv([1,0],conv([1,4],[1,5]));sys=tf(num,den);w=logspace(0,4,50);bode(sys,w);grid;[Gm,Pm,Wg,Wc]=margin(sys)5、i=0;t=0:0.1:10;forwn=1:5i=i+1;zeta=[0.707];num=[wn^2];den=[1,2*zeta*wn,wn^2];sys=tf(num,den);y(:,i)=step(sys,t);endplot(t,y(:,1),t,y(:,2),t,y(:,3),t,y(:,4),t,y(:,5));gridxlabel("time(s)")ylabel("stepresponse")
擴(kuò)展閱讀:matlab上機(jī)實(shí)驗(yàn)報(bào)告
上海電力學(xué)院
數(shù)值計(jì)算方法上機(jī)實(shí)習(xí)
院系:專業(yè)年級(jí):學(xué)生姓名:學(xué)號(hào):指導(dǎo)教師:報(bào)告
能源與機(jī)械工程學(xué)院動(dòng)力機(jī)械及工程201*級(jí)張亞杰ys1210124014黃建雄
201*年12月26日1
數(shù)值計(jì)算方法上機(jī)實(shí)習(xí)題
xndx,1.設(shè)In05x1(1)由遞推公式In5In1解:I0=
I=0.1823forn=1:1:20,I=-5*I+1/n;1,從I0的幾個(gè)近似值出發(fā),計(jì)算I20;n105xdx=0.1823
1計(jì)算I20編輯matlab命令如下:fprintf("%.1d%.4f\\n",n,I);end結(jié)果:
(2)粗糙估計(jì)I20,用In111In,計(jì)算I0;55nx20d解:I20=05xx1使用復(fù)合中點(diǎn)公式進(jìn)行積分,相應(yīng)的matlab程序如下:
I=0;forh=0:0.001:1,m=h+0.0005;I=I+0.001*m^20/(5+m);fprintf("%.1d%.4f\\n",m,I);enddisp(I);fork=1:20,n=21-k;I=0.2*(1/n-I);fprintf("%.1d%.4f\\n",n,I);enddisp(I)結(jié)果:
程序結(jié)束時(shí)輸出兩個(gè)I值,第一個(gè)表示I20,第二個(gè)表示I0;分別為I20=0.0082I0=0.1823
(3)分析結(jié)果的可靠性及產(chǎn)生此現(xiàn)象的原因(重點(diǎn)分析原因)。
從上述計(jì)算中分析得到如果先得到I0,再從I0由遞推公式得到I20,I20結(jié)果跟精確值相比誤差很大;如果先估算I20,在從I20有遞推公式得到I0,I0的結(jié)果跟精確值相比近似相等。
原因分析:
如果從I0推I20的近似值,需要用到遞推公式In=-5In-1+1/n,I0本身結(jié)果是有誤差的;經(jīng)過遞推公式計(jì)算20次,就等于誤差被認(rèn)為的放大5的20次方倍,所以得到的I20與其精確值相差甚遠(yuǎn)。
如果從I20推I0的近似值,需要用到In-1=0.2(1/n-In),盡管I20本身有誤差,但是經(jīng)過20次運(yùn)算,其誤差縮小到原來的0.2的20次方倍,所以得到的I0與其精確值比較相近。
2.求方程e10x20的近似根,要求xk1xk5104,并比較計(jì)算量。
x(1)在[0,1]上用二分法;
Matlab程序如下:a=0;b=1;c=b-a;n=0whilec>0.0005,x=(a+b)/2;f=exp(x)+10*x-2;iff>0,b=x;c=b-a;elseiff(2)取初值x00,并用迭代xk12ex;
10采用matlab進(jìn)行迭代的程序如下:
x=0;c=1;n=0;whilec>0.0005,m=x;m=(2-exp(m))/10;c=abs(m-x);x=m;n=n+1;fprintf("%.1d%.4f%.4f\\n",n,x,c);end結(jié)果:
解得x=0.0905
(3)加速迭代的結(jié)果;
采用matlab進(jìn)行迭代的程序如下:x=0;n=0;a=0;b=1;whileabs(a-b)>0.0005,n=n+1;a=x;y=(2-exp(x))/10;z=(2-exp(y))/10;x=x-(y-x)^2/(z-2*y+x);b=x;fprintf("%.1d%.4f%.4f\\n",n,x,abs(a-b));end結(jié)果如下:
(4)取初值x00,并用牛頓迭代法;Matlab程序如下:
x=0;a=1;n=0;whileabs(a)>0.0005,n=n+1;a=(exp(x)+10*x-2)/(exp(x)+10);x=x-a;fprintf("%.1d%.4f%.4f\\n",n,x,abs(a));end運(yùn)行結(jié)果:
(5)分析絕對(duì)誤差。迭代次數(shù)1234567891011二分法X(k)0.50000.25000.12500.06250.09380.07810.08590.08980.09180.09080.0903Erroe0.50000.25000.12500.06250.03130.01560.00780.00390.00200.00100.0005代數(shù)式迭代X(k)0.10000.08950.09060.0905Erroe0.10000.01050.00120.0001加速迭代X(k)0.09050.0905Erroe0.0905牛頓迭代X(k)0.0909Erroe0.09090.00040.00000.0905我們可以看到,在運(yùn)算要求到同一精度的情況下,采用(1)的二分法運(yùn)算了11次,采用(2)的方法運(yùn)算了4次,采用(3)的加速迭代法運(yùn)算了2次,采用(4)的牛頓迭代法也需運(yùn)算2次。也就是說牛頓的迭代的收斂速度與加速迭代速度都是超線性收斂的,而簡單迭代法是線性收斂的。而二分法收斂速度較慢,所以在工程上不經(jīng)常使用。3.鋼水包使用次數(shù)多以后,鋼包的容積增大,數(shù)據(jù)如下:xy10234569.7147101581112131610.4910.5910.6010.810.610.910.76試從中找出使用次數(shù)和容積之間的關(guān)系,計(jì)算均方差。(注:增速減少,用何種模型)解:將使用次數(shù)x與體積y的關(guān)系用matlab采用如下程序繪制在二維坐標(biāo)系:x=[2345678910111213141516];y=[6.428.29.589.59.7109.939.9910.4910.5910.6010.810.610.910.76];plot(x,y,"b*-");96.428.29.589.59.939.99結(jié)果如下:
由數(shù)據(jù)點(diǎn)分布圖可知,擬合曲線y=f(x)隨著x的增加而上升,但上升速度由快到慢,當(dāng)x趨于無窮大時(shí),y趨于某個(gè)常數(shù),故曲線有一水平漸進(jìn)線。根據(jù)上述特征很容易想到用
b/x
Logistic模型來擬合該曲線。設(shè)y=f(x)的形式為y=ae(a>0,b
計(jì)算均方差s,matlab程序如下:y=[6.428.29.589.59.7109.939.9910.4910.5910.6010.810.610.910.76];s=0;forn=2:1:16,a=abs(11.679*exp(-1.1107*n.^(-1))-y(n-1));s=s+a.^2;ends=(s/15).^(1/2);disp(s);運(yùn)算結(jié)果均方差S=0.2438
小結(jié):根據(jù)已給的條件計(jì)算函數(shù)是十分困難的,但通過對(duì)離散點(diǎn)的分析及變化規(guī)律找出其中的規(guī)律,并通過計(jì)算來得到實(shí)際的函數(shù)是十分有用的方法。本題就是這樣做的一個(gè)典型,在n=1/x和m=lny的基礎(chǔ)上找到了它們之間的關(guān)系并通過這種關(guān)系來擬合原函數(shù),并最終驗(yàn)證計(jì)算結(jié)果。
041010014101050141012,b,Axb4.設(shè)A101410501014120601014分析下列迭代法的收斂性,并求xk1xk(1)JACOBI迭代;解matlab計(jì)算程序如下:A=[4-10-100;-14-10-10;0-14-10-1;-10-14-10;0-10-14-1;00-10-14];b=[0;5;-2;5;-2;6];error=1;D=diag(diag(A));L=D-tril(A);U=D-triu(A);X=zeros(size(b));whileerror>0.0001,
82104的近似解及相應(yīng)的迭代次數(shù)。X=D\\(b+L*X+U*X);error=norm(b-A*X)/norm(b);enddisp(x);disp(error);解得X=[0.9999;1.9999;0.9998;1.9999;0.9998;1.9999]error=7.0206e-05
(2)GAUSS-SEIDEL迭代;
A=[4-10-100;-14-10-10;0-14-10-1;-10-14-10;0-10-14-1;00-10-14];b=[0;5;-2;5;-2;6];error=1;D=diag(diag(A));L=D-tril(A);U=D-triu(A);X=zeros(size(b));whileerror>0.0001,X=(D-L)\\(b+U*X);error=norm(b-A*X)/norm(b);enddisp(x);disp(error);解得X=[0.9998;1.9998;0.9998;1.9999;0.9999;1.9999]error=5.5892e-05(3)SOR迭代(1.334,1.95,0.95)。
N=1.334使用matlab求解程序如下:A=[4-10-100;-14-10-10;0-14-10-1;-10-14-10;0-10-14-1;00-10-14];b=[0;5;-2;5;-2;6];error=1;D=diag(diag(A));L=D-tril(A);U=D-triu(A);X=zeros(size(b));whileerror>0.0001,n=1.334;X=(D-n*L)\\[(1-n)*D+n*U]*X+n*[(D-n*L)\\b];error=norm(b-A*X)/norm(b);disp(X);enddisp(error);此循環(huán)得到的X=[0.9999;2.0000;1.0000;1.9999;1.0000;2.0000]error=6.3632e-05
N=1.95使用matlab求解程序如下:A=[4-10-100;-14-10-10;0-14-10-1;-10-14-10;0-10-14-1;00-10-14];b=[0;5;-2;5;-2;6];error=1;D=diag(diag(A));L=D-tril(A);U=D-triu(A);X=zeros(size(b));whileerror>0.0001,n=1.95;X=(D-n*L)\\[(1-n)*D+n*U]*X+n*[(D-n*L)\\b];error=norm(b-A*X)/norm(b);disp(X);enddisp(error);此循環(huán)得到的X=[0.9999;2.0001;0.9999;1.9999;1.0001;1.9999]error=9.0363e-05N=0.95使用matlab求解程序如下:A=[4-10-100;-14-10-10;0-14-10-1;-10-14-10;0-10-14-1;00-10-14];b=[0;5;-2;5;-2;6];error=1;D=diag(diag(A));L=D-tril(A);U=D-triu(A);X=zeros(size(b));whileerror>0.0001,n=0.95;X=(D-n*L)\\[(1-n)*D+n*U]*X+n*[(D-n*L)\\b];error=norm(b-A*X)/norm(b);disp(X);enddisp(error);此循環(huán)得到的X=[0.9997;1.9997;0.9997;1.9998;0.9998;1.9999]error=8.6235e-05
6313A3215.用逆冪迭代法求最接近于11的特征值和特征向量,準(zhǔn)確到10。
111解:matlab程序如下:a=[631;321;111];I=[100;010;001];b0=a-11*I;v0=[1;1;1];m=max(abs(v0));flab=1;whileflab>0.001,u=v0/m;v0=b0\\u;[tv,ti]=max(abs(v0));n=v0(ti);flab=abs(n-m);m=n;endm=1/m+11;disp(m);運(yùn)行結(jié)果如下:
即離11最近的特征值為7.8745;相應(yīng)的特征向量u=[1.0000;0.5503;0.2271]。6.用經(jīng)典R-K方法求解初值問題(1)2y1y22sinxy1,x[0,10],
y2y12y22cosx2sinxy1(0)2;y2(0)3解:采用經(jīng)典R-K公式計(jì)算的MATLAB程序如下:
y1=2;y2=3;forh=0:0.01:10,k1=0.01*(-2*y1+y2+2*sin(h));l1=0.01*(y1-2*y2+2*cos(h)-2*sin(h));k2=0.01*(-2*(y1+0.5*k1)+(y2+0.5*l1)+2*sin(h+0.005));l2=0.01*((y1+0.5*k1)-2*(y2+0.5*l1)+2*cos(h+0.005)-2*sin(h+0.005));k3=0.01*(-2*(y1+0.5*k2)+(y2+0.5*l2)+2*sin(h+0.005));l3=0.01*((y1+0.5*k2)-2*(y2+0.5*l2)+2*cos(h+0.005)-2*sin(h+0.005));k4=0.01*(-2*(y1+k3)+(y2+l3)+2*sin(h+0.01));l4=0.01*((y1+k3)-2*(y2+l3)+2*cos(h+0.01)-2*sin(h+0.01));y1=y1+1/6*(k1+2*k2+2*k3+k4);y2=y2+1/6*(l1+2*l2+2*l3+l4);ifh==fix(h);fprintf("%.1d%.4f%.4f\\n",h,y1,y2);elseendend結(jié)果如下所示:
2y1y22sinxy1(2),x[0,10],
y998y999y999cosx999sinx122y1(x)2exsinx和精確解比較,分析結(jié)論。xy2(x)2ecosxMatlab程序如下:y1=2;y2=3;forh=0:0.00001:10,k1=0.00001*(-2*y1+y2+2*sin(h));l1=0.00001*(998*y1-999*y2+999*cos(h)-999*sin(h));y1(0)2。y(0)32k2=0.00001*(-2*(y1+0.5*k1)+(y2+0.5*l1)+2*sin(h+0.000005));l2=0.00001*(998*(y1+0.5*k1)-999*(y2+0.5*l1)+999*cos(h+0.000005)-999*sin(h+0.000005));k3=0.00001*(-2*(y1+0.5*k2)+(y2+0.5*l2)+2*sin(h+0.005));l3=0.00001*(998*(y1+0.5*k2)-999*(y2+0.5*l2)+999*cos(h+0.000005)-999*sin(h+0.000005));k4=0.00001*(-2*(y1+k3)+(y2+l3)+2*sin(h+0.00001));l4=0.00001*(998*(y1+k3)-999*(y2+l3)+999*cos(h+0.00001)-999*sin(h+0.00001));y1=y1+1/6*(k1+2*k2+2*k3+k4);y2=y2+1/6*(l1+2*l2+2*l3+l4);ifh==fix(h),fprintf("%.1d%.4f%.4f\\n",h,y1,y2);elseendend結(jié)果如下:
精確解:
forx=0:1:10,y1=2*exp(-x)+sin(x);y2=2*exp(-x)+cos(x);fprintf("%.1d%.4f%.4f\\n",x,y1,y2);end結(jié)果;
結(jié)果分析:
四階RungeKutta方法得到的結(jié)果已很接近精確解,證明這種迭代方法精確度很好,是一種有效的算法。但是要注意龍格-庫塔公式的推導(dǎo)基于泰勒展開方法,因而它要求所求的的解具有較好的光滑性質(zhì)。反之,如果解得光滑性差,那么,使用四階龍格-庫塔求得的數(shù)值解精度就不是太高,此種情況可以采用縮小步長來解決,比如上述計(jì)算。7.用有限差分法求解邊值問題(h=0.1):
y(1x2)y0.y(1)y(1)1y(1x2)y0微分方程式可以變?yōu)橛糜邢薏罘址╩atlab程序如下:
y(1)y(1)1h=0.1;n=2/0.1-1;g(1)=1/(h.^2);g(n)=1/(h^2);fori=2:1:18,
g(i)=0;endg=[g(1);g(2);g(3);g(4);g(5);g(6);g(7);g(8);g(9);g(10);g(11);g(12);g(13);g(14);g(15);g(16);g(17);g(18);g(19)];disp(g);fori=1:1:19,forj=1:1:19,ifi==1,H(1,1)=2/(h.^2)+(1+(-1+0.1*i).^2);H(1,2)=-1/(h.^2);elseifi==19,H(19,18)=-1/(h.^2);H(19,19)=2/(h.^2)+(1+(-1+0.1*i).^2);elseifj==i,H(i,j)=2/(h.^2)+(1+(-1+0.1*i).^2);elseifj==i-1,H(i,i-1)=-1/(h.^2);elseifj==i+1;H(i,i+1)=-1/(h.^2);elseH(i,j)=0;endendendenddisp(H);y=H\\g;fori=1:1:19,fprintf("%.4f%.4f\\n",-1+0.1*i,y(i,1))end
運(yùn)算結(jié)果為:
g10000000000000000000100
H矩陣為:
Y在各點(diǎn)的近似值為:
XY8.用函數(shù)y=asin(bx)擬合數(shù)據(jù).
x0.10.20.30.40.50.60.70.8y0.61.11.61.82.01.91.71.3Matlab上機(jī)程序?yàn)椋篺unction[err,a,b]=nlfitb(x,y)ifnarginbeta0=[11]";beta=nlinfit(x,y,@mymodel,beta0);fprintf("Thenonlinearleastsquarefittingy=a*sin(b*x)fordata\\n\\n");fprintf("%6.1f",x);fprintf("%6.1f",y);fprintf("\\n\\nis\\n\\nty=%7.4f*sin(%7.4f*x)\\n\\n",beta);z=linspace(x(1),x(end),100);plot(x,y,"ro",z,beta(1)*sin(beta(2)*z),"b-.");functionyb=mymodel(beta,xb)yb=beta(1)*sin(beta(2)*xb);計(jì)算結(jié)果:
9.?dāng)M合形如f(x)≈(a+bx)/(1+cx)的函數(shù)的一種快速方法是將最小二乘法用于下列問題:f(x)(1+cx)≈(a+bx),試用這一方法擬合表4-4給出的中國人口數(shù)據(jù)。表4-4次序年份人口(億)a)b)c)d)e)19535.8219646.95198210.08199011.34201*12.66解:把f(x)(1+cx)≈(a+bx)變成f(x)≈a+bx-cxf(x)則近似看成基函數(shù)是1,x,-x*f(x)而數(shù)據(jù)是(xi,f(xi))的最小二乘擬合問題,程序如下:
function[a,b,c]=ex41x=[1953196419821990201*]";y=[5.826.9510.0811.3412.66]";A=[ones(5,1)x-x.*y];Z=A\\y;a=Z(1)b=Z(2)c=Z(3)z=linspace(1953,201*,100);plot(x,y,"ro",z,(a+b*z)./(1+c*z),"b-.");結(jié)果:
所以fx
2.94560.0014x。
10.0004956x10.已知20世紀(jì)美國人口的統(tǒng)計(jì)數(shù)字如表5-12(單位:百萬):
表5-12美國人口統(tǒng)計(jì)數(shù)據(jù)年份1900191019201930194019501960197019801990人口76.092.0106.5123.2131.7150.7179.3204.0226.5251.4試分別用兩點(diǎn)公式和三點(diǎn)公式計(jì)算美國人口20世紀(jì)的年增長率。
1采用兩點(diǎn)公式,設(shè)t(i)表示年份,y(i)表示該年份的美國人口;兩點(diǎn)公式為:解;○
l(i)兩點(diǎn)公式matlab編程;
t=[1900:10:1990];y(i1)y(i)
t(i1)t(i)y=[76.092.0106.5123.2131.7150.7179.3204.0226.5251.4];fori=1:10,ifi2采用三點(diǎn)公式,設(shè)t(i)表示年份,y(i)表示該年份的美國人口;三點(diǎn)公式為:○
l(i)三點(diǎn)公式matlab程序如下:t=[1900:10:1990];3y(i)4y(i1)y(i2)
t(i2)t(i)y=[76.092.0106.5123.2131.7150.7179.3204.0226.5251.4];fori=1:10,ifi
友情提示:本文中關(guān)于《MATLAB上機(jī)實(shí)驗(yàn)報(bào)告》給出的范例僅供您參考拓展思維使用,MATLAB上機(jī)實(shí)驗(yàn)報(bào)告:該篇文章建議您自主創(chuàng)作。
來源:網(wǎng)絡(luò)整理 免責(zé)聲明:本文僅限學(xué)習(xí)分享,如產(chǎn)生版權(quán)問題,請(qǐng)聯(lián)系我們及時(shí)刪除。