數(shù)值計算課程設(shè)計選編

1、數(shù)值計算課程設(shè)計說明書題目： 典型數(shù)值算法的C+語言程序設(shè)計 院 （系）： 理學(xué)院 專業(yè)班級： 信息131 學(xué) 號： 201312030120 學(xué)生姓名： 郝騰宇 指導(dǎo)教師： 劉海峰 2015年 7月 10 日陜 西 科 技 大 學(xué)數(shù)值計算課程設(shè)計任務(wù)書理學(xué)院 信息與計算科學(xué) 專業(yè) 信息131 班級 學(xué)生： 郝騰宇 題目：典型數(shù)值算法的C+語言程序設(shè)計 課程設(shè)計從 2015 年 5 月 20 日起到 2015 年 7月 10 日1、課程設(shè)計的內(nèi)容和要求（包括原始數(shù)據(jù)、技術(shù)要求、工作要求等）：每人需作10個算法的程序、必做6題、自選4題。對每個算法要求用C+語言進行編程。必選題：1、高斯列主元法

2、解線性方程組2、牛頓法解非線性方程組3、經(jīng)典四階龍格庫塔法解一階微分方程組4、三次樣條插值算法（壓緊樣條）用C+語言進行編程計算 依據(jù)計算結(jié)果，用Matlab畫圖并觀察三次樣條插值效果。5、龍貝格求積分算法6、M次多項式曲線擬合，據(jù)計算結(jié)果，用Matlab畫圖并觀察擬合效果。自選題：自選4道其他數(shù)值算法題目.每道題目重選次數(shù)不得超過5次.2、對課程設(shè)計成果的要求包括圖表、實物等硬件要求：2.1 提交課程設(shè)計報告按照算法要求，應(yīng)用C+語言設(shè)計和開發(fā)算法程序，提交由：每個算法的原理與公式說明；與算法相應(yīng)的程序設(shè)計說明（程序中的主要變量語義說明，變量的數(shù)據(jù)類型，數(shù)據(jù)在內(nèi)存中組織和存儲結(jié)構(gòu)說明，各函數(shù)

3、模塊的主要流程圖，函數(shù)功能說明，函數(shù)的形參說明，函數(shù)的調(diào)用方法說明）；程序調(diào)試與實例運行記錄 (包括程序調(diào)試和修改記錄、測試結(jié)論、運行結(jié)果記錄)，每個算法的源程序代碼編入附錄構(gòu)成的課程設(shè)計報告。2.2 課程設(shè)計報告版式要求目錄的要求：居中打印目錄二字，（四號黑體，段后1行），字間空一字符；章、節(jié)、小節(jié)及其開始頁碼（字體均為小四號宋體）。節(jié)向右縮進兩個字符（漢字），小節(jié)及以后標題均向右縮進四個字符（漢字）。目錄中應(yīng)包含正文部分每個算法章節(jié)標題、設(shè)計體總結(jié)、無序號的“參考文獻資料”，目錄的最后一項是“附錄”正文的要求：算法說明論述清楚，公式符號撰寫規(guī)范，流程圖圖符規(guī)范， 計算正確，文字簡練通順，插

4、圖簡潔規(guī)范，書寫整潔。文中圖、表按制圖要求繪制，程序調(diào)試和運行情況記錄詳實。打印版面要求：A4紙，頁邊距：上2cm，下2cmcm、右2cm；字體：正文宋體、小四號；行距：固定值20；頁眉1.5cm ，頁腳1.75cm；頁碼位于頁腳居中打?。黄鏀?shù)頁頁眉“數(shù)值計算課程設(shè)計”，偶數(shù)頁頁眉“具體算法名”，頁眉宋體小5號；段落及層次要求：每節(jié)標題以四號黑體左起打?。ǘ吻岸魏蟾?.5行），節(jié)下為小節(jié)，以小四號黑體左起打印（段前段后各0.5行）。換行后以小四號宋體打印正文。章、節(jié)、小節(jié)編號分別以1、1.1、1.1.1格式依次標出，空一字符后接各部分的標題。每一章的標題都應(yīng)出現(xiàn)在本章首頁的第一行上。當課程設(shè)計

