簡(jiǎn)介廣東工業(yè)大學(xué)數(shù)值分析課程實(shí)驗(yàn)四的Python實(shí)驗(yàn)報(bào)告完整覆蓋超定方程組最小二乘解與最小二乘擬合兩大核心實(shí)驗(yàn)。報(bào)告為docx文檔共1個(gè)文件、112KB包含題目三、題目六及石油產(chǎn)量預(yù)測(cè)三道典型問題均附有可運(yùn)行的Python源碼、算法介紹和結(jié)果分析。題目三基于最小二乘解的充要條件構(gòu)造正規(guī)方程組ATAxATb借助numpy的線性代數(shù)求解器高效算出兩組超定方程的解題目六以1和1/x為基函數(shù)建立線性組合模型調(diào)用sklearn的LinearRegression擬合觀察數(shù)據(jù)得到形如y2.01431.00931/x的擬合函數(shù)石油產(chǎn)量預(yù)測(cè)則演示了將數(shù)值方法應(yīng)用于實(shí)際數(shù)據(jù)建模的完整流程。讀者可從中學(xué)習(xí)最小二乘問題的矩陣推導(dǎo)、正規(guī)方程構(gòu)造、Python庫調(diào)用及結(jié)果分析方法適合正在學(xué)習(xí)數(shù)值分析、需要參考完整實(shí)驗(yàn)報(bào)告或復(fù)習(xí)相關(guān)算法的高校學(xué)生。目前已有212人學(xué)習(xí)該資源。1. 某高校Python數(shù)值分析實(shí)驗(yàn)四為什么代碼三分鐘跑完作業(yè)卻要寫一整晚數(shù)值分析這門課真正動(dòng)手做起來和聽課時(shí)完全是兩個(gè)世界。某高校的Python數(shù)值分析實(shí)驗(yàn)四表面上是插值和擬合兩個(gè)經(jīng)典主題——拉格朗日插值、牛頓插值、最小二乘擬合代碼用NumPy就能寫完跑通不過三分鐘。但真正消耗時(shí)間的是把課本上的公式變成能交代得清楚的程序把結(jié)果畫成能放進(jìn)實(shí)驗(yàn)報(bào)告的圖表以及處理那些課本沒寫、IDE也不會(huì)告訴你的邊界情況。這篇筆記按我實(shí)際完成實(shí)驗(yàn)四的順序來寫從任務(wù)拆解、環(huán)境準(zhǔn)備、算法實(shí)現(xiàn)到避坑排查每一步都給可復(fù)現(xiàn)的代碼和參數(shù)說明。適合第一次做數(shù)值分析實(shí)驗(yàn)的學(xué)生也適合想快速拾起插值與擬合實(shí)現(xiàn)的工程師。2. 拆開實(shí)驗(yàn)四的作業(yè)要求任務(wù)清單、環(huán)境準(zhǔn)備與報(bào)告結(jié)構(gòu)2.1 先搞清實(shí)驗(yàn)四考什么兩個(gè)主線任務(wù)和三類題型數(shù)值分析實(shí)驗(yàn)四在多數(shù)高校的教學(xué)安排里主題都落在“插值法與曲線擬合”這一章。拆開來看主線任務(wù)基本是固定的用Python實(shí)現(xiàn)拉格朗日插值和牛頓插值再用最小二乘法做曲線擬合最后把結(jié)果和理論做對(duì)比。不同年份、不同老師的題目形式會(huì)有變化但底層要求逃不出三類題型。第一類是給定函數(shù)或數(shù)據(jù)表要求計(jì)算指定點(diǎn)的插值結(jié)果第二類是給一組觀測(cè)數(shù)據(jù)要求用多項(xiàng)式擬合出趨勢(shì)第三類是綜合題比如對(duì)某個(gè)病態(tài)函數(shù)做插值并討論誤差或者用擬合預(yù)測(cè)下一組數(shù)據(jù)點(diǎn)。先確認(rèn)自己拿到的是哪類題再?zèng)Q定代碼怎么組織能省下大量返工時(shí)間。我見過不少同學(xué)上來就寫拉格朗日的循環(huán)寫完才發(fā)現(xiàn)題目要求的是擬合。更常見的坑是題目給了十個(gè)數(shù)據(jù)點(diǎn)但要求分別用三次、五次、七次多項(xiàng)式做擬合這時(shí)候如果沒提前把擬合邏輯封裝成函數(shù)后面每一個(gè)次數(shù)都要復(fù)制粘貼改參數(shù)實(shí)驗(yàn)報(bào)告的代碼部分會(huì)顯得很亂查錯(cuò)也麻煩。所以第一步不是寫代碼是把題目要求翻譯成一個(gè)清單要實(shí)現(xiàn)的算法有幾個(gè)、輸入輸出是什么、要求畫幾張圖、誤差分析要算哪些量。2.2 實(shí)驗(yàn)環(huán)境準(zhǔn)備conda創(chuàng)建虛擬環(huán)境與依賴清單實(shí)驗(yàn)四的代碼依賴面不寬但版本不一致會(huì)在后面圖表上翻車。我一般會(huì)先建一個(gè)干凈的虛擬環(huán)境避免和機(jī)器上其他項(xiàng)目共用NumPy和SciPy。用conda創(chuàng)建虛擬環(huán)境的命令如下conda create -n numeric-analysis python3.11 conda activate numeric-analysis pip install numpy matplotlib scipy選Python 3.11而不是最新版本是因?yàn)镹umPy和SciPy對(duì)最新Python版本的支持往往滯后幾個(gè)月裝完才發(fā)現(xiàn)某個(gè)依賴編譯失敗再回退更浪費(fèi)時(shí)間。NumPy負(fù)責(zé)數(shù)組計(jì)算和多項(xiàng)式運(yùn)算SciPy在擬合部分會(huì)用到curve_fitMatplotlib負(fù)責(zé)畫圖。這三個(gè)庫覆蓋面足夠pandas在這個(gè)實(shí)驗(yàn)里用不上裝不裝看個(gè)人習(xí)慣。依賴安裝完以后建議順手驗(yàn)證一下版本號(hào)是否匹配尤其是NumPy和SciPy之間的版本兼容。我遇到過NumPy 1.x和SciPy 1.11搭配時(shí)某個(gè)接口行為變化導(dǎo)致polyfit返回結(jié)果格式不同。驗(yàn)證方式很簡(jiǎn)單在Python交互環(huán)境里分別import并打印版本號(hào)確認(rèn)能正常加載再開始寫實(shí)驗(yàn)代碼。如果用的是機(jī)房預(yù)裝的環(huán)境版本往往偏舊建議先確認(rèn)能不能裝新包不能裝就基于現(xiàn)有版本調(diào)整寫法別到最后一刻才暴露問題。2.3 實(shí)驗(yàn)報(bào)告的結(jié)構(gòu)算法原理、程序?qū)崿F(xiàn)、結(jié)果分析三位一體實(shí)驗(yàn)報(bào)告怎么寫決定了這份作業(yè)的最終成績(jī)。數(shù)值分析實(shí)驗(yàn)的評(píng)分點(diǎn)通常有三個(gè)算法原理是否說清楚、代碼是否能跑出正確結(jié)果、結(jié)果分析是否有深度。三者缺一不可但很多同學(xué)只重視代碼原理部分抄課本分析部分寫“圖像符合預(yù)期”這等于把最容易拿分的部分放棄了。我的習(xí)慣是先在報(bào)告里畫一個(gè)整體框架再按照框架去補(bǔ)內(nèi)容。第一部分寫算法原理包括公式推導(dǎo)和算法步驟不需要長篇大論但要寫清楚公式中每個(gè)符號(hào)的含義比如插值節(jié)點(diǎn)、基函數(shù)、差商。第二部分寫程序?qū)崿F(xiàn)代碼按功能拆分關(guān)鍵函數(shù)要給注釋。第三部分寫結(jié)果分析這一部分要包含數(shù)值結(jié)果和圖像的對(duì)比以及誤差分析——插值法算出的值和理論值差多少、擬合的殘差有多大、為什么會(huì)有這個(gè)誤差。整份報(bào)告的邏輯是讀者只看報(bào)告不看你跑代碼也能復(fù)現(xiàn)你的結(jié)果這才是一份合格的實(shí)驗(yàn)報(bào)告。3. 拉格朗日插值與牛頓插值從公式到Python實(shí)現(xiàn)3.1 拉格朗日插值法的數(shù)學(xué)原理與代碼實(shí)現(xiàn)插值問題的本質(zhì)是已知n1個(gè)數(shù)據(jù)點(diǎn)找一條通過所有這些點(diǎn)的曲線。拉格朗日插值的思想最直白構(gòu)造一組基函數(shù)每個(gè)基函數(shù)在某個(gè)數(shù)據(jù)點(diǎn)上取值為1在其他數(shù)據(jù)點(diǎn)上取值為0然后把所有基函數(shù)加權(quán)疊加。權(quán)重就是對(duì)應(yīng)數(shù)據(jù)點(diǎn)的函數(shù)值。這樣做的好處是公式對(duì)稱、好理解代碼實(shí)現(xiàn)也不需要額外存儲(chǔ)中間量。import numpy as np def lagrange_interp(x_points, y_points, x_eval): 拉格朗日插值計(jì)算單個(gè)點(diǎn) x_eval 處的插值結(jié)果 x_points: 插值節(jié)點(diǎn)數(shù)組例如 [1, 2, 3] y_points: 對(duì)應(yīng)的函數(shù)值數(shù)組例如 [1, 4, 9] x_eval: 需要求插值的點(diǎn)可以是標(biāo)量或數(shù)組 n len(x_points) if n 0: raise ValueError(插值節(jié)點(diǎn)不能為空) total 0.0 for i in range(n): term y_points[i] for j in range(n): if j ! i: term * (x_eval - x_points[j]) / (x_points[i] - x_points[j]) total term return total這段代碼的邏輯是兩層循環(huán)外層遍歷每個(gè)數(shù)據(jù)點(diǎn)內(nèi)層累乘計(jì)算對(duì)應(yīng)基函數(shù)在當(dāng)前求值點(diǎn)x_eval處的值。內(nèi)層循環(huán)里分母是x_points[i] - x_points[j]當(dāng)x_eval恰好等于某個(gè)插值節(jié)點(diǎn)時(shí)對(duì)應(yīng)項(xiàng)的分母會(huì)精確抵消成一其他項(xiàng)因?yàn)榉肿映霈F(xiàn)零而整個(gè)為0最終返回該節(jié)點(diǎn)的函數(shù)值這是拉格朗日插值最基本的性質(zhì)。代碼里加了個(gè)空節(jié)點(diǎn)檢查防止傳入空數(shù)組時(shí)報(bào)出難以理解的numpy錯(cuò)誤。用的時(shí)候注意兩個(gè)參數(shù)語義x_points是插值節(jié)點(diǎn)而非插值區(qū)間端點(diǎn)它決定多項(xiàng)式的次數(shù)n個(gè)節(jié)點(diǎn)對(duì)應(yīng)n-1次多項(xiàng)式。x_eval可以是單個(gè)數(shù)值也可以是numpy數(shù)組函數(shù)內(nèi)部會(huì)自動(dòng)做廣播。實(shí)測(cè)中對(duì)數(shù)組求值時(shí)這個(gè)函數(shù)會(huì)在每個(gè)元素上獨(dú)立執(zhí)行兩層循環(huán)數(shù)據(jù)量上千點(diǎn)時(shí)速度會(huì)明顯下降所以如果只需要在特定幾個(gè)點(diǎn)求值不要傳一個(gè)很長的密集數(shù)組這屬于代碼效率和精度無關(guān)但和體驗(yàn)相關(guān)的問題。3.2 牛頓插值法差商表與遞推實(shí)現(xiàn)牛頓插值和拉格朗日插值數(shù)學(xué)上等價(jià)但實(shí)現(xiàn)思路不同。牛頓插值通過構(gòu)造差商表來逐步構(gòu)建多項(xiàng)式它的優(yōu)勢(shì)在于當(dāng)新增一個(gè)插值節(jié)點(diǎn)時(shí)不需要重寫所有基函數(shù)只需在原有多項(xiàng)式上加一項(xiàng)。這個(gè)性質(zhì)在實(shí)驗(yàn)里體現(xiàn)得最明顯的是題目要求先取前三個(gè)點(diǎn)插值看結(jié)果再取前四個(gè)點(diǎn)插值看結(jié)果用牛頓插值可以復(fù)用之前的計(jì)算。def newton_divided_differences(x_points, y_points): 計(jì)算差商表返回對(duì)角線上的差商系數(shù)數(shù)組 n len(x_points) coef np.array(y_points, dtypefloat) for j in range(1, n): for i in range(n - 1, j - 1, -1): coef[i] (coef[i] - coef[i - 1]) / (x_points[i] - x_points[i - j]) return coef def newton_interp(coef, x_points, x_eval): 用差商系數(shù)計(jì)算牛頓插值多項(xiàng)式在 x_eval 處的值 coef: newton_divided_differences 的返回值 n len(coef) result coef[0] for i in range(1, n): term coef[i] for j in range(i): term * (x_eval - x_points[j]) result term return result差商表計(jì)算是原地更新的coef數(shù)組在每一輪外層循環(huán)中被覆蓋最終coef[i]存放的是第i階差商。注意內(nèi)層循環(huán)從后向前遍歷這保證計(jì)算新值時(shí)用到的是上一輪的結(jié)果而不是本輪已經(jīng)覆蓋的值順序反了整個(gè)差商表就是錯(cuò)的。牛頓插值的求值過程用到了嵌套乘法從最高次項(xiàng)開始逐層提取公因子這樣比直接展開多項(xiàng)式運(yùn)算量更小也更符合數(shù)值分析教材里對(duì)秦九韶算法的介紹。實(shí)操中差商表函數(shù)只在構(gòu)造階段調(diào)用一次之后對(duì)多個(gè)x_eval求值時(shí)只需要反復(fù)調(diào)用newton_interp。如果你要插值的點(diǎn)有幾百上千個(gè)這個(gè)做法的優(yōu)勢(shì)非常明顯——差商構(gòu)造是O(n2)而每個(gè)點(diǎn)求值是O(n)。拉格朗日插值對(duì)每個(gè)求值點(diǎn)都要重新計(jì)算所有基函數(shù)同樣是O(n2)但乘在點(diǎn)的數(shù)量上數(shù)據(jù)量一大差距就拉開了。3.3 用matplotlib畫出插值效果對(duì)比圖實(shí)驗(yàn)報(bào)告里最直觀的內(nèi)容是一張插值效果圖原始數(shù)據(jù)點(diǎn)、拉格朗日插值曲線、牛頓插值曲線在同一坐標(biāo)系中展示。畫圖代碼不長但有兩個(gè)細(xì)節(jié)要注意插值曲線的橫坐標(biāo)要足夠密不然曲線會(huì)顯得很生硬縱坐標(biāo)范圍要留出余量避免曲線在邊緣處超出畫布。import matplotlib.pyplot as plt x_points np.array([0.0, 1.0, 2.0, 3.0, 4.0]) y_points np.array([1.0, 2.0, 1.5, 3.0, 2.8]) x_dense np.linspace(0, 4, 200) fig, ax plt.subplots(figsize(8, 5)) ax.plot(x_points, y_points, o, label原始數(shù)據(jù)點(diǎn)) ax.plot(x_dense, lagrange_interp(x_points, y_points, x_dense), -, label拉格朗日插值) coef newton_divided_differences(x_points, y_points) ax.plot(x_dense, newton_interp(coef, x_points, x_dense), --, label牛頓插值) ax.set_xlabel(x) ax.set_ylabel(y) ax.legend() ax.grid(True, linestyle--, alpha0.6) plt.show()畫圖時(shí)x_dense用np.linspace(0, 4, 200)生成200個(gè)均勻分布的點(diǎn)具體點(diǎn)數(shù)取決于插值多項(xiàng)式的次數(shù)和區(qū)間跨度一般200到500足夠。如果兩條插值曲線完全重合說明兩個(gè)算法在精度上一致這是正常的如果出現(xiàn)細(xì)微差異優(yōu)先檢查差商表的迭代順序是否正確。對(duì)比圖存在的意義是驗(yàn)證兩個(gè)算法結(jié)果一致順帶檢查曲線的形狀是否符合直覺。4. 最小二乘曲線擬合把公式擬合到數(shù)據(jù)點(diǎn)上的完整流程4.1 最小二乘原理與numpy.polyfit最小實(shí)現(xiàn)擬合和插值的核心區(qū)別在于插值要求曲線嚴(yán)格穿過每個(gè)數(shù)據(jù)點(diǎn)擬合只要求曲線在整體上最接近數(shù)據(jù)點(diǎn)。最小二乘的思想是找一組參數(shù)讓模型計(jì)算值與觀測(cè)值的平方誤差之和最小。平方誤差比絕對(duì)誤差更敏感而且方便求導(dǎo)所以成為最常用的擬合準(zhǔn)則。import numpy as np x_data np.array([0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0]) y_data np.array([0.1, 0.8, 1.8, 2.9, 5.1, 7.2, 10.0]) degree 2 coeffs np.polyfit(x_data, y_data, degree) p np.poly1d(coeffs) print(擬合系數(shù), coeffs) print(擬合多項(xiàng)式, p)np.polyfit返回的是多項(xiàng)式系數(shù)從最高次到常數(shù)項(xiàng)排列系數(shù)順序和poly1d的輸入是對(duì)應(yīng)的。degree參數(shù)是最關(guān)鍵的調(diào)試項(xiàng)取1時(shí)是線性擬合取2時(shí)是拋物線擬合取大了會(huì)出現(xiàn)過擬合——曲線為了貼近個(gè)別數(shù)據(jù)點(diǎn)而劇烈抖動(dòng)尤其是在數(shù)據(jù)本身不是高階多項(xiàng)式的情況下。判斷次數(shù)是否合適的標(biāo)準(zhǔn)很簡(jiǎn)單看擬合曲線走勢(shì)是否符合數(shù)據(jù)的物理含義再看殘差是否隨機(jī)分布而不是呈現(xiàn)系統(tǒng)性偏差。poly1d是把系數(shù)數(shù)組包裝成可以直接調(diào)用的多項(xiàng)式對(duì)象調(diào)用p(x_eval)就能批量求值。這里有個(gè)小細(xì)節(jié)polyfit內(nèi)部用的是基于SVD的穩(wěn)定算法對(duì)于次數(shù)不超過8的多項(xiàng)式擬合很少出現(xiàn)數(shù)值問題但超過10次就要警惕矩陣病態(tài)結(jié)果可能不可信。所以實(shí)驗(yàn)題里如果要求擬合次數(shù)特別高建議檢查一下擬合系數(shù)是否出現(xiàn)數(shù)量級(jí)異常的大值那通常是條件數(shù)過大導(dǎo)致的數(shù)值性翻車。4.2 自定義函數(shù)擬合scipy.optimize.curve_fit的使用實(shí)驗(yàn)題里有時(shí)候不要求多項(xiàng)式擬合而是指定一個(gè)函數(shù)形式比如指數(shù)衰減模型或者帶偏移的正弦函數(shù)。np.polyfit只支持多項(xiàng)式這時(shí)候要換用scipy的curve_fit它可以擬合任意你寫出的Python函數(shù)。核心用法是先定義一個(gè)函數(shù)形式函數(shù)的第一個(gè)參數(shù)是自變量后面是待擬合參數(shù)。from scipy.optimize import curve_fit def exp_model(x, a, b, c): 指數(shù)衰減加常數(shù)偏移的擬合模型 return a * np.exp(-b * x) c popt, pcov curve_fit(exp_model, x_data, y_data, p0[1.0, 1.0, 0.0]) a_opt, b_opt, c_opt popt print(f擬合參數(shù): a{a_opt:.4f}, b{b_opt:.4f}, c{c_opt:.4f}) x_smooth np.linspace(x_data.min(), x_data.max(), 100) y_fit exp_model(x_smooth, *popt)p0是參數(shù)初值這個(gè)參數(shù)在很多時(shí)候決定了擬合能不能收斂。初值差得太遠(yuǎn)算法可能陷入局部極小值甚至直接報(bào)錯(cuò)。給p0的經(jīng)驗(yàn)法則是先看數(shù)據(jù)點(diǎn)的取值范圍大致估計(jì)一下參數(shù)的量級(jí)。比如數(shù)據(jù)y值在0到10之間c的初值可以設(shè)為數(shù)據(jù)中點(diǎn)或最小值附近a的初值設(shè)為y的跨度b的初值根據(jù)曲線衰減的快慢估一個(gè)正數(shù)。如果你發(fā)現(xiàn)擬合結(jié)果明顯不合理最先該改的是p0而不是模型本身。curve_fit返回的pcov是參數(shù)協(xié)方差矩陣對(duì)角線元素開根號(hào)是參數(shù)的標(biāo)準(zhǔn)差代表擬合參數(shù)的不確定度。實(shí)驗(yàn)報(bào)告里如果想體現(xiàn)嚴(yán)謹(jǐn)性可以把參數(shù)誤差也打印出來。當(dāng)協(xié)方差矩陣?yán)锏闹荡蟮诫x譜時(shí)往往說明模型參數(shù)之間存在強(qiáng)相關(guān)性或者數(shù)據(jù)量不足以確定這么多參數(shù)這時(shí)候的擬合系數(shù)雖然能算出數(shù)但不能當(dāng)作物理解釋使用。這個(gè)點(diǎn)算是一個(gè)隱藏的知識(shí)點(diǎn)寫在報(bào)告的分析部分會(huì)很加分。4.3 擬合效果的量化評(píng)估殘差、R方與可視化擬合做完不評(píng)估就寫報(bào)告等于黑匣子輸出直接交差。評(píng)估擬合質(zhì)量最常用的是三個(gè)量殘差、均方根誤差和R方。殘差是每個(gè)數(shù)據(jù)點(diǎn)的觀測(cè)值減擬合值R方表示擬合曲線解釋了數(shù)據(jù)方差的百分比越接近1越好。residuals y_data - exp_model(x_data, *popt) ss_res np.sum(residuals**2) ss_tot np.sum((y_data - np.mean(y_data))**2) r_squared 1 - ss_res / ss_tot rmse np.sqrt(np.mean(residuals**2)) print(fR方 {r_squared:.4f}) print(f均方根誤差 RMSE {rmse:.4f})注意R方有一個(gè)陷阱當(dāng)模型形式比較復(fù)雜、參數(shù)很多時(shí)R方會(huì)虛高甚至接近1但這不代表模型真的好。一個(gè)十次多項(xiàng)式擬合十個(gè)數(shù)據(jù)點(diǎn)R方必然是1但那完全沒有預(yù)測(cè)意義。所以評(píng)估時(shí)要同時(shí)看殘差圖——如果殘差隨x增大而增大說明模型形式不對(duì)如果殘差呈現(xiàn)波浪形說明漏掉了周期性分量。這些判斷是數(shù)值分析實(shí)驗(yàn)最有價(jià)值的訓(xùn)練比單純算對(duì)拉格朗日插值結(jié)果重要得多。把殘差圖畫出來和擬合效果圖畫在一起報(bào)告的分析部分立刻有了深度。5. 實(shí)驗(yàn)四避坑指南5個(gè)最常見的翻車點(diǎn)與排查方法5.1 拉格朗日插值次數(shù)過高導(dǎo)致龍格現(xiàn)象數(shù)據(jù)越多結(jié)果越離譜現(xiàn)象用拉格朗日插值去逼近某些連續(xù)函數(shù)插值節(jié)點(diǎn)增加到十幾個(gè)甚至更多插值曲線在區(qū)間兩端劇烈震蕩誤差反而比少用幾個(gè)節(jié)點(diǎn)時(shí)更大。原因這是數(shù)值分析教材里反復(fù)強(qiáng)調(diào)的龍格現(xiàn)象。拉格朗日插值多項(xiàng)式在等距節(jié)點(diǎn)下對(duì)龍格函數(shù)、絕對(duì)值函數(shù)這類目標(biāo)做高次逼近時(shí)插值余項(xiàng)在端點(diǎn)附近會(huì)放大多項(xiàng)式次數(shù)越高震蕩越明顯。解決明確實(shí)驗(yàn)四的題目通常只要求4到6個(gè)節(jié)點(diǎn)的插值不需要追求高次。如果題目要求對(duì)比不同節(jié)點(diǎn)數(shù)下的誤差把插值點(diǎn)數(shù)從3逐個(gè)增加到10畫誤差曲線讓讀者直觀看到誤差先降后升的過程這本身就是一份合格的現(xiàn)象分析。如果實(shí)際工程中遇到類似需求應(yīng)該改用分段插值或樣條插值而不是硬上高次多項(xiàng)式。5.2 圖表的橫縱坐標(biāo)出現(xiàn)方塊字中文標(biāo)注全部亂碼現(xiàn)象matplotlib畫圖時(shí)標(biāo)題和圖例里的中文顯示成一個(gè)個(gè)小方塊英文和數(shù)字正常實(shí)驗(yàn)報(bào)告截圖里一眼就能看出問題。原因matplotlib默認(rèn)字體不包含中文字形需要顯式指定中文字體。Windows系統(tǒng)一般有SimHei或Microsoft YaHeiLinux和macOS的字體名不同用系統(tǒng)自帶字體名逐個(gè)嘗試。解決在畫圖代碼最前面加上字體設(shè)置讓matplotlib全局使用中文字體如下所示import matplotlib.pyplot as plt plt.rcParams[font.sans-serif] [SimHei, Microsoft YaHei, PingFang SC] plt.rcParams[axes.unicode_minus] False第二行去掉坐標(biāo)軸負(fù)號(hào)的特殊處理是因?yàn)榍袚Q了中文字體之后負(fù)號(hào)有可能顯示成一個(gè)奇怪的方框。這個(gè)問題屬于程序能跑但報(bào)告不能看的隱形翻車查錯(cuò)成本高提前設(shè)置好就能完全避免。5.3 差商表計(jì)算順序?qū)懛磁ｎD插值結(jié)果對(duì)不上拉格朗日現(xiàn)象牛頓插值和拉格朗日插值結(jié)果不一致差距可能很小也可能很大越小規(guī)模的數(shù)據(jù)差距越小很多同學(xué)因此沒發(fā)現(xiàn)。原因newton_divided_differences函數(shù)里內(nèi)層循環(huán)的方向錯(cuò)了。差商表必須從后向前原地更新如果寫成從前向后計(jì)算新的差商時(shí)用的是已經(jīng)被本輪覆蓋過的值導(dǎo)致差商全部錯(cuò)誤。數(shù)據(jù)點(diǎn)少時(shí)誤差不明顯數(shù)據(jù)點(diǎn)多時(shí)完全對(duì)不上。解決內(nèi)層循環(huán)固定使用for i in range(n - 1, j - 1, -1)從后往前遍歷。排查思路也很簡(jiǎn)單用三個(gè)點(diǎn)的手工差商計(jì)算結(jié)果和Python輸出對(duì)比如果差商對(duì)不上問題就在循環(huán)方向。這類錯(cuò)誤的特點(diǎn)是代碼能跑、結(jié)果貼近但不對(duì)最容易在實(shí)驗(yàn)報(bào)告的誤差分析部分暴露靠直覺很難發(fā)現(xiàn)。5.4 擬合結(jié)果出現(xiàn)巨大系數(shù)多項(xiàng)式次數(shù)太高導(dǎo)致矩陣病態(tài)現(xiàn)象polyfit指定degree9或更高時(shí)擬合系數(shù)出現(xiàn)10的6次方以上的大數(shù)擬合曲線在數(shù)據(jù)點(diǎn)之間劇烈波動(dòng)部分區(qū)間甚至超出數(shù)據(jù)范圍好幾倍。原因高次多項(xiàng)式擬合的等效設(shè)計(jì)矩陣條件數(shù)隨次數(shù)指數(shù)增長數(shù)據(jù)點(diǎn)x跨度大時(shí)尤其嚴(yán)重。矩陣接近奇異SVD分解得到的解被微小舍入誤差放大系數(shù)自然不穩(wěn)定。解決首先判斷題目是否真的要求高次擬合。如果要求做不同次數(shù)擬合的對(duì)比把次數(shù)范圍控制在1到6之間重點(diǎn)觀察擬合誤差變化趨勢(shì)。如果必須擬合高次多項(xiàng)式考慮先對(duì)x數(shù)據(jù)做歸一化把數(shù)據(jù)平移到零均值并縮放到單位方差再擬合最后把系數(shù)還原回原始坐標(biāo)。歸一化這步在實(shí)驗(yàn)報(bào)告里是加分項(xiàng)能體現(xiàn)你對(duì)數(shù)值穩(wěn)定性的理解。5.5 求值點(diǎn)落在插值節(jié)點(diǎn)范圍之外結(jié)果完全失控現(xiàn)象用插值函數(shù)計(jì)算區(qū)間外的點(diǎn)計(jì)算結(jié)果非常大和真實(shí)函數(shù)值差出幾個(gè)數(shù)量級(jí)有的同學(xué)以為是代碼精度不夠。原因插值多項(xiàng)式只在節(jié)點(diǎn)覆蓋的區(qū)間內(nèi)有可靠逼近能力。區(qū)間外沒有任何數(shù)據(jù)點(diǎn)約束多項(xiàng)式形狀高次項(xiàng)起主導(dǎo)作用結(jié)果必然外溢。這是插值方法的固有屬性不是bug。解決求值前先檢查x_eval是否在x_points的范圍內(nèi)。如果題目故意要求計(jì)算外推值那么在報(bào)告里要明確記錄這個(gè)現(xiàn)象用誤差對(duì)比說明外推不可靠這屬于正面的實(shí)驗(yàn)結(jié)果。如果在做擬合預(yù)測(cè)用polyfit擬合出的多項(xiàng)式做短距離預(yù)測(cè)是可行的做遠(yuǎn)距離預(yù)測(cè)同樣會(huì)失效。邊界意識(shí)是從數(shù)值分析實(shí)驗(yàn)里能帶走的最重要工程習(xí)慣。6. 用一份可視化腳本把實(shí)驗(yàn)結(jié)果整理成答辯圖表實(shí)驗(yàn)四答辯或提交報(bào)告時(shí)圖表質(zhì)量往往直接影響印象分。同樣是解決同一個(gè)問題一個(gè)用默認(rèn)繪圖風(fēng)格、圖例位置遮擋曲線、坐標(biāo)軸沒標(biāo)單位和另一個(gè)圖面干凈、曲線粗細(xì)分明、每個(gè)子圖都有標(biāo)題和注釋給閱卷人的感受完全不同。整理圖表這一步建議獨(dú)立寫一個(gè)腳本不要和算法代碼混在一起輸出為高分辨率圖片再插入報(bào)告而不是直接截圖。import matplotlib.pyplot as plt import numpy as np fig, axes plt.subplots(1, 3, figsize(14, 4)) # 第一個(gè)子圖插值結(jié)果對(duì)比 axes[0].plot(x_points, y_points, o, label數(shù)據(jù)點(diǎn)) axes[0].plot(x_dense, lagrange_interp(x_points, y_points, x_dense), label拉格朗日) coef newton_divided_differences(x_points, y_points) axes[0].plot(x_dense, newton_interp(coef, x_points, x_dense), --, label牛頓) axes[0].set_title(插值結(jié)果對(duì)比) axes[0].legend() # 第二個(gè)子圖擬合效果 axes[1].plot(x_data, y_data, o, label觀測(cè)數(shù)據(jù)) axes[1].plot(x_smooth, y_fit, -, label擬合曲線) axes[1].set_title(最小二乘擬合) # 第三個(gè)子圖殘差分布 axes[2].plot(x_data, residuals, o-, label殘差) axes[2].axhline(0, colorgray, linestyle--) axes[2].set_title(擬合殘差分布) axes[2].legend() fig.tight_layout() plt.savefig(experiment4_summary.png, dpi200, bbox_inchestight) plt.show()一處細(xì)節(jié)是plt.savefig里的參數(shù)dpi至少200才夠打印清晰bbox_inchestight可以把空白邊緣裁掉避免圖片四周留白過多。如果你是準(zhǔn)備提交電子版報(bào)告300dpi更好圖片文件大一點(diǎn)沒關(guān)系打印效果勻稱。另一個(gè)細(xì)節(jié)是每個(gè)子圖的標(biāo)題要能自解釋也就是說讀者不看正文只看圖就知道這張圖在表達(dá)什么這是實(shí)驗(yàn)報(bào)告圖表規(guī)范的基本要求多數(shù)同學(xué)沒意識(shí)到所以認(rèn)真做反而會(huì)留下印象。我自己的習(xí)慣是在畫圖腳本開頭放一段注釋記錄這份圖表對(duì)應(yīng)的題目要求、數(shù)據(jù)文件、生成日期方便復(fù)查時(shí)快速定位。數(shù)值分析實(shí)驗(yàn)四做完回頭看最大的收獲往往不是那幾個(gè)算法本身而是理解了“數(shù)據(jù)點(diǎn)給出的信息量是有限的任何試圖從有限數(shù)據(jù)外推的做法都必須說明假設(shè)和邊界”。這個(gè)意識(shí)在后面的課程和工程項(xiàng)目里反復(fù)被驗(yàn)證。希望這篇筆記能幫你把實(shí)驗(yàn)四從頭到尾走通拿到一個(gè)能進(jìn)報(bào)告、能講清楚、也能經(jīng)得起提問的結(jié)果。本文還有配套的精品資源點(diǎn)擊獲取