【正文】 本科畢業(yè)設(shè)計(jì)論文 33 參考文獻(xiàn) [1] William H. Press. NUMERICAL RECIPES, Third Edition, Cambridge University Press[J]. 2020 [2] Ardence Company, RTX SDK Release Notes, US. 2020. [3] GSL Reference Manual Edition . 2020 [4]飛行器系統(tǒng)仿真 [M] . 2020 [5] Kenh C. Johnson .Multidimensional Interpolation Methods . 2020 。特別感謝我的導(dǎo)師，感謝許老師在學(xué)習(xí)和研究上對(duì)我 的悉心指導(dǎo)和不斷鼓勵(lì)，這使我得到了許多的鍛煉和提高。但是這些都是比較特殊的情況，在一般的仿真應(yīng)用中，本文中算法的效率可以滿(mǎn)足需求。 當(dāng)然，本文也存在很多問(wèn)題，比如在 5 維以上的數(shù)據(jù)中插值時(shí)，除了線(xiàn)性插值外，其他算法的運(yùn)行速度將會(huì)急劇下降。插值類(lèi)提供了常見(jiàn)的幾種插值算法， 5 維以下時(shí)插值速度較快，接口設(shè)計(jì)簡(jiǎn)單， xml 格式插值表具有很強(qiáng)的可擴(kuò)展性。 向量和矩陣類(lèi)提供了很多的方法，可以非常靈活的使用向量和矩陣進(jìn)行數(shù)據(jù)處理。對(duì)相同的數(shù)據(jù)，多項(xiàng)式插值和樣條插值的速度遠(yuǎn)比線(xiàn)性插值慢。 線(xiàn)性 插值 拉格朗日 三次樣條 1 階 2 階 3 階 4 階 1 維 7874150 5774730 4122784 2769474 1877545 1465980 2 維 3663710 2021627 836020 412075 222817 269202 3 維 1969490 823054 209880 77002 47331 59627 4 維 1097230 357223 59643 24421 14818 18854 表 41 插值算法性能對(duì)比 根據(jù)以上測(cè)試數(shù)據(jù)可知，對(duì)網(wǎng)格數(shù)據(jù)使用線(xiàn)性插值速度非常快，可以達(dá)到每秒鐘百萬(wàn)次。 測(cè)試算例 2 此算例對(duì)插值類(lèi)進(jìn)行測(cè)試，對(duì)拉格朗日和三次樣條測(cè)試分別進(jìn)行統(tǒng)計(jì)，插值數(shù)據(jù)為網(wǎng)格數(shù) 據(jù)。由于是采用 C++編寫(xiě)，并且是單線(xiàn)程運(yùn)算， 20 秒的測(cè)試結(jié)果已經(jīng)比較滿(mǎn)意了。 本科畢業(yè)設(shè)計(jì)論文 30 double pute_time= end_time start_time。 res = ()。 Matrixdouble mtx(1000,1000), res。然后使用 Matrix 類(lèi)的 inv 方法求逆，并統(tǒng)計(jì)運(yùn)行時(shí)間。 測(cè)試算例 1 此算例使用矩陣類(lèi) Matrix 來(lái)進(jìn)行矩陣求逆。 RTX 提供的關(guān)機(jī)處理函數(shù)使得應(yīng)用程序可以在 windows 出現(xiàn)錯(cuò)誤時(shí)繼續(xù)運(yùn)行，某些函數(shù)可以在沒(méi)有操作系統(tǒng)的情況下依然正常執(zhí)行。 3. 并且 RTX 提供了附加的向?qū)Ш驼{(diào)試工具，可以無(wú)縫的嵌入到微軟的標(biāo)準(zhǔn)開(kāi)發(fā)工具中。 2. 具有工業(yè)標(biāo)準(zhǔn)和有效的開(kāi)發(fā)調(diào)試工具。在包括工業(yè)自動(dòng)化，電信，測(cè)量，醫(yī)藥，政府，軍事，模擬等各方面企業(yè)級(jí)應(yīng)用中，使基于 windows 的應(yīng)用程序擁有嵌入式程序所具備的硬實(shí)時(shí)特征。 }。 InterpMethod method。 //可否外插 , 否則出錯(cuò)提示 ,外插使用的是線(xiàn)性插值 vectorExtrapControlamp。 //應(yīng)變量 , dependent variable const vectordoubleamp。 //自變量 , independent variable const vector vectordouble amp。 關(guān)于外插的選擇，使用如下枚舉量： /*外插使用線(xiàn)性插值 */ enum ExtrapControl { /**不允許外插 */ NEITHER, /** 當(dāng)輸入數(shù)據(jù) x x_min 時(shí)允許外插 */ MINEX, /** 當(dāng)輸入數(shù)據(jù) x x_max 時(shí)允許外插 */ MAXEX, /** 當(dāng)輸入數(shù)據(jù) x x_min 或 x x_max 時(shí)允許外插 */ 本科畢業(yè)設(shè)計(jì)論文 28 BOTH }。接口不再羅列。39。39。 三次樣條插值 三次樣條插值的優(yōu)點(diǎn)是可以使插值曲線(xiàn)在插值區(qū)間內(nèi)部和端點(diǎn)出一階導(dǎo)數(shù)光滑，二階導(dǎo)數(shù)連續(xù)。 //連乘 double Add(size_t next)。 jlo)。 private: virtual double RawInterp(const vectordoubleamp。 interp_data)。如下 圖所示，其中 mulRes 為 n個(gè)連乘的結(jié)果。 拉格朗日插值 拉格朗日插值屬于多項(xiàng)式插值，其一維 n 點(diǎn)插值公式如下： ? ???? ??? nknkiikik i yxxxxyF it0 0 (32) n 維 n 點(diǎn)插值公式如下： ????????? ?????????? ????? ??? ???個(gè)個(gè)nnkii ikinnkii ikinkiikikinnknknknnnn nnnnn xxxxxxxxyxxxxy F i t????????????????????????????????????? ??????????????? 002010002222 2221111 1112211 (33) 上述 n 維插值公式翻譯成代碼的難點(diǎn)是連乘和連加的個(gè)數(shù)不確定，是個(gè)變量，用for 循環(huán)是很難完成的。 圖 35 一維線(xiàn)性插值 圖 36 一維線(xiàn)性 插值誤差 本科畢業(yè)設(shè)計(jì)論文 25 可以看到線(xiàn)性插值只需要兩個(gè)點(diǎn)即可進(jìn)行插值，計(jì)算量很小，速度是所有插值算法中最快的，可是誤差卻比較大，而且插值函數(shù)一階導(dǎo)數(shù)不連續(xù)，插值曲線(xiàn)不光滑。插值公式如下 101 0010 1 yxx xxyxx xxy F it ?????? (31) 比如對(duì)函數(shù)在區(qū)間 ]1,1[? 中選取等距的 6 個(gè)點(diǎn)進(jìn)行采樣，然后使用采樣點(diǎn)和采樣值進(jìn)行線(xiàn)性插值。 }。 idpVar, const vectorsize_tamp。 size)。 i)。 //在第 i 維自變量中尋找 x 所對(duì)應(yīng)的下標(biāo) , 使用加速器 ,局部二分法 size_t Hunt(const doubleamp。 x, const size_tamp。 itable)。 yy。 xx。 // 總區(qū)間分段 }ACCEL。 //是否使用加速器 size_t cache。比較好的情況是在 6 的附近，比如 8的話(huà)，只需要 2 步即可覆蓋，再使用二分法，一步就可確定位置。下圖是一個(gè)示例。還需要對(duì)搜索區(qū)間進(jìn)行分段。下面對(duì)相關(guān)性搜索進(jìn)行詳述。所以在如果數(shù)據(jù)量很多，每次都是從頭到尾的二分法搜索的話(huà)，會(huì)執(zhí)行很多重復(fù)搜索步驟，效率也不會(huì)讓人滿(mǎn)意。 經(jīng)驗(yàn)告訴我們，在很多情況下，相鄰幾次的插值調(diào)用都是在某一維的一個(gè)小區(qū)間內(nèi)進(jìn)行的。但是注意這是建立在有序表的基礎(chǔ)上。二分法的時(shí)間復(fù)雜度是 O( N2log )。 另一種簡(jiǎn)單的搜索算法是二分法，利用有序表的有序特征進(jìn)行對(duì)半搜索。最直接的方法是從第一個(gè)值開(kāi)始一個(gè)一個(gè)比較，直到滿(mǎn)足上述條件。注意我們對(duì)輸入的插值表是有要求的，必須是有序遞增或遞減的，而且在同一維中不能有相同值。所以這就是基類(lèi) CBaseInterp 要完成的工作了。之后就可以從這個(gè)基類(lèi)縱派生出想要的派生類(lèi)，并在派生類(lèi)中完成具體的插值算 本科畢業(yè)設(shè)計(jì)論文 22 法，從而得到不同的插值類(lèi)，并且共享基類(lèi)的接口。它們之間的關(guān)系如下面的結(jié)構(gòu)圖所示。插值基類(lèi)并不是給用戶(hù)使用的，而是給派生類(lèi)使用。插值法就是尋求近似函數(shù)的方法之一 [4]。雖然其函數(shù)關(guān)系 )(xfy? 在某個(gè)區(qū)間 ],[ ba 上式客觀(guān)存在的，但是卻不知道具體的解析表達(dá)式，只能通過(guò)觀(guān)察、測(cè)量或?qū)嶒?yàn)得到函數(shù)在區(qū)間 ],[ ba 上一些離散的函數(shù)值，導(dǎo)數(shù)值等，因此，希望對(duì)這樣的函數(shù)用一個(gè)比較簡(jiǎn)單的函數(shù)表達(dá)式來(lái)近似地給出整體上的描述。 本科畢業(yè)設(shè)計(jì)論文 20 。數(shù)學(xué)上被證明為最好的算法是那些基于 ranlux 的算法。 } printf(\n)。 // [0,1)區(qū)間的浮點(diǎn)數(shù) int int_rand = ()。 i 10。 // use default seed printf (generator is %s\n, ())。 } 使用示例 include include using namespace gslcpp。 現(xiàn)在一切都非常明了，構(gòu)造函數(shù)的具體實(shí)現(xiàn)也是非常簡(jiǎn)單 CRandom::CRandom(RNG_ORDER type /*= RNG_TAUS*/) { if (type () || type RNG_ORDER(0)) { gsl_error(type of generator not exists, __FILE__, __LINE__, GSL_ENOMEM)。 } 所以 type_stack 的初始化可寫(xiě)成如下形式 const CRandom::TypeArray CRandom::type_stack = CRandom::TypeInit()。 (gsl_rng_cmrg)。 (gsl_rng_borosh13)。 /*type_stack 為我們構(gòu)造的堆棧 , 由于是此類(lèi)對(duì)象所共享的變量， 所以設(shè)為靜態(tài) */ static const TypeArray type_stack。 下面考慮成員變量。 /*省略代碼 */ }。 unsigned long int Min()。 /*打印所有的算法及其對(duì)應(yīng)的枚舉量 */ void printAllNames()。 /*返回 [min,max]區(qū)間的一個(gè)整數(shù) , b_pos 為 true 時(shí) ,返回非零值 */ unsigned long int Rand(bool b_pos = false)。 virtual ~CRandom()。至于使用哪種算法可以在構(gòu)造函數(shù)里使用枚舉變量進(jìn)行選擇，同時(shí)提供一個(gè)改變算法的成員函數(shù)已實(shí)現(xiàn)運(yùn)行時(shí)改變算法。但是用戶(hù)使用的時(shí)候，不必關(guān)心狀態(tài)到底是怎么變的，這是算法內(nèi)部要考慮的事。需要注意的是，隨機(jī) 數(shù)生成器并不是正弦或余弦等函數(shù)，對(duì)隨機(jī)數(shù)生成器的連續(xù)調(diào)用會(huì)產(chǎn)生不同的值，這正是我們想要的特性。我要做的就是使用 GSL 提供的算法，給出隨機(jī)數(shù)生成器的 C++統(tǒng)一接口，并且在運(yùn)行時(shí)可以改變隨機(jī)數(shù)算法和種子。 對(duì)隨機(jī)數(shù)算法的測(cè)試通常是使用數(shù)值仿真來(lái)做統(tǒng)計(jì)測(cè)試，用隨機(jī)數(shù)來(lái)估計(jì)一些能用概率論求解出精確值的量，然后和精確值做比較即可測(cè)量該算法的隨機(jī)性。 隨機(jī)數(shù) 一個(gè)好的隨機(jī)數(shù)生成算法必須能夠在理論上和統(tǒng)計(jì)上同時(shí)滿(mǎn)足所需特性。 矩陣類(lèi)的構(gòu)架 由于矩陣和向量的需求分析有很多相似的地方，比如都要求有視圖類(lèi)，所以可以使用向量類(lèi)的構(gòu)架方法。 4． cols_，列數(shù)。 2． startCol_，起始列。 }。 size_t rows_。 class mslice { public: ………… . private: size_t startRow_。但是如何把它當(dāng)做二維數(shù)據(jù) 來(lái)讀取是個(gè)問(wèn)題。 矩陣 GSL 提供的矩陣類(lèi)型使用起來(lái)非常麻煩，由于是對(duì) c 語(yǔ)言中的數(shù)組進(jìn)行包裝，所以用戶(hù)時(shí)刻都要注意內(nèi)存的變化，花費(fèi)大量精力減少內(nèi)存泄露 ,并且每種數(shù)據(jù)類(lèi)型都有不同的函數(shù)版本，且不支持用戶(hù)自定義類(lèi)型，所以我決定使用 GSL的算法，重寫(xiě)構(gòu)架矩陣類(lèi)，并且給出 C++的接口，支持多種類(lèi)型，包括未知的用戶(hù)自定義類(lèi)型，盡量減少用戶(hù)對(duì)程序本身的維護(hù)工作量，使用戶(hù)能專(zhuān)注于算法本 本科畢業(yè)設(shè)計(jì)論文 14 身。 另外，由于使用標(biāo)準(zhǔn)庫(kù)的 valarray模板，所以在所有這些類(lèi)的實(shí)現(xiàn)中，甚至沒(méi)有主動(dòng)使用 C++標(biāo)準(zhǔn)內(nèi)存分配和刪除操作符 new， delete，一個(gè)也沒(méi)有出現(xiàn)。通常來(lái)說(shuō)，初始化列表必須是在編譯時(shí)可確定的，但用上面的方式，初始化列表里的 元素可以使變量，甚至是一個(gè)函數(shù)的返回值，只要能轉(zhuǎn)換為 Vector模板對(duì)應(yīng)的實(shí)例化類(lèi)型就可以了。 輸入數(shù)據(jù)的個(gè)數(shù)必須和對(duì)象的大小一致，否則會(huì)觸發(fā)庫(kù)中的出錯(cuò)處理函數(shù)。 v , , , 1, 。 為了解決這個(gè)問(wèn)題，我使用了 操作符，把我的類(lèi)對(duì)象當(dāng)做標(biāo)準(zhǔn)輸出。 // Error – invalid 不幸的是這不符合 C++的語(yǔ)法規(guī)則，變通一下，使用如下方式： Vectordouble v((double[]) { , , , 1, }, 5)。 但這種方法的不便之處是使用了局部變量。當(dāng)然也可以用下面這種