5、報告結(jié)構(gòu)復(fù)雜，小節(jié)以下的標題，左起頂格書寫，編號依次用（1）、（2）或1）、2）順序表示。字體為小四號宋體。對條文內(nèi)容采用分行并敘時，其編號用（a）、（b）或a）、b）順序表示，如果編號及其后內(nèi)容新起一個段落，則編號前空兩個中文字符。曲線圖表要求：所有曲線、圖表、線路圖、流程圖、程序框圖、示意圖等不準徒手畫，必須按國家規(guī)定標準或工程要求繪制（應(yīng)采用計算機繪圖）。課程設(shè)計說明書（報告）中圖表、公式要求如下：（a）圖：圖的名稱采用中文，中文字體為五號宋體，圖號圖名在圖片下面。引用圖應(yīng)在圖題右上角標出文獻來源。圖號以章為單位順序編號。格式為：圖1-1，空一字符后接圖名，比如第1章第5個圖是關(guān)于高斯列

6、主元法解方程組算法圖，圖的下方的圖號圖名應(yīng)為：圖1-5 高斯列主元法解方程組算法圖。（b）表格：表的名稱及表內(nèi)文字采用中文，中文字體為五號宋體，表號表名在表格上面。表號以章為單位順序編號，表內(nèi)必須按規(guī)定的符號標注單位。格式為：表1-1，空一字符后接表格名稱。比如第4章第1個表是關(guān)于三次樣條插值的插值點列表，表的上方表號表名則應(yīng)為：表4-1 已知插值點的列表。（c）公式：公式書寫應(yīng)在文中另起一行，居中排列。公式序號按章順序編號。字體為五號宋體，序號靠頁面右側(cè)。比如第3章第1個公式其編號則應(yīng)為：（3-1）。設(shè)計體會及今后的改進意見：設(shè)計總結(jié)要寫出算法理解，編程經(jīng)驗等技術(shù)性、學(xué)術(shù)性總結(jié)；體會要簡潔、

7、真實、深刻，切忌空話、大話，客套話和矯揉造作之詞。改進意見要合理、中肯。參考文獻的要求：另起一頁，居中打印參考文獻四字（四號黑體，段前段后1行），字間空一字符；另起一行，按論文中參考文獻出現(xiàn)的先后順序用阿拉伯數(shù)字連續(xù)編號（參考文獻編號應(yīng)在正文中標注出）；參考文獻中每條項目應(yīng)齊全（字體均為小四號宋體）。（格式：編號作者.論文或著作名稱.期刊名或出版社.出版時間）。（期刊應(yīng)注明第幾期、起止頁數(shù)（包括論著）。參考文獻中條目要符合科技文獻引用文獻條目書寫的國家標準規(guī)范。2.3 設(shè)計報告裝訂順序與規(guī)范封面數(shù)值計算課程設(shè)計任務(wù)書目錄數(shù)值計算課程設(shè)計報告正文設(shè)計體會及今后的改進意見參考文獻(無需加目錄序號)

8、附錄（無需加目錄序號)左邊緣裝訂3、課程設(shè)計工作進度計劃：時間設(shè)計任務(wù)及要求第18周編寫和調(diào)試程序并按要求撰寫設(shè)計報告 目 錄1 高斯列主元法解線性方程組11.1 高斯列主元法解線性方程組算法思想11.2 高斯列主元法解線性方程組流程圖11.3 高斯列主元法解線性方程組運行結(jié)果22 牛頓法解非線性方程組32.1 牛頓法解非線性方程組算法思想32.2 牛頓法解非線性方程組流程圖42.3 牛頓法解非線性方程組運行結(jié)果52.4 牛頓法解非線性方程組手工計算63 經(jīng)典四階龍格庫塔法解一階微分方程組73.1 經(jīng)典四階龍格庫塔法解一階微分方程算法思想73.2 經(jīng)典四階龍格庫塔法解一階微分方程流程圖83.3

