計算方法綜合實踐課件

1、1 應(yīng)用自己熟悉的算法語言編寫程序，使之盡可能具有通用性。2 上機前充分準備，復(fù)習(xí)有關(guān)算法，寫出計算步驟，反復(fù)檢查，調(diào)試程序。（注：在練習(xí)本上寫，不上交）3 完成計算后寫出實驗報告，內(nèi)容包括:算法步驟敘述，變量說明，程序清單，輸出計算結(jié)果，結(jié)構(gòu)分析和小結(jié)等。（注：具體題目具體分析，并不是所有的題目的實驗報告都包含上述內(nèi)容?。? 獨立完成，如有雷同，一律判為零分！5 上機期間不允許做其他任何與課程設(shè)計無關(guān)的事情，否則被發(fā)現(xiàn)一次扣10分，被發(fā)現(xiàn)三次判為不及格！非特殊情況，不能請假。曠課3個半天及以上者，直接判為不及格。目 錄一、基本技能訓(xùn)練41、誤差分析42、求解非線性方程43、插值44、數(shù)值積分

2、4二、提高技能訓(xùn)練41、42、4三、本課程設(shè)計的心得體會（500字左右）4一、基本技能訓(xùn)練1、誤差分析實驗 1.2誤差傳播與算法穩(wěn)定性實驗?zāi)康模后w會穩(wěn)定性在選擇算法中的地位。誤差擴張的算法是不穩(wěn)定的，是我們所不期望的；誤差衰減的算法是穩(wěn)定的，是我們努力尋求的，這是貫穿本課程的目標。問題提出：考慮一個簡單的由積分定義的序列顯然。當(dāng)時，而對于時，利用分部積分易得另一方面，我們有實驗內(nèi)容：由以上遞推關(guān)系，我們可得到計算序列的兩種方法。（I）（II）syms n In5 In6 In7;In5=vpa(exp(-1),5);In6=vpa(exp(-1),6);In7=vpa(exp(-1),7);f

3、printf(%.5f %.6f %.7fn,eval(In5),eval(In6),eval(In7);for n=2:10 In5=vpa(1-n*In5),5); In6=vpa(1-n*In6),6); In7=vpa(1-n*In7),7); fprintf(%.5f %.6f %.7fn,eval(In5),eval(In6),eval(In7);end五位 六位 七位0.36788 0.367879 0.36787940.26424 0.264241 0.26424110.20728 0.207277 0.20727660.17089 0.170893 0.17089340.14

4、553 0.145533 0.14553290.12680 0.126802 0.12680240.11238 0.112384 0.11238350.10093 0.100932 0.10093200.09161 0.091612 0.09161230.08388 0.083877 0.0838771syms nEn;En5=vpa(0,5);En6=vpa(0,6);En7=vpa(0,7);fprintf(%.5f %.6f %.7fn,eval(En5),eval(En6),eval(En7);for n=10:-1:2 En5=vpa(1-En5)/n),5); En6=vpa(1-

5、En6)/n),6); En7=vpa(1-En7)/n),7); fprintf(%.5f %.6f %.7fn,eval(En5),eval(En6),eval(En7);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.

6、20727680.26424 0.264241 0.26424110.36788 0.367879 0.3678795（1） 由所得數(shù)據(jù)可以了解，兩種算法隨著n的增大，誤差越來越大，而第二種算法隨著n的減小，數(shù)據(jù)越來越精確，且三種有效數(shù)字結(jié)果大致相同，所以第二種更精確。（2） 算法一E1的誤差e1，由于En=1-n*En-1，n=2,3,4.，通過推算可得誤差=，n越大，其誤差就越大，那么最后算出的結(jié)果是e1的n!倍；算法二的誤差，由于，推倒可得誤差=，比縮小了（N-n)!倍，則n減少，所以就越小。所以算法一越往后算，一步步的誤差會使誤差越大，而算法二由于是從后面遞推回來，其誤差會被縮小。所以

7、算法二更優(yōu)。（3） 算法一中，n增大，誤差增大，算法二中，n減小，誤差減小。所以當(dāng)某一步發(fā)生誤差 后隨著n變大，算法一的誤差越來越大，而算法二由于是往回推，所以誤差變小。（4） 通過以上可一推出，算法一隨著n變大，容易使誤差越來越大，而算法二會使誤差變小，所以算法二更加穩(wěn)定。2、求解非線性方程3、插值實驗?zāi)康模赫莆誏agrange插值法和Newton插值法問題提出：，已知的函數(shù)值表如下x00.10.20.30.40.50000.53980.57930.61790.7554用插值法求和的近似值。實驗內(nèi)容：（1）分別用Lagrange插值法和Newton插值法編程求解；（2）求出插值多項式系數(shù)，對

