簡介北京理工大學(xué)光電學(xué)院《光電成像器件計算機(jī)輔助設(shè)計CAD》課程實驗報告核心內(nèi)容是基于電場高斯定律、環(huán)路定律與電勢差定理配合連續(xù)超張弛迭代算法與拉格朗日插值算法用C語言實現(xiàn)像管內(nèi)電勢計算的CAD程序設(shè)計。報告完整記錄了從網(wǎng)格劃分、有限差分方程建立到超張弛迭代收斂、邊界封閉處理的數(shù)值求解流程并給出了圓柱形像管內(nèi)任意點(diǎn)電位計算、對稱軸電位拉格朗日插值的實現(xiàn)方法。資源共1個PDF文檔壓縮包743KB適合光電工程、微光夜視技術(shù)方向的高年級本科生與研究生閱讀也可供從事電子光學(xué)系統(tǒng)仿真及數(shù)值計算工作的工程技術(shù)人員參考。已有171人學(xué)習(xí)下載。通過這份實驗報告讀者可掌握將電磁場理論模型轉(zhuǎn)化為可運(yùn)行程序的方法理解CAD軟件在像管電場分布計算與等勢線繪制中的具體應(yīng)用為像管設(shè)計、電子軌跡分析與像質(zhì)評估打下工程實踐基礎(chǔ)。1. 像管電場計算的一種落地解法用C語言把拉普拉斯方程變成可迭代的CAD工具拿到一份像管設(shè)計圖陰極、陽極、熒光屏電極形狀不規(guī)則手寫拉普拉斯方程去求解析解幾乎是不可能完成的任務(wù)。解析法只對少數(shù)對稱理想模型有效實驗測量法又受制于探頭尺寸和場擾動精度很難上去。這時候用C語言寫一個數(shù)值求解器在網(wǎng)格上把偏微分方程離散成差分方程通過連續(xù)超松弛迭代在毫秒級逼近真實電場分布再用拉格朗日插值求任意位置的電勢就是一條非常務(wù)實的工程路徑。這個項目正是一個完整的最小實現(xiàn)代碼量不算大卻把有限差分、連續(xù)超松弛迭代、等勢線掃描繪制、疊加定理驗證串成了一條完整的計算機(jī)輔助設(shè)計工具鏈。它不僅算出了像管內(nèi)部的電勢分布還把迭代因子、殘差、等勢線坐標(biāo)全部落盤到文本文件非常適合想搞懂“電場數(shù)值計算到底怎么落地成軟件”的開發(fā)者尤其是對電磁場數(shù)值解法剛上手的研究生或工程師。2. 理論基礎(chǔ)與整體架構(gòu)從三大靜電場定理到五段式模塊設(shè)計2.1 靜電場三大定理為什么是數(shù)值算法的“免死金牌”在著手寫代碼之前必須先把物理邊界條件想清楚。靜電場有三大基本定理這不僅是理論推導(dǎo)的基石更是你調(diào)試程序時判斷結(jié)果對錯的標(biāo)尺。第一個是唯一性定理。如果給定系統(tǒng)中所有電極的形狀、排列以及每個電極的電位那么由這些電極產(chǎn)生的靜電場將由拉普拉斯方程唯一確定。這句話的工程含義很直接只要你的網(wǎng)格足夠密、邊界條件給得正確不管你用迭代法還是直接法求解得到的電位分布應(yīng)該收斂到同一個值。如果你在程序中嘗試了多種初始猜測值最后迭代出來的結(jié)果不一樣那一定是你的邊界封閉處理或者差分格式寫錯了而不是算法發(fā)散的問題。第二個是相似性定理。若各電極的電位都增大為K倍空間各點(diǎn)的電位也增大為K倍等勢面形狀不變。這個定理在驗證程序時非常有用你把輸入文件里所有電極電壓都乘2算出來的每個網(wǎng)格點(diǎn)電壓也應(yīng)該乘2等勢線坐標(biāo)則應(yīng)該完全相同。如果做不到這一點(diǎn)說明你的程序中存在某些“硬編碼”的常數(shù)影響了結(jié)果的線性比例關(guān)系。第三個是多電極系統(tǒng)的電位疊加定理。當(dāng)各電極形狀、相對位置確定后總電場分布等于各電極單獨(dú)作用其余電極接地時電場的線性疊加。這個定理是程序中“疊加定理驗證模塊”的核心依據(jù)。實現(xiàn)上就是跑n遍單電極的迭代把結(jié)果加起來再跟一次跑完所有電極的完整結(jié)果做逐點(diǎn)對比。這既是對算法的驗證也是對邊界封閉處理正確性的交叉檢查。2.2 不等距五點(diǎn)差分格式從拉普拉斯方程到線性方程組在無空間電荷的真空中靜電場電位滿足拉普拉斯方程。解析求解時這是偏微分方程但數(shù)值求解時需要把它離散成差分方程。程序采用旋轉(zhuǎn)對稱場下的“十”字形不等距五點(diǎn)差分格式。假設(shè)網(wǎng)格點(diǎn)坐標(biāo)為(z, r)電位為φ0相鄰上下左右四個點(diǎn)的間距分別為h1、h2、h3、h4取絕對值則差分公式可以寫成φ0 (c1*φ1 c2*φ2 c3*φ3 c4*φ4) / c0其中系數(shù)c1到c4由相鄰點(diǎn)的間距決定c0是它們的加權(quán)和。對于軸上點(diǎn)r0由于徑向方向的間距存在0點(diǎn)需要用洛必達(dá)法則處理c3、c4的計算公式。這部分是程序里最容易寫錯的地方因為一旦網(wǎng)格劃分不均勻系數(shù)的計算就不能用簡單的“等間距差分模板”套必須根據(jù)實際網(wǎng)格坐標(biāo)動態(tài)計算各點(diǎn)的間距和權(quán)重。實現(xiàn)時我的習(xí)慣是單獨(dú)寫一個系數(shù)計算函數(shù)輸入是當(dāng)前網(wǎng)格點(diǎn)的四個鄰居坐標(biāo)輸出是c0到c4。這樣等間距網(wǎng)格和不等距網(wǎng)格可以共用同一套差分迭代邏輯后續(xù)如果要換網(wǎng)格劃分方式就不需要動主迭代循環(huán)了。2.3 程序模塊劃分輸入、初始化、迭代、繪圖、驗證五段式這個程序的整體架構(gòu)是典型的“串行流水線”風(fēng)格模塊邊界非常清晰非常像早期C語言工程的標(biāo)準(zhǔn)組織方式。全程序一共分成五個功能模塊輸入模塊負(fù)責(zé)讀入像管參數(shù)、初始化部分變量并原樣輸出參數(shù)到日志文件。初始化模塊負(fù)責(zé)劃分網(wǎng)格步長初始化電位矩陣用線性插值法封閉像管邊界。超松弛迭代模塊這是整個CAD程序的心臟負(fù)責(zé)對像管內(nèi)部電勢分布進(jìn)行數(shù)值計算。繪圖與輸出模塊進(jìn)行行列掃描得到等位線坐標(biāo)用EasyX插件繪制電極和等位線圖形再進(jìn)行拉格朗日插值繪制軸上電位變化曲線。疊加定理驗證模塊將每個電極單獨(dú)置為有效電位其余電極為0調(diào)用迭代模塊計算單位電位分布函數(shù)最后疊加并與完整計算結(jié)果對比。這種五段式架構(gòu)的優(yōu)勢在于“可替換性”。比如你不想用EasyX繪圖只需要改繪圖模塊其他模塊完全不用動你想把迭代算法換成共軛梯度法也只需要替換超松弛迭代模塊本身輸入輸出格式和行為完全兼容。3. 核心代碼拆解SOR迭代器與等勢線掃描的實現(xiàn)細(xì)節(jié)3.1 差分方程系數(shù)計算與邊界封閉的C語言實現(xiàn)這里給出系數(shù)計算的簡化版本。需要注意的是實際工程代碼通常會預(yù)先計算好所有網(wǎng)格點(diǎn)上的系數(shù)并存入數(shù)組中避免每次迭代都重復(fù)計算因為這是典型的“內(nèi)存換時間”優(yōu)化策略。// 計算不等距五點(diǎn)差分格式的系數(shù) // 輸入當(dāng)前網(wǎng)格點(diǎn)坐標(biāo) (z, r) 以及上下左右四個鄰居的坐標(biāo) // 輸出差分系數(shù) cc[0] 到 cc[4]其中 cc[0] 是中心點(diǎn)權(quán)重分母 void calc_coeff(double z, double r, double h1, double h2, double h3, double h4, double *cc) { double c1, c2, c3, c4; // 軸向差分系數(shù)直接按等距公式計算 c1 2.0 / (h1 * (h1 h2)); c2 2.0 / (h2 * (h1 h2)); // 徑向差分系數(shù)需要處理 r0 的奇異性 if (r 1e-9) { // 軸上點(diǎn)使用洛必達(dá)法則后的極限公式 c3 4.0 / (h3 * (h3 h4)); c4 4.0 / (h4 * (h3 h4)); } else { // 非軸上點(diǎn)常規(guī)公式 c3 (2.0 * r h3) / (h3 * r * (h3 h4)); c4 (2.0 * r - h4) / (h4 * r * (h3 h4)); } cc[0] c1 c2 c3 c4; cc[1] c1; cc[2] c2; cc[3] c3; cc[4] c4; }代碼里的關(guān)鍵點(diǎn)在于徑向系數(shù)的來源。在柱坐標(biāo)下拉普拉斯算子的徑向部分帶有1/r項直接做差分會在r0處產(chǎn)生奇點(diǎn)。程序的處理方式是對軸上點(diǎn)單獨(dú)推導(dǎo)極限公式這能保證對稱軸上的電位計算精度不會因為坐標(biāo)奇異性而明顯劣化。實際調(diào)試時如果發(fā)現(xiàn)軸上電位曲線出現(xiàn)毛刺或不光滑先檢查這一段的系數(shù)計算是否正確。邊界封閉的處理同樣重要。像管電極并不是完全封閉的電極之間存在空白區(qū)域這些區(qū)域的邊界電位是未知的。如果不做處理迭代計算就缺少邊界值無法收斂。程序采用最簡單的線性插值法在兩個相鄰電極的邊界之間按照距離比例給空白區(qū)域內(nèi)的邊界點(diǎn)賦予電位值。這樣做帶來的誤差在實際像管設(shè)計中是可以接受的因為遠(yuǎn)離電極的邊界區(qū)域?qū)S上電位的影響本身就很小。3.2 連續(xù)超松弛迭代核心循環(huán)從初始猜測到動態(tài)調(diào)整連續(xù)超松弛迭代的核心思想是每次遍歷網(wǎng)格時用當(dāng)前點(diǎn)的四個鄰居更新值來計算新的電位再引入一個松弛因子ω來加速收斂。當(dāng)1 ω 2時收斂過程會明顯加快但ω過大會導(dǎo)致震蕩甚至發(fā)散。程序中的迭代策略分為三個階段第一階段用ω1進(jìn)行純高斯-賽德爾迭代作為基準(zhǔn)啟動第二階段用固定的ω1.375迭代12次之后每12次迭代作為一個“輪次”在輪次結(jié)束時重新計算最優(yōu)松弛因子并更新全局的ω值。// 連續(xù)超松弛迭代核心循環(huán) // phi: 電位矩陣, R: 行數(shù), C: 列數(shù) // eps: 迭代控制精度, max_round: 最大輪次 void sor_solve(double **phi, int R, int C, double eps, int max_round) { double omega 1.0; double max_res 0.0, avg_res 0.0; for (int rnd 1; rnd max_round; rnd) { // 每輪迭代次數(shù)前兩輪特殊處理后續(xù)每輪12次 int iters_per_round (rnd 2) ? (rnd 1 ? 1 : 12) : 12; for (int k 0; k iters_per_round; k) { max_res 0.0; avg_res 0.0; int cnt 0; // 遍歷所有內(nèi)部網(wǎng)格點(diǎn) for (int i 1; i R - 1; i) { for (int j 1; j C - 1; j) { double old phi[i][j]; double new_val 0.0; // 這里需要根據(jù)網(wǎng)格點(diǎn)類型選取對應(yīng)的差分系數(shù) double c0, c1, c2, c3, c4; get_coeff_at(i, j, c0, c1, c2, c3, c4); new_val (c1 * phi[i-1][j] c2 * phi[i1][j] c3 * phi[i][j-1] c4 * phi[i][j1]) / c0; // SOR 更新公式 phi[i][j] (1.0 - omega) * old omega * new_val; double res fabs(phi[i][j] - old); max_res fmax(max_res, res); avg_res res; cnt; } } avg_res / cnt; // 輸出本輪迭代的殘差 printf(第 %d 輪第 %d 次迭代迭代因子: %f 平均殘差%f , 最大殘差%f\n, rnd, k 1, omega, avg_res, max_res); // 滿足精度要求則提前退出 if (max_res eps) { printf(迭代收斂\n); return; } } // 每輪迭代結(jié)束后更新松弛因子 // 實際代碼中omega的計算會用到譜半徑估計這里簡化為遞增策略 if (rnd 1) omega 1.375; else if (rnd 2) omega 1.75; else omega fmin(1.95, omega 0.05); } }這段代碼最需要注意的地方有兩個。第一new_val的更新必須使用“最新值”也就是當(dāng)前輪次已經(jīng)計算過的節(jié)點(diǎn)值這是高斯-賽德爾迭代的基本要求。如果你用上一輪次的舊值做更新那就變成了雅可比迭代收斂速度會大幅下降。第二omega的更新策略需要根據(jù)實際殘差變化做調(diào)整。在原始實驗的迭代日志中可以看到omega從1.0逐步爬到1.85左右整個過程非常平滑。如果你在調(diào)試自己代碼時發(fā)現(xiàn)omega調(diào)整后殘差反而變大說明當(dāng)前問題域需要用更保守的更新策略。3.3 行列掃描繪制等勢線把電位矩陣變成坐標(biāo)點(diǎn)集迭代完成之后你手里有一整張網(wǎng)格的電位值。要畫等勢線需要把“哪幾個點(diǎn)電位相等”這個問題轉(zhuǎn)化為算法可執(zhí)行的操作。程序采用的方法非常經(jīng)典行掃描加列掃描然后對相鄰網(wǎng)格點(diǎn)做線性插值。// 掃描并記錄所有等勢線點(diǎn) // phi: 電位矩陣, V_level: 當(dāng)前要搜索的電位值 // points: 輸出點(diǎn)集, max_points: 點(diǎn)集容量 int scan_equipotential(double **phi, int R, int C, double V_level, Point *points, int max_points) { int count 0; // 行掃描沿著每一行檢查水平方向穿越 for (int i 1; i R - 1; i) { for (int j 1; j C - 1; j) { if (phi[i][j] V_level phi[i][j1] V_level) { // 在線性插值求精確z坐標(biāo) double t (V_level - phi[i][j]) / (phi[i][j1] - phi[i][j]); double z j t * 1.0; // 假設(shè)網(wǎng)格間距為1 double r i * 1.0; if (count max_points) { points[count].z z; points[count].r r; count; } } } } // 列掃描沿著每一列檢查垂直方向穿越 for (int j 1; j C - 1; j) { for (int i 1; i R - 1; i) { if (phi[i][j] V_level phi[i1][j] V_level) { // 在線性插值求精確r坐標(biāo) double t (V_level - phi[i][j]) / (phi[i1][j] - phi[i][j]); double r i t * 1.0; double z j * 1.0; if (count max_points) { points[count].z z; points[count].r r; count; } } } } // 對點(diǎn)集按z坐標(biāo)排序 qsort(points, count, sizeof(Point), compare_z); return count; }行掃描和列掃描都會產(chǎn)生等勢線點(diǎn)但單獨(dú)用其中一種會漏掉那些與網(wǎng)格線平行的等勢線段。行列結(jié)合能最大程度保證等勢線的連續(xù)性。掃描結(jié)束后得到的點(diǎn)集并不是有序的繪圖路徑所以程序還需要對點(diǎn)集按照z坐標(biāo)進(jìn)行排序然后按順序連接。這個排序步驟在原始報告中描述得很清楚先按照z坐標(biāo)增大順序排序連接一次再按照z坐標(biāo)減小順序排序連接一次主要目的是處理等勢線在像管底部發(fā)生回折的情況。4. 把程序跑起來dat輸入文件與EasyX繪圖環(huán)境配置4.1 輸入文件參數(shù)詳解從電極個數(shù)到迭代精度輸入文件是后綴為.dat的文本文件每一行對應(yīng)一類參數(shù)。文件的組織方式幾乎是“無格式”的純序列文本參數(shù)的順序完全固定任何一行讀錯了后面的數(shù)據(jù)全部錯位。參數(shù)含義示例值n電極總數(shù)不包括電位為0V的陰極9δ電極厚度mm占1個步長0.4Δzi相鄰電極之間的距離mm一共n-1個5, 7.1, 7.1...Ni相鄰電極之間要劃分的步長數(shù)10, 12, 12...Vi電極電位V21, 38, 42...Δr1電極內(nèi)孔徑半徑mm16M1Δr1范圍內(nèi)等步長劃分的網(wǎng)格數(shù)25Δr2從電極內(nèi)孔邊沿到封閉邊界處的徑向距離mm10M2Δr2范圍內(nèi)等步長劃分的網(wǎng)格數(shù)16ε迭代控制精度0.00002NST輸出打印空間電位時網(wǎng)格點(diǎn)間隔數(shù)1或2INS軸上電位等距插值時的步長數(shù)80ΔV掃描搜索等勢線電位間隔8m特定等勢線掃描電位個數(shù)任意從輸入格式可以看到程序?qū)W(wǎng)格劃分的約束是M1 M2 50所有電極厚度相同所有電極內(nèi)孔徑半徑相同。這些約束在原始代碼里被寫成了硬編碼的數(shù)組大小如果直接跑實驗數(shù)據(jù)沒問題但想改成像管結(jié)構(gòu)就需要動全局常量。4.2 開發(fā)環(huán)境配置與編譯運(yùn)行步驟程序使用C語言編寫繪圖依賴EasyX插件所以運(yùn)行的推薦環(huán)境是Windows系統(tǒng)加Visual Studio。具體配置步驟為第一步下載并安裝Visual Studio社區(qū)版第二步從官網(wǎng)下載EasyX庫運(yùn)行安裝程序它會自動識別當(dāng)前環(huán)境中的Visual Studio版本并完成集成第三步新建一個空項目把源碼文件加入項目在項目屬性中確認(rèn)“字符集”設(shè)置為“使用多字節(jié)字符集”因為早期代碼習(xí)慣用char*處理文件名如果是Unicode字符集可能會報類型不匹配的編譯錯誤第四步把輸入文件放到與main.c相同的目錄下直接編譯運(yùn)行。運(yùn)行時程序會在控制臺窗口打印迭代過程包括每一輪次、每一次迭代的迭代因子、平均殘差和最大殘差就像附錄一展示的那樣。全部迭代完成后會彈出EasyX繪圖窗口顯示等勢線和軸上電位曲線同時會在輸入文件同目錄下生成一個.txt文件里面包含了完整的迭代結(jié)果和等勢線坐標(biāo)。由于繪圖窗口是圖形界面程序退出后窗口會消失如果看不到繪圖結(jié)果需要在main函數(shù)末尾加上system(pause);或者getchar();來保留窗口。4.3 避坑指南迭代不收斂、邊界封閉失效與繪圖黑屏現(xiàn)象一迭代到第3輪時殘差非但不下降反而突然增大。原因松弛因子ω更新過快跳過了最優(yōu)點(diǎn)導(dǎo)致迭代過程進(jìn)入震蕩區(qū)間。原始日志中ω從1.375跳到1.755時最大殘差從0.103反彈到0.280這是典型的過松弛震蕩。解決把ω的更新步長調(diào)小每次只增加0.02到0.03或者引入自適應(yīng)機(jī)制當(dāng)檢測到最大殘差連續(xù)兩次上升時自動回退ω值?，F(xiàn)象二繪制的等勢線在電極之間的缺口處斷開或者明顯扭曲。原因邊界封閉處理不正確。如果兩個相鄰電極的電位差很大線性插值得到的邊界點(diǎn)電位可能不符合實際的物理衰減規(guī)律導(dǎo)致靠近邊界的網(wǎng)格點(diǎn)電位被錯誤牽引。解決檢查邊界封閉的賦值順序確保左右邊界、上下邊界各自獨(dú)立處理不要讓插值公式跨邊界引用。另外可以嘗試在電極之間增加少量過渡網(wǎng)格點(diǎn)讓插值更平緩?，F(xiàn)象三EasyX繪圖窗口一閃而過根本來不及看結(jié)果。原因程序運(yùn)行完main函數(shù)后進(jìn)程退出圖形窗口被系統(tǒng)自動回收。解決在main函數(shù)的繪圖代碼后增加getchar();或system(pause);讓程序掛起等待鍵盤輸入。如果是在Visual Studio里按F5運(yùn)行也可以改用CtrlF5開始執(zhí)行但不調(diào)試這樣控制臺窗口結(jié)束后不會自動關(guān)閉。現(xiàn)象四修改網(wǎng)格數(shù)量M1或M2后程序直接崩潰或拋出越界錯誤。原因原始代碼的數(shù)組大小是按50x113硬編碼的M1M2超過50時就會越界寫入內(nèi)存。解決建議把這些數(shù)組定義改成動態(tài)分配使用malloc和free或者直接定義足夠大的常量比如200并加上輸入?yún)?shù)校驗在讀取M1M2時如果超過上限就報錯退出而不是等到運(yùn)行到一半才崩潰。5. 數(shù)據(jù)處理與誤差分析殘差日志、疊加定理驗證與實際工程意義5.1 迭代殘差日志解讀平均殘差與最大殘差的玄學(xué)迭代收斂的過程不是一條光滑下降的曲線而是帶著“毛刺”的震蕩衰減。查看附錄一里的殘差日志你會發(fā)現(xiàn)第3輪第1次迭代的平均殘差是0.056243比第2輪最后一次迭代的0.021593還要大。這正是調(diào)整松弛因子后的“回彈”現(xiàn)象。松弛因子的調(diào)整會讓算法暫時偏離原來的路徑但隨后會以更快的速度收斂。解讀殘差日志時重點(diǎn)關(guān)注兩點(diǎn)。第一是“最大殘差”它代表當(dāng)前網(wǎng)格中最不收斂的那個點(diǎn)的誤差大小。如果最大殘差一直不降說明某個區(qū)域可能被錯誤的邊界條件污染了比如邊界封閉線上有某個點(diǎn)的電位值沒被正確初始化會持續(xù)向周圍網(wǎng)格“發(fā)射”錯誤信號。第二是“平均殘差”與“最大殘差”的比值。正常情況下這個比值應(yīng)該在0.2到0.5之間。如果平均殘差很小但最大殘差很大說明整個場分布已經(jīng)基本穩(wěn)定但有幾個孤立點(diǎn)存在異常這種時候通常要回溯檢查那幾個點(diǎn)的坐標(biāo)位置看看是否落在了電極邊緣或邊界線上。5.2 疊加定理驗證模塊把n次SOR結(jié)果疊起來疊加定理的驗證邏輯非常樸素但工程價值很高。具體實現(xiàn)是先將所有電極電位清零只保留第1個電極的電位運(yùn)行一次完整的連續(xù)超松弛迭代模塊把結(jié)果存入一個獨(dú)立的疊加數(shù)組然后清零所有電極只保留第2個電極的電位再運(yùn)行一次迭代把結(jié)果累加到剛才的數(shù)組里重復(fù)這個過程直到所有電極都單獨(dú)處理完成。// 疊加定理驗證模塊的核心偽代碼 // base_phi: 正常計算得到的總電場電位 // super_phi: 疊加定理計算得到的電位 // single_phi: 單電極電場電位數(shù)組 void verify_superposition(double **base_phi, double **single_phi, double **super_phi, int n_electrode, int R, int C) { // 初始化疊加數(shù)組為零 for (int i 0; i R; i) for (int j 0; j C; j) super_phi[i][j] 0.0; // 循環(huán)處理每個電極 for (int k 1; k n_electrode; k) { // 設(shè)置電極電位只保留第k個電極的電壓其余電極置零 set_electrode_potential(k); // 調(diào)用SOR迭代模塊 sor_solve(single_phi, R, C, EPS, MAX_ROUND); // 累加到疊加數(shù)組 for (int i 0; i R; i) for (int j 0; j C; j) super_phi[i][j] single_phi[i][j] * V[k]; } // 對比疊加結(jié)果與直接計算結(jié)果 double max_diff 0.0, avg_diff 0.0, rms_diff 0.0; for (int i 0; i R; i) { for (int j 0; j C; j) { double diff fabs(base_phi[i][j] - super_phi[i][j]); max_diff fmax(max_diff, diff); avg_diff diff; rms_diff diff * diff; } } avg_diff / (R * C); rms_diff sqrt(rms_diff / (R * C)); printf(最大絕對殘差: %f, 平均殘差: %f, 均方根殘差: %f\n, max_diff, avg_diff, rms_diff); }這段代碼里有一個隱含的“坑”set_electrode_potential(k)函數(shù)需要重新初始化電位矩陣把非當(dāng)前電極的電壓全部清零。如果只是把電極電位值改掉而沒有清零邊界條件之前輪次的殘余電位會污染下一次迭代的初值導(dǎo)致疊加結(jié)果完全錯誤。我一般會在set_electrode_potential里顯式執(zhí)行一次全矩陣初始化把邊界條件重新賦值一遍確保每次迭代都在同一套邊界設(shè)置下進(jìn)行。疊加定理驗證最大的意義在于它能以極低成本檢查整個程序的邊界處理和迭代邏輯是否自洽。如果疊加驗證的最大絕對殘差超過預(yù)設(shè)的迭代精度一個數(shù)量級那幾乎可以肯定邊界封閉或者迭代更新順序里有隱藏的bug這種bug在單一組的迭代結(jié)果里很難被肉眼發(fā)現(xiàn)。5.3 誤差來源拆解截斷誤差、精度誤差和網(wǎng)格離散誤差任何一種數(shù)值計算方法的誤差都不能被完全消除只能被控制在可接受范圍內(nèi)。這個程序里的誤差來源主要有三類。第一類是計算機(jī)本身的截斷誤差由于CPU位數(shù)有限浮點(diǎn)數(shù)在存儲時只能保留有限的有效位數(shù)每一次加減乘除都會丟失最末尾幾個bit的精度。這類誤差會隨著迭代次數(shù)累積但累積速度很慢在精度要求為0.00002V的數(shù)量級下完全可以忽略。第二類是給定精度帶來的截斷誤差迭代結(jié)束條件是前后兩次殘差小于預(yù)設(shè)的ε這意味著所有小于ε的電位波動都被視為“已收斂”而忽略即使它可能是一個真實存在的微弱物理信號。第三類是網(wǎng)格離散誤差這是數(shù)值方法的核心誤差來源。網(wǎng)格永遠(yuǎn)是離散的電場在網(wǎng)格點(diǎn)之間的連續(xù)變化被線性近似替代網(wǎng)格越粗誤差越大。程序用行列掃描加線性插值來畫等勢線本身就是對離散誤差的一種補(bǔ)救但插值階數(shù)只有一階在曲率大的地方比如電極附近仍然會有明顯的位置偏移。6. 二次開發(fā)思路從圓柱形像管到任意旋轉(zhuǎn)對稱結(jié)構(gòu)的改造拿到這個程序后大部分人第一反應(yīng)就是“怎么改代碼能自定義電極形狀”。原始程序假設(shè)了圓柱形像管所有電極都是完整的環(huán)形電極內(nèi)孔徑半徑相同厚度相同。如果要做更復(fù)雜的結(jié)構(gòu)和參數(shù)可以從以下三個方向入手。第一個方向是網(wǎng)格分布優(yōu)化。把軸向網(wǎng)格改成非均勻分布在電極邊緣和陰極附近加密網(wǎng)格在其他區(qū)域適當(dāng)放寬網(wǎng)格間距。這需要修改初始化模塊中的步長數(shù)組計算邏輯并且讓差分系數(shù)計算函數(shù)支持任意疏密分布的網(wǎng)格。你可以對比一下在相同最大網(wǎng)格數(shù)的條件下均勻網(wǎng)格和非均勻網(wǎng)格算出來的軸上電位曲線誰更平滑。第二個方向是算法升級。把連續(xù)超松弛迭代換成多重網(wǎng)格法。原始報告參考文獻(xiàn)里也提到了多重網(wǎng)格法的應(yīng)用它在三維像管電場計算中能明顯加速收斂。實現(xiàn)思路是先在最細(xì)的網(wǎng)格上迭代幾次然后限制到粗網(wǎng)格上繼續(xù)迭代得到修正項后再延拓回細(xì)網(wǎng)格。這個過程循環(huán)往復(fù)能大幅度消除低頻誤差分量。相比單純調(diào)ω多重網(wǎng)格的收斂速度提升是指數(shù)級的。第三個方向是插值算法替換。程序目前用拉格朗日三項二次插值來計算對稱軸上任意點(diǎn)的電位優(yōu)點(diǎn)是代碼簡潔、容易實現(xiàn)缺點(diǎn)是精度受限。你可以把它替換成三次樣條插值在已知網(wǎng)格點(diǎn)電位和電極邊界值的前提下樣條插值能保證二階導(dǎo)數(shù)連續(xù)軸上電位曲線會更光滑計算電子軌跡時得到的偏導(dǎo)數(shù)也更可信。從投入產(chǎn)出比來看我建議你優(yōu)先做第一項改造即非均勻網(wǎng)格因為它的工程量最小但對計算精度和速度的提升立竿見影。網(wǎng)格分布優(yōu)化的代碼改造點(diǎn)集中在初始化模塊里只要把等間距的循環(huán)改為按比例函數(shù)分布即可。從那以后我每次把物理模型轉(zhuǎn)成代碼時都會強(qiáng)制走一遍“網(wǎng)格獨(dú)立性驗證”——用稀疏網(wǎng)格和加密網(wǎng)格各算一遍如果結(jié)果差異在可接受范圍內(nèi)才敢把程序交給別人用。希望幫到你。本文還有配套的精品資源點(diǎn)擊獲取