Matlab上機實驗報告

實驗目的

學會用數(shù)學軟件matlab和蒙特卡洛方法估計積分值，并于其中應用概率論中的概率密度等知識點。實驗問題

估計,,的值，將估計值與真實值進行比較；估計,的值，并對誤差進行估計。實驗要求

針對要估計的積分選擇適當?shù)母怕史植荚O計蒙特卡洛估計算法；利用計算機產(chǎn)生所選分布的隨機樣本的估計積分值；通過計算平均誤差對估計結果進行評價。

實驗過程（x為運行結果平均值，y為樣本方差）估計程序如下：clc;clear;m=10;n=10000;d=0;e=0;fori=1:md=0;

a=rand(1,n);forj=1:nb=a(j)+2;c=b^2;e(i)=d+c/n;d=e(i);end

fprintf("e=%.8f\\n",e(i))end

p=sum(e)/m;forj=1:m;

s(j)=(e(j)-p)^2;end

q=sum(s);

fprintf("x=%.8f\\ny=%.8f\\n",p,q);結果為：

e=6.34879520e=6.34068140e=6.35081124e=6.31353632e=6.35586630e=6.33058791e=6.32419121e=6.33707454e=6.30357011

e=6.35063255x=6.33557468y=0.00304243

2.估計程序變動如下：b=a(j)*3;c=b*sin(b)*3;結果為：

e=3.12211717e=3.11373037e=3.07484948e=3.08660758e=3.10052243e=3.10475698e=3.13762746e=3.16481618e=3.1155201*e=3.09615989x=3.11167076y=0.00669223

3.估計程序變動如下：b=exp(-a(j)^2/2);c=b/n*(2*pi)^0.5;e(i)=d+c/2;結果為：

e=0.88617655e=0.88538972e=0.88635209e=0.88575809e=0.88653705e=0.88606366e=0.88634011e=0.88613926e=0.88573325e=0.88644154x=0.88609313y=0.00000132

4.估計程序變動如下：b=exp(a(j)^2);c=b/n;e(i)=d+c;

2

結果為：

e=1.46211146e=1.46154792e=1.46327379e=1.46256348e=1.46318297e=1.46235828e=1.46241378e=1.46316145e=1.46203052e=1.46280489x=1.46254485y=0.00000320

5.估計程序變動如下：b=a(i)*4;

c=1/((1+b^2)^0.5);e(i)=d+c*4/n;結果為：

e=1.98511173e=1.02167881e=1.26713031e=0.98770837e=1.14216662e=1.75642022e=1.97055988e=1.96227794e=1.83229787e=1.06190231x=1.49872541y=1.90228258

實驗結果

通過對實驗所得平均值與真實值的比較，可以看出實驗結果與真實值相比非常接近，而且樣本方差很小，從而說明概率分布的選取比較適當，計算機實驗很準確。

擴展閱讀：matlab上機實驗報告

上海電力學院

數(shù)值計算方法上機實習

院系：專業(yè)年級：學生姓名：學號：指導教師：報告

能源與機械工程學院動力機械及工程201*級張亞杰ys1210124014黃建雄

201*年12月26日1

數(shù)值計算方法上機實習題

xndx，1．設In05x1（1）由遞推公式In5In1解：I0=

I=0.1823forn=1:1:20,I=-5*I+1/n;1，從I0的幾個近似值出發(fā)，計算I20；n105xdx=0.1823

1計算I20編輯matlab命令如下：fprintf("%.1d%.4f\\n",n,I);end結果：

（2）粗糙估計I20，用In111In，計算I0；55nx20d解：I20=05xx1使用復合中點公式進行積分，相應的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)結果：

程序結束時輸出兩個I值，第一個表示I20，第二個表示I0；分別為I20=0.0082I0=0.1823

（3）分析結果的可靠性及產(chǎn)生此現(xiàn)象的原因(重點分析原因)。

從上述計算中分析得到如果先得到I0，再從I0由遞推公式得到I20，I20結果跟精確值相比誤差很大；如果先估算I20，在從I20有遞推公式得到I0，I0的結果跟精確值相比近似相等。

原因分析：

如果從I0推I20的近似值，需要用到遞推公式In=-5In-1+1/n，I0本身結果是有誤差的；經(jīng)過遞推公式計算20次，就等于誤差被認為的放大5的20次方倍，所以得到的I20與其精確值相差甚遠。

如果從I20推I0的近似值，需要用到In-1=0.2（1/n-In），盡管I20本身有誤差，但是經(jīng)過20次運算，其誤差縮小到原來的0.2的20次方倍，所以得到的I0與其精確值比較相近。

