數(shù)值計(jì)算課程7090307龍凱翔

1、東北大學(xué)秦皇島分校數(shù)值計(jì)算課程設(shè)計(jì)報(bào)告系 別信息與計(jì)算科學(xué)專 業(yè)學(xué) 號(hào)7090307姓 名龍凱翔指導(dǎo)教師張建波 姜玉山成 績教師評(píng)語：指導(dǎo)教師簽字： 2011年7月13日1 緒 論數(shù)值分析是研究分析用計(jì)算機(jī)求解數(shù)學(xué)計(jì)算問題的數(shù)值計(jì)算方法及其理論的學(xué)科，是數(shù)學(xué)的一個(gè)分支，它以數(shù)字計(jì)算機(jī)求解數(shù)學(xué)問題的理論和方法為研究對(duì)象。MATLAB是主要面對(duì)科學(xué)計(jì)算、可視化以及交互式程序設(shè)計(jì)的高科技計(jì)算環(huán)境。它將數(shù)值分析、矩陣計(jì)算、科學(xué)數(shù)據(jù)可視化以及非線性動(dòng)態(tài)系統(tǒng)的建模和仿真等諸多強(qiáng)大功能集成在一個(gè)易于使用的視窗環(huán)境中，為科學(xué)研究、工程設(shè)計(jì)以及必須進(jìn)行有效數(shù)值計(jì)算的眾多科學(xué)領(lǐng)域提供了一種全面的解決方案。2 用M