9、 經(jīng)典四階龍格庫塔法解一階微分方程運行結(jié)果94 三次樣條插值算法（壓緊樣條）104.1 三次樣條插值算法說明104.2 三次樣條插值算法流程圖 114.3 三次樣條插值運行結(jié)果124.4 Matlab畫圖并觀察三次樣條插值效果135 龍貝格求積分算法14龍貝格算法說明14龍貝格算法流程圖15龍貝格算法運行結(jié)果166 M次多項式擬合176.1 M次多項式擬合算法思想176.2 M次多項式擬合流程圖176.3 程序運行結(jié)果186.4 MATLAB擬合結(jié)果187 雅克比迭代法197.1 雅克比迭代法算法說明197.2 雅克比迭代法流程圖197.3 雅克比迭代法程序運行結(jié)果208 高斯-賽德爾迭代法2

10、18.1 算法說明218.2 高斯-賽德爾迭代法流程圖228.3 高斯-賽德爾迭代法流運行結(jié)果239 復(fù)化辛普森公式求積分249.1 復(fù)化辛普森公式求積分算法說明249.2 復(fù)化辛普森公式求積分流程圖259.3 復(fù)化辛普森公式運行結(jié)果2610 最小二乘線性擬合2710.1 最小二乘線性擬合算法說明2710.2 最小二乘法進行線性擬合流程圖2710.3 最小二乘法進行線性擬合運行結(jié)果2811 設(shè)計體會及改進意見29參 考 文 獻30附 錄31附錄一 高斯列主元法程序代碼31附錄二 牛頓法解非線性方程組的程序代碼35附錄三 經(jīng)典四階龍格庫塔法39附錄四 三次樣條插值算法（壓緊樣條）41附錄五 龍貝

11、格算法程序代碼44附錄六 M次多項式曲線擬合程序代碼46附錄七 雅克比迭代程序代碼50附錄八 高斯賽德爾迭代程序代碼52附錄九 復(fù)化辛普森求積公式程序代碼55附錄十 最小二乘法線性擬合程序代碼571 高斯列主元法解線性方程組問題提出：求解線性方程組1.1 高斯列主元法解線性方程組算法思想1） 將線性方程組寫成增廣矩陣的形式；2） 對增廣矩陣進行行變換，對元素，在第i列中及以下的元素選取絕對值最大的元素，將所有元素中最大的所在的行與第i行交換，然后采用高斯消元法使得新得到的第i行以下的元素均為零。一直重復(fù)上述過程直到。從而得到上三角矩陣。3）對上三角矩陣進行回代求解，即可以得到方程組的解。1.2

12、 高斯列主元法解線性方程組流程圖輸入增廣矩陣Start 輸入未知量的個數(shù)n高斯列主元法變換為上三角矩陣矩陣是否退化YN回代求解End 圖1-1 高斯列主元法流程圖1.3 高斯列主元法解線性方程組運行結(jié)果圖1-2 高斯列主元法程序運行界面圖2 牛頓法解非線性方程組問題提出：求解非線性方程組設(shè)初值，用牛頓法計算。2.1 牛頓法解非線性方程組算法思想前提：設(shè)已知1）計算函數(shù)（2-1）2）計算雅可比矩陣（2-2）(2-3)3）求線性方程組的解。(2-4)4）計算下一點重復(fù)上述過程，直到達到要求的精度，停止計算。2.2 牛頓法解非線性方程組流程圖Start 輸入精度要求和最大迭代次數(shù)數(shù) 計算計算雅可比矩

13、陣 =- 或者精度不夠Y EndN圖2-1 牛頓法解非線性方程組流程圖2.3 牛頓法解非線性方程組運行結(jié)果圖2-2 牛頓法解非線性方程組程序運行圖2.4 牛頓法解非線性方程組手工計算求解非線性方程組設(shè)初值，用牛頓法計算。函數(shù)向量和雅克比矩陣為(2-5)在點（2.00,0.25）處的函數(shù)值為(2-6)由牛頓法可得出下面的線性方程組 解得(2-7)迭代得(2-8)手工計算與程序運行結(jié)果一致。3 經(jīng)典四階龍格庫塔法解一階微分方程組問題提出：求解線性微分方程組 3.1 經(jīng)典四階龍格庫塔法解一階微分方程算法思想龍格-庫塔法的基本思想是：用在幾個不同點的數(shù)值加權(quán)平均來代替的值，而使截斷誤差的階數(shù)盡可能高。

