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