1、1 .應用自己熟悉的算法語言編寫程序，使之盡可能具有通用性2 .上機前充分準備，復習有關算法，寫出計算步驟，反復檢查，調試程序。（注：在練習本上寫，不上交）3 .完成計算后寫出實驗報告， 內容包括：算法步驟敘述， 變量說明，程序清單， 輸 出計算結果，結構分析和小結等。（注：具體題目具體分析，并不是所 有的題目的實驗報告都包含上述內容！）4 .獨立完成，如有雷同，一律判為零分！5 .上機期間不允許做其他任何與課程設計無關的事情，否則被發現一次扣10分，被發現三次判為不及格！非特殊情況，不能請假。曠課3個半天及以 上者，直接判為不及格。一、基本技能訓練 31、誤差分析 32、求解非線性方程 33

2、、插值 54、數值積分 9二、提高技能訓練 111、 112、 12三、本課程設計的心得體會（500字左右） 14一基本技能訓練1、誤差分析實驗1? 2誤差傳播與算法穩定性實驗目的：體會穩定性在選擇算法中的地位。誤差擴張的算法是不穩定的，是我們所不期望的；誤差衰減的算法是穩定的，是我們努力尋求的，這是貫穿本課程的目標。問題提出：考慮一個簡單的由積分定義的序列iiIn 二 xnex dx,n =1,2,3 JI I顯然In0,n = 1,2,3,|。當n =1時，j而對于n 2時，利用分部積0 xe d分易得dxIn = oxnex'dx = x nex-'。nxn'ex

3、'dx = 1 - n 1 心，n = 2,3,|l 另一方面，我們有In實驗內容：由以上遞推關系，我們可得到計算序列m的兩種方法。1)I1, In = 1 一 nIn n 二 2,3,HIe1 - E(h)En =0,巳_1N,n 二 N,N -1,N -2,|,3,2nsyms n In5 In6 In7 ;In 5=vpa(exp(-1),5);In 6=vpa(exp(-1),6);,eval(In5),eval(In6),eval(In7);In 7=vpa(exp(-1),7); fprintf( %5f %.6f %.7fn' for n=2:10In 5=vpa

4、(1- n*ln 5),5);In 6=vpa(1- n*ln 6),6);fprintf('%.5f %.6f %.7fn'In 7=vpa(1- n*ln 7),7);,eval(I n5),eval(l n6),eval(l n7);end0.367880.264240.207280.170890.145530.126800.112380.100930.091610.083880.367 8790.2642410.2072770.1708930.1455330.126802五位六位七位0.36787940.26424110.20727660.17089340.145532