14、我們用四個不同點上的函數(shù)值的線性組合，將精度提高到四階就可以得到四階龍格-庫塔公式。常用的標準四階龍格-庫塔公式為： (3-1)其中 (3-2) 經(jīng)典四階龍格庫塔法解一階微分方程流程圖 Start Y N輸出結(jié)果End圖3-1 四階龍格庫塔法程序流程圖3.3 經(jīng)典四階龍格庫塔法解一階微分方程運行結(jié)果圖3-2 四階龍格庫塔法程序運行界面4 三次樣條插值算法（壓緊樣條）4.1 三次樣條插值算法說明設(shè)已知及，插值函數(shù)在每個小區(qū)間上是不超過3次的多項式且具有2階連續(xù)導(dǎo)數(shù)，則稱為三階樣條插值。具體的說，三階樣條插值是滿足下列條件的分段三次多項式：（1） 插值條件：(4-1)； (2)連接條件：(4-2)

15、壓緊樣條公式(4-3)三階樣條插值的邊界條件為（1）已知兩端的一階導(dǎo)數(shù)；（2）兩端的二階導(dǎo)數(shù)已知；（3）默認:自然邊界條件。4.2 三次樣條插值算法流程圖 Start輸入所有點坐標及邊界條件計算三次樣條插值系數(shù)輸出結(jié)果End圖4-1 三次樣條插值程序流程圖4.3 三次樣條插值運行結(jié)果圖4-2 三次樣條差值運行界面圖4.4 Matlab畫圖并觀察三次樣條插值效果圖4-3 三次樣條插值Matlab畫圖5 龍貝格求積分算法問題提出:用龍貝格積分法求解以下積分：(5-1)步1 輸入a,b及精度；步2 步3 對分區(qū)間a,b，并計算(5-2)步4 若不滿足終止條件，作循環(huán)：(5-3)(5-5)(5-4)

16、計算 對計算終止條件一般取為手工計算過程如下：開始 輸入下限a、上限b、精度。 ,，計算 得到積分值結(jié)束圖5-1 龍貝格積分程序流程圖圖5-2 龍貝格積分程序運行圖6 M次多項式擬合問題提出：根據(jù)4個數(shù)據(jù)點（-3,3）（0,1）（2,1）（4,3），求解最小二乘拋物線。6.1 M次多項式擬合算法思想（6-1）本算法基于最小二乘原理來通過數(shù)據(jù)點來構(gòu)造M解最小二乘多項式。 以為例考慮求解最小二乘的拋物線擬合。設(shè)有個點，橫坐標的確定的。最小二乘拋物線的系數(shù)表示為（6-2）求解的線性方程組為（6-5）（6-4）（6-3）6.2 M次多項式擬合流程圖解線性方程組得到擬合多項式系數(shù)Start 輸入N個數(shù)據(jù)

17、點坐標和擬合多項式系數(shù)M用x、y的對應(yīng)次方構(gòu)造方程 End圖6-1 多項式擬合程序流程圖6.3 程序運行結(jié)果圖6-2 多項式擬合運行結(jié)果6.4 MATLAB擬合結(jié)果圖6-3 MATLAB擬合結(jié)果7 雅克比迭代法7.1 雅克比迭代法算法說明步1 取初始點，精度要求，最大迭代次數(shù)N,置步2 由計算步3 若則停算，輸出作為方程組的近似解；步4 若則停算，輸出迭代失敗信息；否則置，轉(zhuǎn)步2.7.2 雅克比迭代法流程圖開始 請輸入精度要求和最大迭代次數(shù) 計算判斷 否 是結(jié)束 圖7-1 雅克比迭代程序流程圖7.3 雅克比迭代法程序運行結(jié)果圖7-2 雅克比迭代運行結(jié)果8 高斯-賽德爾迭代法8.1 算法說明步1

