【文章內(nèi)容簡(jiǎn)介】 86018427387904(5761533065688605*(x 900)*(x 1910)*(x1920)*(x1930))/73786976294838206464+(5016235413793213*(x1900)*(x1910)*(x1920)*(x1930)*(x1940))/1180591620717411303424(5436283353602131*(x1900)*(x1910)*(x1920)*(x1930)*(x1940)*(x1950))/37778931862957161709568 + (3912362197104435*(x1900)*(x1910)*(x1920)*(x1930)*(x 1940)*(x 1950)*(x 1960))/1208925819614629174706176 (3789694605011005*(x 1900)*(x 1910)*(x 1920)*(x 1930)*(x 1940)*(x 1950)*(x 1960)*(x 1970))/77371252455336267181195264 + (8429618423594317*(x 1900)*(x 1910)*(x 1920)*(x 1930)*(x 1940)*(x 1950)*(x 1960)*(x 1970)*(x 1980))/19807040628566084398385987584 + (7234635122020997*(x 1900)*(x 1910)*(x 1920)*(x 1930)*(x 1940)*(x 1950)*(x 1960)*(x 1970)*(x 1980)*(x 1990))/5070602400912917605986812821504 (6099905420751575*(x 1900)*(x 1910)*(x 1920)*(x 1930)*(x 1940)*(x 1950)*(x 1960)*(x 1970)*(x 1980)*(x 1990)*(x 2020))/40564819207303340847894502572032 591927/200 N0 = Newton 插值 2020 年的產(chǎn)量 N01 = 計(jì)算方法實(shí)驗(yàn)報(bào)告 20 b) 三次樣條插值實(shí)驗(yàn)結(jié)果 三次樣條插值 2020 年的產(chǎn)量及插值函數(shù) f = (67507170809569*t)/21990232555520+((6226409313348779*t)/2161727821137838080417169423994368193/72057594037927936)*(t2020)^212853306597599315/2199023255552 f0 = 三次樣條插值 2020 年的產(chǎn)量 f1 = 結(jié)果分析 由實(shí)驗(yàn)結(jié)果可以看出， 在對(duì) 2020 年產(chǎn)量計(jì)算的時(shí)候， Newton 插值出現(xiàn)了很明顯的偏差，而三次樣條插值表現(xiàn)的很好。與 Newton 多項(xiàng)式插值相比較 三次樣條 插值顯然得到的結(jié)果更理想，這是由于多項(xiàng)式插值，當(dāng) n 較大時(shí)，插值結(jié)果可能出現(xiàn)不一致收斂的 現(xiàn)象，并且高次插值計(jì)算復(fù)雜，此外三次樣條插值還具有光滑性的優(yōu)點(diǎn)，因此插值計(jì)算時(shí)應(yīng)盡量應(yīng)該三次樣條的方法。 計(jì)算方法實(shí)驗(yàn)報(bào)告 21 問題 5 五、 假定某天的氣溫變化記錄如下表所示，試用最小二乘法找出這一天的氣溫變化的規(guī)律 ,試計(jì)算這一天的平均氣溫。 時(shí)刻 0 1 2 3 4 5 6 7 8 9 10 11 12 平均氣溫 15 14 14 14 14 15 16 18 20 20 23 25 28 時(shí)刻 13 14 15 16 17 18 19 20 21 22 23 24 平均氣溫 31 32 31 29 27 25 24 22 20 18 17 16 考慮使用二次函數(shù)、三次函數(shù)、四次函數(shù)以及指數(shù)型的函數(shù) 2()b t cC ae??? ,計(jì)算相應(yīng)的系數(shù)，估算誤差，并作圖比較各種函數(shù)之間的區(qū)別。 問題分析 從整體上考慮近似函數(shù) 同所給數(shù)據(jù)點(diǎn) (i=0,1,? ,m) 誤差 (i=0,1,? ,m) (i=0,1, ? ,m) 絕對(duì)值的最大值 ，即誤差 向量的∞ —范數(shù)；二是誤差絕對(duì)值的和 ，即誤差向量 r 的 1—范數(shù)；三是誤差平方和 的算術(shù)平方根，即誤差向量 r 的 2—范數(shù)；前兩種方法簡(jiǎn)單、自然，但不便于微分運(yùn)算 ，后一種方法相當(dāng)于考慮 2—范數(shù)的平方，因此在曲線擬合中常采用誤差平方和 來 度量誤差 (i=0， 1，?， m)的整計(jì)算方法實(shí)驗(yàn)報(bào)告 22 數(shù)據(jù)擬合的具體作法是：對(duì)給定數(shù)據(jù) (i=0,1,?， m)，在取定的函數(shù)類 中 ,求 ,使誤差 (i=0,1,? ,m)的平方和最小，即 = 從幾何意義上講，就是尋求與給定點(diǎn) (i=0,1,? ,m)的距離平方和為最。函數(shù) 稱為擬合 函數(shù)或最小二乘解，求擬合函數(shù)的方法稱為曲線擬合的最小二乘法。 在曲線擬合中，函數(shù)類 可有不同的選取方法 ，如多項(xiàng)式擬合算法。 假設(shè)給定數(shù)據(jù)點(diǎn) (i=0,1,? ,m)， 為所有次數(shù)不超過 的多項(xiàng)式構(gòu)成的函數(shù)類，現(xiàn)求一 ,使得 (1) 當(dāng)擬合函數(shù)為多項(xiàng)式時(shí)，稱為多項(xiàng)式擬合，滿足式（ 1）的 稱為最小二乘擬合多項(xiàng)式。特別地，當(dāng) n=1 為 的多元函數(shù)，因此上述問題即為求 的極值 問題。由多元函數(shù)求極值的必要條件，得 (2) 即 (3) （ 3）是關(guān)于 的線性方程組，用矩陣表示為 計(jì)算方法實(shí)驗(yàn)報(bào)告 23 (4) 式（ 3）或式（ 4 可以證明，方程組（ 4）的系數(shù)矩陣是一個(gè)對(duì)稱正定矩陣，故存在唯一解。從式（ 4）中解出 (k=0,1,?， n) (5) 可以證明，式（ 5）中的 滿足式（ 1），即 為所求的擬合多項(xiàng)式。我們把 稱為最小二乘擬合多項(xiàng)式 由式 (2) (6) 求解 的方法主要有下面兩種 ， 本算法 并分別 采用 兩種方法求解。 i. 求解法方程法 如果 直接利用（ 4）進(jìn)行求解解方程，解法 比較簡(jiǎn)單， 但是一方面運(yùn)算量較大，并且 由于條件數(shù)的影響 可能會(huì)有較大誤差。 本設(shè)計(jì)采用了高斯消去法，減少了誤差。 ii. 正交化法 對(duì)生成的矩陣進(jìn)行 QR 分解， 上式（ 4）可寫成矩陣形式， GTGa = GTy，利用矩陣的 QR 分解，可以 G = QT (RO)，（ QR 分解算法 具體實(shí)現(xiàn)過程見附錄中的Matlab 程序 ）， Qy = (h1h2)，因此 Ra=h1，利用解方程的回代公式可以求得向量 a。計(jì)算方法實(shí)驗(yàn)報(bào)告 24 具體實(shí)現(xiàn)程序見附錄 5. 用 2()b t cC ae??? 擬合時(shí)，可以考慮兩邊取對(duì)數(shù) log(C) = log(a) ?bt2 +2bct?c2，然后對(duì) (c,log(C))按二次多項(xiàng)式擬合。 Matlab 求解程序見附錄五。 實(shí)驗(yàn)結(jié)果 a) 法方程解法實(shí)驗(yàn)結(jié)果 擬合二次函數(shù)如下及相應(yīng)的系數(shù) p = Px = (1238*x^2)/13455 + (2882931353438435*x)/1125899906842624 + 986/117 估算誤差 e = 擬合三次函數(shù)如下及相應(yīng)的系數(shù) p = Px=(4727727282891159*x^3)/576460752303423488+(6016*x^2)/29601 (487295332134639*x)/2251799813685248 + 13072/975 估算誤差 e = 擬合四次函數(shù)如下及相應(yīng)的系數(shù) p = Px=(8277009104528411*x^4)/9223372036854775808 (7389688649119097*x^3)/144115188075855872+ (3869637877579055*x^2)/4503599627370496 (8004500866520783*x)/2251799813685248 + 94306/5655 估算誤差 e = 計(jì)算方法實(shí)驗(yàn)報(bào)告 25 擬合指數(shù)函數(shù)如下的系數(shù) p = b = c = a = 估算誤差 e = b) 正交分解法 擬合二次函數(shù)如下及相應(yīng)的系數(shù) p1 = Px1 = (1238*x^2)/13455 + (5765862706876871*x)/2251799813685248 + 986/117 估算誤差 e1 = 擬合三次函數(shù)如下及相應(yīng)的系數(shù) p2 = Px2=(2363863641445577*x^3)/288230376151711744+(6016*x^2)/29601 計(jì)算方法實(shí)驗(yàn)報(bào)告 26 (1949181328538565*x)/9007199254740992 + 13072/975 估算誤差 e2 = 擬合四次函數(shù)如下及相應(yīng)的系數(shù) p3 = Px3 = (8277009104528155*x^4)/9223372036854775808 (3694844324559451*x^3)/72057594037927936+ (7739275755157927*x^2)/9007199254740992 (8004500866520577*x)/2251799813685248 + 94306/5655 估算誤差 e3 = 擬合指數(shù)函數(shù)的系數(shù)如下 p4 = b = c = a = 估算誤差 e4 = 計(jì)算方法實(shí)驗(yàn)報(bào)告 27 結(jié)果分析 由結(jié)果可以看出，隨著擬合次數(shù)的增多誤差有明顯減少的趨勢(shì)，指數(shù)函數(shù)擬合跟二次函數(shù)擬合曲線很接近，但相對(duì)二次函數(shù)擬合誤差減少了，四次函數(shù)擬合誤差最小。 計(jì)算方法實(shí)驗(yàn)報(bào)告 28 附錄 附錄一 問題 1Matlab 求解 %計(jì)算 S=1/(1^2)+1/(2^2)+1/(3^2)+...+1/(100000^2),誤差小于 1*10e6 clc。 clear。 sum1=0。 a=1。 %估算截取的項(xiàng)數(shù)，使誤差小于 1*10e6 while ((1/a1/100000)(1/(a+1)1/100001))/21e6 a=a+1。 end %限定誤差 disp(39。最大誤差 :39。) e1=((1/(a1)1/100000)(1/(a)1/100001))/2 %截取的項(xiàng)數(shù) disp(39。截取的項(xiàng)數(shù) :39。) b=a1 %按從小到大的順序相加 for k=b:1:1 sum1=sum1+1/(k.^2)。 end sum1=sum1+((1/(b)1/100000)+(1/(b+1)1/100001))/2。 disp(39。估算值 :39。) fprintf(39。%\n\n39。,sum1) %準(zhǔn)確值，為了證明方法的可靠性 sum=0。 for k=100000:1:1 sum=sum+1/(k.^2)。 end disp(39。準(zhǔn)確值 :39。) fprintf(39。%\n39。,sum) %實(shí)際誤差 disp(39。實(shí)際誤差 :39。) ess=abs((sumsum1))/sum 計(jì)算方法實(shí)驗(yàn)報(bào)告 29 附錄二 問題 2Matlab 求解 %對(duì)非奇異矩陣 a 作 lu 分解 clc clear disp(39。請(qǐng)輸入矩陣 A39。) A=input(39。請(qǐng)輸入系數(shù)矩陣 A=:39。) %判斷被分解矩陣的維數(shù) [n,n]=size(A)。 %給出 U 的第一行 L 的第一列 for j=1:n U(1,j)=A(1,j)。 L(j,1)=A(j,1)/U(1,1)。 end i=1。 for k=2:n for j=k:n for m=1:(k1) %計(jì)算公式中對(duì)應(yīng)的 l*u 的和 uu(m)=L(k,m)*U(m,j)。 end %u 的分解結(jié)果 U