5、90.12680240.11238350.10093200.09161230.0838771syms nEn;En 5=vpa(0,5);En 6=vpa(0,6);En 7=vpa(0,7);fprintf( %5f %.6f %.7fn' ,eval(En5),eval(En6),eval(En7);for n=10:-1:2En 5=vpa(1-E n5)/n),5);En 6=vpa(1-E n6)/n),6);En 7=vpa(1-E n7)/n),7);fprintf( %.5f %.6f %.7fn',eval(E n5),eval(E n6),eval(E n7

6、);end五位六位七位0.00000 0.000000 0.00000000.10000 0.100000 0.10000000.10000 0.100000 0.10000000.11250 0.112500 0.11250000.12679 0.126786 0.12678570.14554 0.145536 0.14553570.17089 0.170893 0.17089290.20728 0.207277 0.20727680.26424 0.264241 0.26424110.36788 0.367879 0.3678795(1)由所得數據可以了解，兩種算法隨著n的增大，誤差越來越

7、大，而第二種算法隨著的減小，數據越來越精確，且三種有效數字結果大致相同，所以第二種更精確。(2) 算法一 E1 的誤差 e1,由于 En=1-n*En-1 ,n=2,3,4 .,通過推算可得誤差 eA_(-1 ) nn!e1算法EN的誤差；N,由于"”席導誤差,3,2(N -n)!;n比；N縮小了( N-n)!倍，則n減少，所以-就越小。所以算法越往后算，一步步的誤差會使誤差越大，而算法二由于是從后面遞推回來，其誤差會被縮小所以算法二更優。(3) 算法一中，n增大，誤差增大，算法二中，n減小，誤差減小。所以當某一步發生誤差后隨著n變大，算法一的誤差越來越大，而算法二由于是往回推，所以

8、誤差變小。(4) 通過以上可一推由，算法一隨著n變大，容易使誤差越來越大，而算法二會使誤差變小，所以算法二更加穩定。2、求解非線性方程3、插值實驗目的：掌握 Lagrange插值法和Newton插值法問題提由：設f(x)=eAdt,已知f(x)的函數值表如下丁 2兀-00x00.10.20.30.4f(x)0.50000.53980.57930.61790.7554用插值法求f(0.13)和f(0.36)的近似值實驗內容：(1)分別用Lagrange插值法和Newton插值法編程求解(2)求生插值多項式系數，對比計算結果。拉格朗日差值fun cti on yh=lagra nge(x,y,xh

9、) tic n = len gth(x);m = len gth(xh);x = x(:);y = y(：);xh = xh(:);yh = zeros(m,1); cl = ones(1,n-1);c2 = ones(m,1); for i=1:n,xp = x(1:i-1 i+1:n);yh = yh + y(i) * prod(xh*c1-c2*xp')./(c2*(x(i)*c1-xp'),2);endtocx=0 0.1 0.2 0.3 0.4;y=0.5000 0.5398 0.5793 0.6179 0.7554; xh=0.13 0.36;lagrange(x,

10、y,xh)時間已過 0.026678秒。ans =0.55370.6780多項式為(5*xh - 1)*(5*xh)/2 - 1)*(10*xh - 1)*(10*xh)/3 - 1)/2 + (5793*xh*(5*xh - 2)*(10*xh- 1)*(10*xh - 3)/2000 - (6179*xh*(5*xh - 1/2)*(10*xh - 2)*(10*xh - 4)/3000 +(3777*xh*(5*xh - 1)*(10*xh - 3)*(10*xh)/3 - 1/3)/2000 - (2699*xh*(5*xh - 3/2)*(10*xh - 2)*(10*xh)/3-

11、4/3)/500作圖x0=0 0.1 0.2 0.3 0.4;y0=0.5000 0.5398 0.5793 0.6179 0.7554; xh=0:0.01:0.4;y=lagrange(x0,y0,xh);時間已過0.000158秒。- > plot(xh,y);hold on;牛頓差值function f=Newton(x,y,x0,x1)ticsyms t ;if (length(x)=length(y) n=length(x); c(1:n)=0.0;elsedisp( 'x 和 y 的維數不相等！' ); return ; end f=y(1); y1=0;l

12、 =1;for (i=1:n-1) for (j=i+1:n)y1(j)=(y(j)-y(i)/(x(j)-x(i)； endc(i)=y1(i+1); l=l*(t-x(i); f=f+c(i)*l;y=y1; end f=simplify(f); g=subs(f, 't' ,x0) g1=subs(f, 't' ,x1)A=zeros(n,n-1); A=y',A; for j=2:n for i=j:nA(i,j)=(A(i,j-1)-A(i-1,j-1)/(x(i)-x(i+1-j); end end disp(' 差商表為')

13、; disp(A); tocx0=0 0.1 0.2 0.3 0.4;y0=0.5000 0.5398 0.5793 0.6179 0.7554; Newton(x0,y0,0.13,0.36)g = 249341921721956909228481/450359962737049600000000 g1 =1192694200291039391741/1759218604441600000000差商表為1.0e+0400.00000.00000.00000.00420 00.0004-0.0004-0.00010.041900-0.00410.00160.210100 0000.019000

14、.69481.6896時間已過 0.368255 秒ans =(251* 護 4)/6-(2269814212194733997M3)/90071992547409920(12474970967816329501*八2)/2702159776422297600(659777345409770407*t)/4503599627370496000 + 1/2作圖>> t=0:0.01:0.4;> > f=;> >f=(251.*t. A4)./6.-(2269814212194733997.*t. A3)./90071992547409920.+(1247497

15、0967816329501 氏八 2)./2702159776422297600.+(659777345409770407 氏)./4503 599627370496000.+ 1/2;> > plot(t,f);hold on;兩種方法作圖曲線為454036302620 -10 -ui i i i i i i0005010.15020 25030 350.4拉格朗日插值法時間復雜度為O(n),牛頓插值法時間復雜度為0(nT),實際上拉格朗日算法計算時間為0.000158,牛頓算法計算時間 0.368255 。說明拉格朗日算法時間復雜度低。4、數值積分實驗目的：掌握Romberg積

16、分法實驗內容：1 dx用Romberg積分法計算下列定積分二 4j cos(s in (x)dx32 xx|n(3x)dxfun cti on R,y=Romberg(a,b ,n) k=1;while 1T=zeros(k+1,k+1);T(1,1)=1/2*(b-a)*(f(a)+f(b); for i=1:k h=(b-a)/2Ai;s=0;for j=1:2A(i-1)s=s+f(a+(2*j-1)*h);endT(i+1,1)=T(i,1)/2+h*s;endfor j=1:kc=1/(4Aj-1);for m=j:kT(m+1,j+1)=T(m+1,j)+c*(T(m+1,j)-T

17、(m,j); endendif abs(T(k,k)-T(k+1,k+1)<nbreak ;endk=k+1;endy=T(k,k);R=T;輸入積分公式(1)function y=f(x)y=1/(1+x a4);R,y=Romberg(0,1,10)R =0.7500 00.8456 0.8775 y =0.7500(2) function y=f(x) y=cos(sin(x);R,y=Romberg(0,pi/4,10)R =0.6912 00.7099 0.71610.6912(3)function y=f(x)y=exp(x)/x;R,y=Romberg(1,2,10)R =

18、3.2064 03.0971 3.06073.2064(4)function y=f(x) y=x4x*log(3+x);R,y=Romberg(2,3,10)R =27.4076 022.1271 20.3669y =27.4076龍貝格計算量較小二、提高技能訓練1、(2)實驗目的與要求(1) 了 解 Kepler 方程；(2)能夠用牛頓迭代方法計算Kepler方程；(3) 通過調節軌道偏心率 e查看運算收斂情況o實驗內容及數據來源編寫解kepler方程的方法，給定偏心率為 0.01、平近點角為32度時計算由偏近點角 Eclear;x0=32;syms E;f=i nline('E-

19、0.01*si n( E)-32');g=i nli ne( '1-0.01*cos(E)');x1=x0-f(x0)/g(x0);y=abs(x1-x0);x0=x1;i=i+1;while y>10 A(-5)x1=x0-f(x0)/g(x0);y=abs(x1-x0);x0=x1;i=i+1;endE=vpa(x1)結果：E =32.005560568887858607922680675983（3）通過調節軌道偏心率e查看運算收斂情況；e 0.001 1E 32.00055導 0.99916可以看由， 后，發散.2、1.832.9999133.535591.

20、013181.93911e在01之間，收斂；在2.21030.8838131.35097-0.8958211.8之間，-8.97892發散；在1.82.2左右，收斂;2.2之（日照時間分布）已知某地區在不同月份的平均日照時間的觀測數據如下表所示月份123456789101112日照/（h/月）80.967廠267.150. 532.033.635廠646.8；52 362.064.171.2試分析日照時間的變化規律Ox=1:1:12;y=80.9 67.2 67.1 50.5 32.0 33.6 36.6 46.8 52.3 62.0 64.1 71.2; x1=1:0.01:12;y1=i

21、nterp1(x,y,x1);y2=in terp1(x,y,x1, 'spli ne');plot(x,y, 'go' ,x1,y1, 'r-',x1,y2, 'b')其中紅線是分段插值得到的，因此為折線，綠色是三次樣條曲線，0 '處是插值節點。從圖中可以看生，日照時間從一月開始減少至5月，之后又開始增加至 12月。（水道測量問題）水深數據如下表所示，設船的吃水深度為 1.5m,問在70,200-100,150 m 2的區域內，哪些地方船只避免進入，畫船只避免進入的區域。x129.0140.0108.588.0185.5

22、195.5105.5y7.5141.528.0147.022.5137.585.5z1.222.441.832.441.832.442.44x157.5107.577.081.0162.0162.0117.5y-6.5-81.53.056.5-66.584-38.5z2.742.742.442.442.741.222.74x=129.0 140.0 108.5 88.0 185.5 195.5 105.5 157.5 107.5 77.0 81.0 162.0 162.0 117.5; y=7.5141.5 28.0 147.0 22.5 137.5 85.5 -6.5 -81.5 3.0 56.5 -66.5 84 -38.5;z=-1.22 2.44 1.83 2.44 1.83 2.44 2.44 2.74 2.74 2.44 2.44 2.74 1.22 2.74;x1=li nspace( min( x),max(x),40);y1=li nspace(mi n(y) ,max(y)