2、ATLAB進(jìn)行插值計(jì)算題目下列數(shù)據(jù)點(diǎn)的插值X01491625364964Y012345678可以得到平方根函數(shù)的近似，在區(qū)間0, 64上作圖。(1) 用這9個(gè)點(diǎn)作8次多項(xiàng)式插值.(2) 用三次樣條（第一邊界條件）程序求.從得到的結(jié)果看在0, 64上，哪個(gè)插值更精確；在區(qū)間0, 1上，兩種插值哪個(gè)更精確？2.1 程序文件建立新的m文件unit2.m：代碼如下：clc; clear; x = 0, 1, 4, 9, 16, 25, 36, 49, 64; y = 0:8; xi = 0:64; p = polyfit(x, y, 8); yi = polyval(p, xi); subplot(2

3、, 2, 1); plot(x, y, 'x', xi, yi, 'k')xlabel('x軸'); ylabel('y軸'); legend('插值點(diǎn)', '8次差值多項(xiàng)式')yi = interp1(x, y, xi, 'cubic'); subplot(2, 2, 2); plot(x, y, 'x', xi, yi, 'k')xlabel('x軸'); ylabel('y軸'); legend('插值點(diǎn)&

4、#39;, '三次樣條插值')xi = 0:0.01:1; p = polyfit(x, y, 8); yi = polyval(p, xi); subplot(2, 2, 3); plot(xi, yi, 'k')xlabel('x軸'); ylabel('y軸'); legend('8次差值多項(xiàng)式')yi = interp1(x, y, xi, 'cubic'); subplot(2, 2, 4); plot(xi, yi, 'k')xlabel('x軸'); y

5、label('y軸'); legend('三次樣條插值') 圖12.2 圖樣運(yùn)行文件得到圖12.3 分析比較兩個(gè)函數(shù)的圖像可知，在區(qū)間0, 64上三次樣條插值函數(shù)要更精確一些，在區(qū)間0, 1上拉格朗日插值函數(shù)仍然不如三次樣條插值函數(shù)更精確。3 多項(xiàng)式曲線擬合題目由實(shí)驗(yàn)給出數(shù)據(jù)表X0.00.10.20.30.50.81.0Y1.00.410.500.610.912.022.46試求3次、4次多項(xiàng)式的曲線擬合，再根據(jù)曲線形狀，求一個(gè)另外函數(shù)的擬合曲線，用圖示數(shù)據(jù)曲線及相應(yīng)的三種擬合曲線。3.1程序文件建立新的m文件unit3.m：代碼如下：clc; clear; x

6、 = 0.0, 0.1, 0.2, 0.3, 0.5, 0.8, 1.0; y = 1.0, 0.41, 0.50, 0.61, 0.91, 2.02, 2.46; a = polyfit(x, y, 3); b = polyfit(x, y, 4); xi = 0.0:0.01:1.0; aa = polyval(a, xi); bb = polyval(b, xi); subplot(1, 2, 1); plot(x, y, 'x', xi, aa, 'k')xlabel('x軸'); ylabel('y軸'); legend

7、('插值點(diǎn)', '3次曲線擬合')subplot(1, 2, 2); plot(x, y, 'x', xi, bb, 'k')xlabel('x軸'); ylabel('y軸'); legend('插值點(diǎn)', '4次曲線擬合')poly2str(a, 'x')poly2str(b, 'x') 3.2 圖形和結(jié)果運(yùn)行程序結(jié)果如下：ans = - 6.6221 x3 + 12.8147 x2 - 4.6591 x + 0.92659ans

8、= 2.8853 x4 - 12.3348 x3 + 16.2747 x2 - 5.2987 x + 0.94272圖像如下：圖23.2 5次多項(xiàng)式曲線擬合下面進(jìn)行的：代碼如下：clc; clear; x = 0.0, 0.1, 0.2, 0.3, 0.5, 0.8, 1.0; y = 1.0, 0.41, 0.50, 0.61, 0.91, 2.02, 2.46; a = polyfit(x, y, 5); xi = 0.0:0.01:1.0; aa = polyval(a, xi); plot(x, y, 'x', xi, aa, 'k')xlabel(&#

9、39;x軸'); ylabel('y軸'); legend('插值點(diǎn)', '5次曲線擬合')poly2str(a, 'x')運(yùn)行程序結(jié)果如下：ans = - 79.3261 x5 + 195.4554 x4 - 172.7104 x3 + 69.0498 x2 - 11.0044 x + 0.99547圖像如下：圖3由圖像可知，次數(shù)越高，擬合越準(zhǔn)確。4 數(shù)值積分題目用不同數(shù)值方法計(jì)算積分.(1)取不同的步長h.分別用復(fù)合梯形及復(fù)合辛普森求積計(jì)算積分，給出誤差中關(guān)于h的函數(shù)，并與積分精確值比較兩個(gè)公式的精度，是否存在一個(gè)最小

10、的h，使得精度不能呢個(gè)再被改善？(2)用龍貝格求積計(jì)算完成問題(1).(3)用自適應(yīng)辛普森積分，使其精度達(dá)到4.1 復(fù)合梯形公式代碼建立新的m文件unit41.m：代碼如下：clc, clear; x = 0.0000001:0.1:1; y = x.(1 / 2). * log(x); Fx = trapz(x, y)error = Fx + 4 / 9運(yùn)行文件結(jié)果為：Fx = - 0.4123error = 0.0321取不同步長h:h = 0.01，F(xiàn)x = - 0.4431，error = - 0.0014；h = 0.001，F(xiàn)x = - 0.4444，error = - 5.487

11、5e - 005；h = 0.0001，F(xiàn)x = - 0.4444，error = - 2.0338e - 006；h = 0.00001，F(xiàn)x = - 0.4444，error = - 6.4496e - 008.可見，沒有一個(gè)最小的h使精度不能再改變。4.2 復(fù)合辛普森公式代碼建立新的m文件unit42.m：代碼如下：Fx = quad(inline('x.(1 / 2). * log(x)'), 0.000001, 1，0.0001)error = Fx + 4 / 9運(yùn)行文件結(jié)果為：Fx = - 0.4440error = 4.3273e - 004取不同的步長h：h

12、= 0.00001，F(xiàn)x = - 0.4444，error = - 6.3537e - 005；h = 0.000001，F(xiàn)x = - 0.4444，error = - 2.9938e - 006；h = 0.0000001，F(xiàn)x = - 0.4444，error = - 3.0657e - 007；h = 0.0000000000001，F(xiàn)x = - 0.4444，error = - 9.6548e - 009；h = 0.00000000000001，F(xiàn)x = - 0.4444，error = - 9.6548e - 009?？梢姡嬖谝粋€(gè)最小的h使精度不能在改善，且h在0.0000000

13、000001與0.00000000000001之間。4.3龍貝格算法建立新的m文件unit43.m：代碼如下：function q, n = unit43(f, a, b); M = 1; abs0 = 10; k = 0; T = zeros(1, 1); h = b - a; T(1, 1) = (h / 2) * (subs(f, a) + subs(f, b); while abs0>0.0001 k = k + 1; h = h / 2; p = 0; for i = 1:M x = a + h * (2 * i - 1); p = p + subs(f, x); end T(

14、k + 1, 1) = T(k, 1) / 2 + h * p; M = 2 * M; for j = 1:k T(k + 1, j + 1) = (4j) * T(k + 1, j) - T(k, j) / (4j - 1); end abs0 = abs(T(k + 1, j + 1) - T(k, j); endq = T(k + 1, k + 1); n = k; 在命令窗口輸入：Fx, n = unit43('sqrt(x) * log(x)', 10( - 8), 1)輸出結(jié)果為：Fx = - 0.4444n = 94.4 自適應(yīng)辛普森算法Fx = quad(inl

15、ine('x.(1 / 2). * log(x)'), 0.000001, 1)輸出結(jié)果為：Fx = - 0.4444。5 解線性方程題目線性方程組Ax = b的A及b為A = ，b = ，則解x = .用MATLAB內(nèi)部函數(shù)求detA及A的所有特征值和cond(A) .若令A(yù) + A = ，求解(A + A)(x + x) = b, 輸出向量x和.從理論結(jié)果和實(shí)際計(jì)算兩方面分析線性方程組Ax = b解的相對(duì)誤差及A的相對(duì)誤差的關(guān)系.5.1 求detA及A的所有特征值和cond(A)輸入程序：A = 10 7 8 7; 7 5 6 5; 8 6 10 9; 7 5 9 10de

16、t(A)V, D = eig(A)norm(A, 2)輸出結(jié)果為：A = 10 7 8 7 7 5 6 5 8 6 10 9 7 5 9 10ans = 1V = 0.5016 - 0.3017 0.6149 0.5286 - 0.8304 0.0933 0.3963 0.3803 0.2086 0.7603 - 0.2716 0.5520 - 0.1237 - 0.5676 - 0.6254 0.5209D = 0.0102 0 0 0 0 0.8431 0 0 0 0 3.8581 0 0 0 0 30.2887ans = 30.2887則得到detA = 1A的特征值為0.0101 0.

17、8431 3.8581 30.2887 = 30.28876 解線性方程組的迭代法題目給出線性方程組其中系數(shù)矩陣為希爾伯矩陣：，i = 1, 2n.假設(shè)，若取n = 6, 8, 10，分別用雅克比迭代及SOR迭代(w = 1, 1.25, 1.5)求解.比較計(jì)算結(jié)果.6.1 用雅可比迭代法求解建立新的m文件unit61.m：代碼如下：n = 6; H = hilb(n); b = H * ones(n, 1); e = 0.00001; for i = 1:n if H(i, i) = = 0 '對(duì)角元為零，不能求解' return endendx = zeros(n, 1);

18、 k = 0; kend = 10000; r = 1; while k< = kend&r>e x0 = x; for i = 1:n s = 0; for j = 1:i - 1s = s + H(i, j) * x0(j); end for j = i + 1:ns = s + H(i, j) * x0(j); end x(i) = b(i) / H(i, i) - s / H(i, i); end r = norm(x - x0, inf); k = k + 1; end if k>kend'迭代不收斂，失敗' else'求解成功

19、9;x k end輸出結(jié)果為：ans = 迭代不收斂，失敗6.2 用SOR迭代法求解建立新的m文件unit62.m：代碼如下：function s = unit62 (n, w); H = hilb(n); b = H * ones(n, 1); e = 0.00001; for i = 1:n if H(i, i) = = 0 '對(duì)角元為零，不能求解' return endendx = zeros(n, 1); k = 0; kend = 10000; r = 1; while k< = kend&r>e x0 = x; for i = 1:n s = 0

20、; for j = 1:i - 1s = s + H(i, j) * x(j); end for j = i + 1:ns = s + H(i, j) * x0(j); end x(i) = (1 - w) * x0(i) + w / H(i, i) * (b(i) - s); end r = norm(x - x0, inf); k = k + 1; end if k>kend'迭代不收斂，失敗' else '求解成功'xk end輸入：unit62(6, 1)輸出結(jié)果為：ans = 求解成功x = 0.99931.01310.95371.03741.0

21、2960.9662k = 997ans = 0.6487當(dāng)n = 6，w = 1.25時(shí)，ans = 求解成功，x = 0.9992, 1.0131, 0.9558, 1.0375, 1.0197, 0.9741，k = 3121，ans = 0.6480；當(dāng)n = 6，w = 1.5時(shí)，ans = 求解成功，x = 0.9995, 1.0081, 0.9726, 1.0232, 1.0123, 0.9839，k = 5635，ans = 0.6471；當(dāng)n = 8，w = 1時(shí)，ans = 求解成功，x = 0.9996, 1.0039, 0.9938, 0.9824, 1.0071，1.0

22、200，1.0093, 0.9790 ，k = 3808，ans = 0.6601；當(dāng)n = 8，w = 1.25時(shí)，ans = 求解成功x = 0.9998, 1.0004, 1.0103, 0.9685, 1.0109, 1.0201, 1.0100, 0.9795，k = 3949，ans = 0.6601；當(dāng)n = 8，w = 1.5時(shí)，ans = 求解成功，x = 1.0000, 0.9975, 1.0231, 0.9469, 1.0275, 1.0118, 1.0158, 0.9771，k = 4031，ans = 0.6602；當(dāng)n = 10，w = 1時(shí)，ans = 求解成功，

23、x = 1.0001, 0.9952, 1.0275, 0.9677, 0.9816, 1.0093, 1.0241, 1.0208, 1.0018, 0.9713，k = 2673，ans = 0.6677, ；當(dāng)n = 10，w = 1.25時(shí)，ans = 求解成功，x = 1.0004, 0.9904, 1.0466, 0.9421, 0.9886, 1.0123, 1.0272, 1.0214, 1.0010, 0.9694，k = 2562，ans = 0.6677；當(dāng)n = 10，w = 1.5時(shí)，ans = 求解成功，x = 1.0006, 0.9857, 1.0688, 0.9

24、028, 1.0177, 0.9993, 1.0377, 1.0169, 1.0030, 0.9669 ，k = 2463，ans = 0.66797 用MATLAB求解非線性方程的根題目求下列方程的實(shí)根(1) ; (2) .要求：(1) 設(shè)計(jì)一種不動(dòng)點(diǎn)迭代法，要使迭代序列收斂，然后再用斯特芬森加速迭代，計(jì)算到為止.(2) 用牛頓迭代，同樣計(jì)算到.輸出迭代初值及各次迭代值和迭代次數(shù)k，比較方法的優(yōu)劣.7.1 不動(dòng)點(diǎn)迭代法建立新的m文件unit71.m：代碼如下：function y, n = unit71(x, eps)if nargin = = 1 eps = 1.0e - 8; else

25、if nargin<1 error returnendx1 = gg(x); n = 1; while (norm(x1 - x)> = 1e - 8)&&(n< = 10000) x = x1; x1 = gg(x); n = n + 1; endy = x; 若是function f = gg(x)f(1) = (exp(x) - 2 - x2) / ( - 3); 在在命令窗口輸入：unit71(1)輸出結(jié)果為：n = 15ans = 0.2575若是function f = gg(x)f(1) = - (x3 + 2 * x2 - 20) / 10;

26、在命令窗口輸入：unit71(1)輸出結(jié)果為：n = 10001ans = 0.54897.2 斯特芬森加速迭代法建立新的m文件unit72.m：代碼如下：function y = unit72(x0, eps)if nargin<2 eps = 1.0e - 8; enda = g(x0); b = g(a); x1 = x0 - (a - x0)2 / (b - 2 * a + x0); n = 1; while (abs(x1 - x0)> = eps)&(n< = 10000) x0 = x1; a = g(x0); b = g(a); x1 = x0 - (a - x0)2 / (b - 2 * a + x0); n = n + 1; endx1n若是function y = g(x); y = (exp(x) - 2 - x2) / ( - 3); 在命令窗口輸入：unit72(1)輸出結(jié)果為：x1 = 0.2575n = 4若是function y = g(x); y = - (x3 + 2 * x2 - 20) / 10在命令窗口輸入：unit72(1)輸出結(jié)果為：x1 = 1.3688n = 57.3 牛頓迭代法建立新的m文件unit73.m：代碼如下：function y = unit73(x0, eps)if narg