18、 取初始點，精度要求，最大迭代次數(shù)N,置步2計算步3 若則停算，輸出作為方程組的近似解；否則轉(zhuǎn)步 4.步4若則停算，輸出迭代失敗信息；否則置，轉(zhuǎn)步2.8.2 高斯-賽德爾迭代法流程圖開始 請輸入精度要求和最大迭代次數(shù) 結(jié)束圖8-1 高斯賽德爾迭代流程圖8.3 高斯-賽德爾迭代法流運行結(jié)果圖8-2 高斯賽德爾迭代運行結(jié)果9 復(fù)化辛普森公式求積分9.1 復(fù)化辛普森公式求積分算法說明將積分區(qū)間a,b分成n等分，分點為其中記區(qū)間的中點為，在每個小區(qū)間上用辛普森公式，則可得到復(fù)化辛普森公式9.2 復(fù)化辛普森公式求積分流程圖開始定義輸入n,a,b輸出結(jié)束圖9-1 復(fù)化辛普森流程圖9.3 復(fù)化辛普森公式運行

19、結(jié)果圖9-2 復(fù)化辛普森運行結(jié)果10 最小二乘線性擬合10.1 最小二乘線性擬合算法說明（10-1）設(shè)采集到一組N個樣本數(shù)據(jù)，對這組數(shù)據(jù)進行線性擬合，求出擬合值。由每個樣本數(shù)據(jù)做誤差求和計算：，err值由值確定，因此，誤差和分別做偏微分計算，當每個偏微分值為0時，誤差和最小。故有：（10-2）（10-3） 根據(jù)每個偏微分等于零，可以得到兩個關(guān)于的一次方程，求解該方程即可以分別得到的值。10.2 最小二乘法進行線性擬合流程圖Start 輸入自變量x與因變量y，數(shù)據(jù)點的坐標用x、y的對應(yīng)次方構(gòu)造方程解線性方程組得到a，bEnd圖10-1 最小二乘進行線性擬合流程圖10.3 最小二乘法進行線性擬合運

20、行結(jié)果圖10-2 最小二乘進行線性擬合運行結(jié)果11 設(shè)計體會及改進意見這次課程設(shè)計是我受益匪淺，掌握了很多的東西，不僅僅包括對于算法的進一步理解，還有對于c+程序編程能力的提高，以及對于工作軟件word的操作能力。在上課期間，由于理論知識比較空泛、抽象，對于常用的算法沒有理解透徹。所以，這次課程設(shè)計對我的幫助很大，通過編程更深一步加深對算法的理解與掌握，提高了我的上機實踐能力。在實驗的編程過程中遇到了很大的困惑，由于對算法理論的不熟悉以及C+學(xué)習(xí)的知識有限，編寫的程序少且僅僅是一些小的算法。所以當數(shù)值算法與C+語言結(jié)合在一起的時候顯得有些吃力。在最后的排版中也是非常的吃力，由于平時對于word

21、操作不熟練導(dǎo)致這次課設(shè)進度緩慢。在這次課程設(shè)計編寫的10個程序中存在一個問題：有的程序在解決數(shù)學(xué)問題時，是在程序中錄入方程或者數(shù)據(jù)，而不是在運行界面輸入。這樣就會導(dǎo)致在更改問題之后，就必須在程序代碼中修改，以至于程序存在一定的局限性。這點是需要改進的。通過這次課程設(shè)計我深深地體會到數(shù)學(xué)方法和數(shù)學(xué)邏輯思維在計算機編程上的重要性，學(xué)好編程不僅僅是如何巧妙地實現(xiàn)問題，數(shù)學(xué)邏輯培養(yǎng)也是不容忽視的。參 考 文 獻1 譚浩強.C+程序設(shè)計.M.北京：清華大學(xué)出版社,2005.2 譚浩強.C+程序設(shè)計與上機指導(dǎo)(第三版).M.北京：清華大學(xué)出版社,2005.3數(shù)值方法（MATLAB版）馬昌鳳，林偉川編著.北

22、京：科學(xué)出版社，2008.附 錄附錄一 高斯列主元法程序代碼#include#include #includeusing namespace std;#define m 10/double Amm,bm,Augmm+1; /全局數(shù)組void equ(double Amm+1,double xm,int n) int i,i1,j,k; double Augmm+1,maxele,Temp,l,s;for(i1=0;i1n-1;i1+)maxele=fabs(Augi1i1);k=i1;for(i=i1;in;i+) /找主元if(maxelefabs(Augii1)maxele=fabs(Au