8、比計算結(jié)果。拉格朗日差值function yh=lagrange(x,y,xh)ticn = length(x);m = length(xh);x = x(:);y = y(:);xh = xh(:);yh = zeros(m,1); c1 = 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);endtoc x=0 0.1 0.2 0.3 0.4;y=0.5000 0.5398 0.5793 0.6179 0.755

9、4;xh=0.13 0.36;lagrange(x,y,xh)時間已過 0.026678 秒。ans = 0.5537 0.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

10、 - 3/2)*(10*xh - 2)*(10*xh)/3 - 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;else disp(x和y的維數(shù)不相等！);return;endf=y(1)

11、;y1=0;l =1;for(i=1:n-1)for(j=i+1:n) y1(j)=(y(j)-y(i)/(x(j)-x(i);end c(i)=y1(i+1); l=l*(t-x(i); f=f+c(i)*l; y=y1;endf=simplify(f);g=subs(f,t,x0)g1=subs(f,t,x1)A=zeros(n,n-1);A=y,A;for j=2:nfor i=j:n A(i,j)=(A(i,j-1)-A(i-1,j-1)/(x(i)-x(i+1-j);endenddisp(差商表為);disp(A);toc x0=0 0.1 0.2 0.3 0.4;y0=0.5000

12、 0.5398 0.5793 0.6179 0.7554;Newton(x0,y0,0.13,0.36)g =249341921721956909228481/450359962737049600000000g1 =1192694200291039391741/1759218604441600000000差商表為 1.0e+04 * 0 0 0 0 0 0.0000 0.0004 0 0 0 -0.0000 -0.0004 -0.0041 0 0 -0.0000 -0.0001 0.0016 0.0190 0 0.0042 0.0419 0.2101 0.6948 1.6896時間已過 0.3

13、68255 秒。ans =(251*t4)/6 - (2269814212194733997*t3)/90071992547409920 + (12474970967816329501*t2)/2702159776422297600 + (659777345409770407*t)/4503599627370496000 + 1/2作圖 t=0:0.01:0.4; f=; f=(251.*t.4)./6.-(2269814212194733997.*t.3)./90071992547409920.+ (12474970967816329501.*t.2)./270215977642229760

14、0.+(659777345409770407.*t)./4503599627370496000.+ 1/2; plot(t,f);hold on;兩種方法作圖曲線為拉格朗日插值法時間復(fù)雜度為O(n),牛頓插值法時間復(fù)雜度為O(n2),實際上拉格朗日算法計算時間為 0.000158 ， 牛頓算法計算時間 0.368255 。說明拉格朗日算法時間復(fù)雜度低。4、數(shù)值積分實驗?zāi)康模赫莆誖omberg積分法實驗內(nèi)容：用Romberg積分法計算下列定積分：function R,y=Romberg(a,b,n)k=1;while 1 T=zeros(k+1,k+1); T(1,1)=1/2*(b-a)*(f

15、(a)+f(b); for i=1:k h=(b-a)/2i; s=0; for j=1:2(i-1) s=s+f(a+(2*j-1)*h); end T(i+1,1)=T(i,1)/2+h*s; end for j=1:k c=1/(4j-1); for m=j:k T(m+1,j+1)=T(m+1,j)+c*(T(m+1,j)-T(m,j); end end if abs(T(k,k)-T(k+1,k+1)10(-5) x1=x0-f(x0)/g(x0); y=abs(x1-x0); x0=x1; i=i+1;endE=vpa(x1)結(jié)果：E =32.0055605688878586079

16、22680675983（3） 通過調(diào)節(jié)軌道偏心率 e 查看運算收斂情況；e0.001 1 1.8 2.2 10E32.0005532.9999133.5355930.8838131.35097導(dǎo)0.999161.013181.93911-0.89582-8.97892 可以看出，e在01之間，收斂；在11.8之間，發(fā)散；在1.82.2左右，收斂；2.2之后，發(fā)散.2、（日照時間分布）已知某地區(qū)在不同月份的平均日照時間的觀測數(shù)據(jù)如下表所示月份123456789101112日照/（h/月）80.967.267.150.532.033.635.646.852.362.064.171.2試分析日照時間

17、的變化規(guī)律。x=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=interp1(x,y,x1); y2=interp1(x,y,x1,spline);plot(x,y,go,x1,y1,r-,x1,y2,b)其中紅線是分段插值得到的，因此為折線，綠色是三次樣條曲線，o 處是插值節(jié)點。從圖中可以看出，日照時間從一月開始減少至5月，之后又開始增加至12月。（水道測量問題）水深數(shù)據(jù)如下表所示，設(shè)船的吃水深度為1.5m,問在的區(qū)域內(nèi)，哪些地方船只避免進入，畫船只避免進入的區(qū)域。

18、x129.0140.0108.588.0185.5195.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.5 141.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=linspace(min(x),max(x),40); y1=linspace(min(y),max(y),40);Xi,Yi=meshgrid(x1,y1);Xi,Yi,Zi=griddata(x,y,z,Xi,Yi,cu