2．求方程e10x20的近似根，要求xk1xk5104，并比較計算量。

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進行迭代的程序如下：

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結果：

解得x=0.0905

（3）加速迭代的結果；

采用matlab進行迭代的程序如下：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結果如下：

（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運行結果：

（5）分析絕對誤差。迭代次數(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我們可以看到，在運算要求到同一精度的情況下，采用（1）的二分法運算了11次，采用（2）的方法運算了4次，采用（3）的加速迭代法運算了2次，采用（4）的牛頓迭代法也需運算2次。也就是說牛頓的迭代的收斂速度與加速迭代速度都是超線性收斂的，而簡單迭代法是線性收斂的。而二分法收斂速度較慢,所以在工程上不經(jīng)常使用。3．鋼水包使用次數(shù)多以后，鋼包的容積增大，數(shù)據(jù)如下：xy10234569.7147101581112131610.4910.5910.6010.810.610.910.76試從中找出使用次數(shù)和容積之間的關系，計算均方差。（注：增速減少，用何種模型）解：將使用次數(shù)x與體積y的關系用matlab采用如下程序繪制在二維坐標系: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結果如下：

由數(shù)據(jù)點分布圖可知，擬合曲線y=f(x)隨著x的增加而上升，但上升速度由快到慢，當x趨于無窮大時，y趨于某個常數(shù)，故曲線有一水平漸進線。根據(jù)上述特征很容易想到用

b/x

Logistic模型來擬合該曲線。設y=f(x)的形式為y=ae（a>0,b

計算均方差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);運算結果均方差S=0.2438

小結：根據(jù)已給的條件計算函數(shù)是十分困難的，但通過對離散點的分析及變化規(guī)律找出其中的規(guī)律，并通過計算來得到實際的函數(shù)是十分有用的方法。本題就是這樣做的一個典型，在n=1/x和m=lny的基礎上找到了它們之間的關系并通過這種關系來擬合原函數(shù)，并最終驗證計算結果。

041010014101050141012，b，Axb4．設A101410501014120601014分析下列迭代法的收斂性，并求xk1xk（1）JACOBI迭代；解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,

8

2104的近似解及相應的迭代次數(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的特征值和特征向量，準確到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);運行結果如下：

即離11最近的特征值為7.8745；相應的特征向量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公式計算的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結果如下所示：

2y1y22sinxy1（2），x[0,10]，

y998y999y999cosx999sinx122y1(x)2exsinx和精確解比較，分析結論。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結果如下：

精確解：

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結果;

結果分析：

四階RungeKutta方法得到的結果已很接近精確解，證明這種迭代方法精確度很好，是一種有效的算法。但是要注意龍格-庫塔公式的推導基于泰勒展開方法，因而它要求所求的的解具有較好的光滑性質(zhì)。反之，如果解得光滑性差，那么，使用四階龍格-庫塔求得的數(shù)值解精度就不是太高，此種情況可以采用縮小步長來解決，比如上述計算。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

運算結果為：

g10000000000000000000100

H矩陣為：

Y在各點的近似值為：

XY8.用函數(shù)y=asin(bx)擬合數(shù)據(jù).

x0.10.20.30.40.50.60.70.8y0.61.11.61.82.01.91.71.3Matlab上機程序為：function[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);計算結果：

9．擬合形如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-.");結果：

所以fx

2.94560.0014x。

10.0004956x10.已知20世紀美國人口的統(tǒng)計數(shù)字如表5-12(單位：百萬)：

表5-12美國人口統(tǒng)計數(shù)據(jù)年份1900191019201930194019501960197019801990人口76.092.0106.5123.2131.7150.7179.3204.0226.5251.4試分別用兩點公式和三點公式計算美國人口20世紀的年增長率。

1采用兩點公式，設t(i)表示年份，y(i)表示該年份的美國人口；兩點公式為：解;○

l(i)兩點公式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采用三點公式，設t(i)表示年份，y(i)表示該年份的美國人口；三點公式為：○

l(i)三點公式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

友情提示：本文中關于《Matlab上機實驗報告》給出的范例僅供您參考拓展思維使用，Matlab上機實驗報告：該篇文章建議您自主創(chuàng)作。

來源：網(wǎng)絡整理 免責聲明：本文僅限學習分享，如產(chǎn)生版權問題，請聯(lián)系我們及時刪除。

《Matlab上機實驗報告》由互聯(lián)網(wǎng)用戶整理提供,轉(zhuǎn)載分享請保留原作者信息,謝謝!
鏈接地址：http://www.weilaioem.com/gongwen/712933.html

• 上一篇：袁越強MATLAB上機實驗報告一
• 下一篇：第二次DSP上機實驗報告 MATLAB實驗報告