23、gii1);k=i; cout檢驗找主元的正確性：主元位于第 k+1 行，第 i1+1 列 endlendl;for(j=i1;jn+1;j+) /實施逐個元素換行Temp=Augi1j;Augi1j=Augkj;Augkj=Temp; /向用戶輸出提示信息 cout交換第 i1+1 行與第 k+1 行,請檢驗換行的正確性 Aug= endlendl; for(i=0;in;i+) /輸出Aug，以供檢驗 for(j=0;jn+1;j+)coutsetw(10)Augij ; coutendl; coutendl;for(k=i1+1;kn;k+) /消元，以第i1行為工具行處理以下各行，下面

24、代碼描述了處理第k行的過程l=-Augki1/Augi1i1;for(j=i1;jn+1;j+) Augkj=Augkj+l*Augi1j;xn-1=Augn-1n/Augn-1n-1; /回帶求解cout在被調(diào)函數(shù) equ 中，回帶計算解X，邊計算，邊輸出endlendl;cout xn=xn-1=0;i-)s=0;for(j=i+1;jn;j+)s=s+Augij*xj;xi=(Augin-s)/Augii; cout xi+1=xi; ; coutendl; /編寫主函數(shù)int main() void equ(double Amm+1,double xm,int n);int i,j,n

25、;static double Amm+1,bm,xm; cout請輸入未知量的個數(shù)nn;if(nm)cout問題規(guī)模太大，需要修改源程序中等符常量m定義endlendl;return 0;for(i=0;in;i+)cout請輸入增廣矩陣A的第i+1行:;for(j=0;jAij; /調(diào)用equ 函數(shù)解 n 元一次線性方程組 equ(A,x,n);/輸出解Xcout在主函數(shù)中輸出解x= ;coutsetiosflags(ios:fixed)setiosflags(ios:right)setprecision(4);for(i=0;in;i+)coutsetw(10)xi , ;coutendl

26、endl;return 0;附錄二 牛頓法解非線性方程組的程序代碼#include#include#define N 2 / 非線性方程組中方程個數(shù)、未知量個數(shù) #define Epsilon 0.0001 / 差向量1范數(shù)的上限#define Max 100 /最大迭代次數(shù)using namespace std;const int N2=2*N;int main()void ff(float xxN,float yyN); /計算向量函數(shù)的因變量向量yyNvoid ffjacobian(float xxN,float yyNN); /計算雅克比矩陣yyNNvoid inv_jacobian(

27、float yyNN,float invNN); /計算雅克比矩陣的逆矩陣invvoid newdundiedai(float x0N, float invNN,float y0N,float x1N); /由近似解向量 x0 計算近似解向量 x1float x0N,y0N,jacobianNN,invjacobianNN,x1N,errornorm;int i,j,iter=0;cout初始近似解向量：endl;for (i=0;ix0i; coutendl;coutendl;doiter=iter+1;cout第 iter 次迭代開始endl; /計算向量函數(shù)的因變量向量 y0ff(x0,

28、y0); /計算雅克比矩陣 jacobianffjacobian(x0,jacobian); /計算雅克比矩陣的逆矩陣 invjacobianinv_jacobian(jacobian,invjacobian); /由近似解向量 x0 計算近似解向量 x1newdundiedai(x0, invjacobian,y0,x1); /計算差向量的1范數(shù)errornormerrornorm=0;for (i=0;iN;i+)errornorm=errornorm+fabs(x1i-x0i);if (errornormEpsilon) break; for (i=0;iN;i+)x0i=x1i; wh

29、ile (iterMax);return 0;void ff(float xxN,float yyN) /調(diào)用函數(shù)float x,y; int i; x=xx0; y=xx1; yy0=x*x-2*x-y+0.5; yy1=x*x+4*y*y-4; /計算初值位置的值 cout向量函數(shù)的因變量向量是： endl; for( i=0;iN;i+) coutyyi ; coutendl; coutendl;void ffjacobian(float xxN,float yyNN) float x,y; int i,j; x=xx0;y=xx1;/jacobian have n*n element

30、/計算函數(shù)雅克比的值yy00=2*x-2;yy01=-1;yy10=2*x;yy11=8*y;cout雅克比矩陣是： endl;for( i=0;iN;i+) for(j=0;jN;j+) coutyyij ; coutendl; coutendl;void inv_jacobian(float yyNN,float invNN)float augNN2,L; int i,j,k; for (i=0;iN;i+) for(j=0;jN;j+) augij=yyij; for(j=N;jN2;j+)if(j=i+N) augij=1;else augij=0;for (i=0;iN;i+) fo

31、r (k=i+1;kN;k+) L=-augki/augii;for(j=i;j0;i-) for (k=i-1;k=0;k-) L=-augki/augii;for(j=N2-1;j=0;j-) augkj=augkj+L*augij; for (i=N-1;i=0;i-) for(j=N2-1;j=0;j-)augij=augij/augii;for (i=0;iN;i+) for(j=N;jN2;j+) invij-N=augij;cout雅克比矩陣的逆矩陣： endl;for (i=0;iN;i+) for(j=0;jN;j+) coutinvij ; void newdundieda

32、i(float x0N, float invNN,float y0N,float x1N)int i,j;float sum=0;for(i=0;iN;i+) sum=0; for(j=0;jN;j+) sum=sum+invij*y0j;x1i=x0i-sum; coutendl; cout近似解向量：endl;for (i=0;iN;i+) coutx1i ; coutendl;coutendl;附錄三 經(jīng)典四階龍格庫塔法#include#includeusing namespace std;double f(double x,double y,double z); /函數(shù)聲明double

33、 g(double x,double y,double z); /函數(shù)聲明int main() double z0,x0,y0; /定義變量 couty0z0; int n; coutn; double a,b,h; /求解的區(qū)間為0,5 a=0;b=5; h=(b-a)/n; cout輸出步長:h=hendl; double f1,f2,f3,f4,g1,g2,g3,g4; double x1,z1,y1; x0=0; for(int k=0;kn;k+) /龍格庫塔解一階線性方程組f1=f(x0,y0,z0);g1=g(x0,y0,z0);f2=f(x0+h/2,y0+h/2*f1,z0+

34、h/2*g1);g2=g(x0+h/2,y0+h/2*f1,z0+h/2*g1);f3=f(x0+h/2,y0+h/2*f2,z0+h/2*g2);g3=g(x0+h/2,y0+h/2*f2,z0+h/2*g2);f4=f(x0+h,y0+h*f3,z0+h*g3);g4=g(x0+h,y0+h*f3,z0+h*g3);x1=x0+h;x0=x1; y1=y0+h/6*(f1+2*f2+2*f3+f4);y0=y1;z1=z0+h/6*(g1+2*g2+2*g3+g4);z0=z1; coutx0=x0t; couty0=y0t; coutz0=z0endl;return 0; double

35、f(double x,double y,double z) /定義方程組第一個函數(shù) double r; r=-0.01*y-99.99*z; return r; double g(double x,double y,double z) /定義方程組第二個函數(shù) double p; p=-100*z; return z; 附錄四 三次樣條插值算法（壓緊樣條）#include#includeusing namespace std;const int max = 50;float xmax, ymax, hmax;float cmax, amax, fxymmax;float f(int x1, in

36、t x2, int x3)float a=(yx3-yx2) / (xx3 - xx2);float b = (yx2 - yx1) / (xx2 - xx1);return (a-b)/(xx3-xx1); /求差分void cal_m(int n) float Bmax;B0=c0 / 2;for(int i=1; in; i+)Bi=ci/(2-ai*Bi-1);fxym0=fxym0/2;for(i=1;i=0;i-)fxymi=fxymi-Bi*fxymi+1;void printout(int n);int main()int n,i; char ch;do coutn;for(i =0;i=n;i+)cout請輸入Xixi; cout請輸入Yiyi; for(i =0; i n; i+) /求步長hi=xi+1-xi;coutt;switch(t) case 1:cout輸入 Y0 Ynf0f1;c0 = 1; an = 1;fxym0 =