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

1、數(shù)值計算課程設計說明書題目： 典型數(shù)值算法的C+語言程序設計 院 系： 理學院 專業(yè)班級： 信息131 學 號： 202112030120 學生姓名： 郝騰宇 指導教師： 劉海峰 2021年 7月 10 日陜 西 科 技 大 學數(shù)值計算課程設計任務書理學院 信息與計算科學 專業(yè) 信息131 班級 學生： 郝騰宇 題目：典型數(shù)值算法的C+語言程序設計 課程設計從 2021 年 5 月 20 日起到 2021 年 7月 10 日1、課程設計的內(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、對課程設計成果的要求包括圖表、實物等硬件要求：2.1 提交課程設計報告按照算法要求，應用C+語言設計和開發(fā)算法程序，提交由：每個算法的原理與公式說明；與算法相應的程序設計說明程序中的主要變量語義說明，變量的數(shù)據(jù)類型，數(shù)據(jù)在內(nèi)存中組織和存儲結(jié)構(gòu)說明，各函數(shù)模塊的主要流程

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

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

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

6、算法圖。b表格：表的名稱及表內(nèi)文字采用中文，中文字體為五號宋體，表號表名在表格上面。表號以章為單位順序編號，表內(nèi)必須按規(guī)定的符號標注單位。格式為：表1-1，空一字符后接表格名稱。比方第4章第1個表是關(guān)于三次樣條插值的插值點列表，表的上方表號表名那么應為：表4-1 插值點的列表。c公式：公式書寫應在文中另起一行，居中排列。公式序號按章順序編號。字體為五號宋體，序號靠頁面右側(cè)。比方第3章第1個公式其編號那么應為：3-1。設計體會及今后的改良意見：設計總結(jié)要寫出算法理解，編程經(jīng)驗等技術(shù)性、學術(shù)性總結(jié)；體會要簡潔、真實、深刻，切忌空話、大話，客套話和矯揉造作之詞。改良意見要合理、中肯。參考文獻的要求：

7、另起一頁，居中打印參考文獻四字四號黑體，段前段后1行，字間空一字符；另起一行，按論文中參考文獻出現(xiàn)的先后順序用阿拉伯數(shù)字連續(xù)編號參考文獻編號應在正文中標注出；參考文獻中每條工程應齊全字體均為小四號宋體。格式：編號作者.論文或著作名稱.期刊名或出版社.出版時間。期刊應注明第幾期、起止頁數(shù)包括論著。參考文獻中條目要符合科技文獻引用文獻條目書寫的國家標準標準。2.3 設計報告裝訂順序與標準封面數(shù)值計算課程設計任務書目錄數(shù)值計算課程設計報告正文設計體會及今后的改良意見參考文獻(無需加目錄序號)附錄無需加目錄序號)左邊緣裝訂3、課程設計工作進度方案：時間設計任務及要求第18周編寫和調(diào)試程序并按要求撰寫設

8、計報告 目 錄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 經(jīng)典四階龍格庫塔法解一階微分方程運行結(jié)果94 三次樣條插值算法壓緊樣條104.1 三次樣條插值算法說明104.

9、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 高斯-賽德爾迭代法218.1 算法說明218.2 高斯-賽德爾迭代法流程圖228.3 高斯-賽德爾迭代法流運行結(jié)果239 復化辛普森公式

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

11、附錄九 復化辛普森求積公式程序代碼55附錄十 最小二乘法線性擬合程序代碼571 高斯列主元法解線性方程組問題提出：求解線性方程組1.1 高斯列主元法解線性方程組算法思想1） 將線性方程組寫成增廣矩陣的形式；2） 對增廣矩陣進行行變換，對元素，在第i列中及以下的元素選取絕對值最大的元素，將所有元素中最大的所在的行與第i行交換，然后采用高斯消元法使得新得到的第i行以下的元素均為零。一直重復上述過程直到。從而得到上三角矩陣。3對上三角矩陣進行回代求解，即可以得到方程組的解。1.2 高斯列主元法解線性方程組流程圖輸入增廣矩陣Start 輸入未知量的個數(shù)n高斯列主元法變換為上三角矩陣矩陣是否退化YN回代

12、求解End 圖1-1 高斯列主元法流程圖1.3 高斯列主元法解線性方程組運行結(jié)果圖1-2 高斯列主元法程序運行界面圖2 牛頓法解非線性方程組問題提出：求解非線性方程組設初值，用牛頓法計算。2.1 牛頓法解非線性方程組算法思想前提：設1計算函數(shù)2-12計算雅可比矩陣2-2(2-3)3求線性方程組的解。(2-4)4計算下一點重復上述過程，直到到達要求的精度，停止計算。2.2 牛頓法解非線性方程組流程圖Start 輸入精度要求和最大迭代次數(shù)數(shù) 計算計算雅可比矩陣 =- 或者精度不夠Y EndN圖2-1 牛頓法解非線性方程組流程圖2.3 牛頓法解非線性方程組運行結(jié)果圖2-2 牛頓法解非線性方程組程序運

13、行圖2.4 牛頓法解非線性方程組手工計算求解非線性方程組設初值，用牛頓法計算。函數(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ù)盡可能高。我們用四個不同點上的函數(shù)值的線性組合，將精度提高到四階就可以得到四階龍格-庫塔公式。常用的標準四階龍格-庫塔公式為： (3-1)其中 (3-2)

14、經(jīng)典四階龍格庫塔法解一階微分方程流程圖 Start Y N輸出結(jié)果End圖3-1 四階龍格庫塔法程序流程圖3.3 經(jīng)典四階龍格庫塔法解一階微分方程運行結(jié)果圖3-2 四階龍格庫塔法程序運行界面4 三次樣條插值算法壓緊樣條4.1 三次樣條插值算法說明設及，插值函數(shù)在每個小區(qū)間上是不超過3次的多項式且具有2階連續(xù)導數(shù)，那么稱為三階樣條插值。具體的說，三階樣條插值是滿足以下條件的分段三次多項式：（1） 插值條件：(4-1)； (2)連接條件：(4-2)壓緊樣條公式(4-3)三階樣條插值的邊界條件為1兩端的一階導數(shù)；2兩端的二階導數(shù)；3默認:自然邊界條件。4.2 三次樣條插值算法流程圖 Start輸入所

15、有點坐標及邊界條件計算三次樣條插值系數(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) 計算 對計算終止條件一般取為手工計算過程如下：開始 輸入下限a、上限b、精度。 ,，計算 得到積分值結(jié)束圖5-1 龍貝格積分程序流程圖圖5-2 龍貝格積分程序運行圖6 M

16、次多項式擬合問題提出：根據(jù)4個數(shù)據(jù)點-3,30,12,14,3，求解最小二乘拋物線。6.1 M次多項式擬合算法思想6-1本算法基于最小二乘原理來通過數(shù)據(jù)點來構(gòu)造M解最小二乘多項式。 以為例考慮求解最小二乘的拋物線擬合。設有個點，橫坐標確實定的。最小二乘拋物線的系數(shù)表示為6-2求解的線性方程組為6-56-46-36.2 M次多項式擬合流程圖解線性方程組得到擬合多項式系數(shù)Start 輸入N個數(shù)據(jù)點坐標和擬合多項式系數(shù)M用x、y的對應次方構(gòu)造方程 End圖6-1 多項式擬合程序流程圖6.3 程序運行結(jié)果圖6-2 多項式擬合運行結(jié)果6.4 MATLAB擬合結(jié)果圖6-3 MATLAB擬合結(jié)果7 雅克比迭

17、代法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 取初始點，精度要求，最大迭代次數(shù)N,置步2計算步3 假設那么停算，輸出作為方程組的近似解；否那么轉(zhuǎn)步 4.步4假設那么停算，輸出迭代失敗信息；否那么置，轉(zhuǎn)步2.8.2 高斯-賽德爾迭代法流程圖開

18、始 請輸入精度要求和最大迭代次數(shù) 結(jié)束圖8-1 高斯賽德爾迭代流程圖8.3 高斯-賽德爾迭代法流運行結(jié)果圖8-2 高斯賽德爾迭代運行結(jié)果9 復化辛普森公式求積分9.1 復化辛普森公式求積分算法說明將積分區(qū)間a,b分成n等分，分點為其中記區(qū)間的中點為，在每個小區(qū)間上用辛普森公式，那么可得到復化辛普森公式9.2 復化辛普森公式求積分流程圖開始定義輸入n,a,b輸出結(jié)束圖9-1 復化辛普森流程圖9.3 復化辛普森公式運行結(jié)果圖9-2 復化辛普森運行結(jié)果10 最小二乘線性擬合10.1 最小二乘線性擬合算法說明10-1設采集到一組N個樣本數(shù)據(jù)，對這組數(shù)據(jù)進行線性擬合，求出擬合值。由每個樣本數(shù)據(jù)做誤差求和

19、計算：，err值由值確定，因此，誤差和分別做偏微分計算，當每個偏微分值為0時，誤差和最小。故有：10-210-3 根據(jù)每個偏微分等于零，可以得到兩個關(guān)于的一次方程，求解該方程即可以分別得到的值。10.2 最小二乘法進行線性擬合流程圖Start 輸入自變量x與因變量y，數(shù)據(jù)點的坐標用x、y的對應次方構(gòu)造方程解線性方程組得到a，bEnd圖10-1 最小二乘進行線性擬合流程圖10.3 最小二乘法進行線性擬合運行結(jié)果圖10-2 最小二乘進行線性擬合運行結(jié)果11 設計體會及改良意見這次課程設計是我受益匪淺，掌握了很多的東西，不僅僅包括對于算法的進一步理解，還有對于c+程序編程能力的提高，以及對于工作軟件

20、word的操作能力。在上課期間，由于理論知識比擬空泛、抽象，對于常用的算法沒有理解透徹。所以，這次課程設計對我的幫助很大，通過編程更深一步加深對算法的理解與掌握，提高了我的上機實踐能力。在實驗的編程過程中遇到了很大的困惑，由于對算法理論的不熟悉以及C+學習的知識有限，編寫的程序少且僅僅是一些小的算法。所以當數(shù)值算法與C+語言結(jié)合在一起的時候顯得有些吃力。在最后的排版中也是非常的吃力，由于平時對于word操作不熟練導致這次課設進度緩慢。在這次課程設計編寫的10個程序中存在一個問題：有的程序在解決數(shù)學問題時，是在程序中錄入方程或者數(shù)據(jù)，而不是在運行界面輸入。這樣就會導致在更改問題之后，就必須在程序

21、代碼中修改，以至于程序存在一定的局限性。這點是需要改良的。通過這次課程設計我深深地體會到數(shù)學方法和數(shù)學邏輯思維在計算機編程上的重要性，學好編程不僅僅是如何巧妙地實現(xiàn)問題，數(shù)學邏輯培養(yǎng)也是不容無視的。參 考 文 獻1 譚浩強.C+程序設計.M.北京：清華大學出版社,2005.2 譚浩強.C+程序設計與上機指導(第三版).M.北京：清華大學出版社,2005.3數(shù)值方法MATLAB版馬昌鳳，林偉川編著.北京：科學出版社，2021.附 錄附錄一 高斯列主元法程序代碼#include#include #includeusing namespace std;#define m 10/double Amm,

22、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(Augii1);k=i; cout檢驗找主元的正確性：主元位于第 k+1 行，第 i1+1 列 endlendl;for(j=i1;jn+1;j+) /實施逐個元素換行Temp=Augi1j;Augi1

23、j=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行為工具行處理以下各行，下面代碼描述了處理第k行的過程l=-Augki1/Augi1i1;for(j=i1;jn+1;j+) Augkj=Augkj+l*Augi1j;xn-1=Augn-1n/Augn-1n-1; /回帶求解

24、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;static double Amm+1,bm,xm; cout請輸入未知量的個數(shù)nn;if(nm)cout問題規(guī)模太大，需要修改源程序中等符常量m定義endlendl;return 0;for(i

25、=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 , ;coutendlendl;return 0;附錄二 牛頓法解非線性方程組的程序代碼#include#include#define N 2 / 非線性方程組中方程個數(shù)、未知量個數(shù) #define Epsilon 0.

26、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(float yyNN,float invNN); /計算雅克比矩陣的逆矩陣invvoid newdundiedai(float x0N, float invNN,float y0N,float x1

27、N); /由近似解向量 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,y0); /計算雅克比矩陣 jacobianffjacobian(x0,jacobian); /計算雅克比矩陣的逆矩陣 invjacobianinv_jacobian(jacobian,invjac

28、obian); /由近似解向量 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; while (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*

29、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 /計算函數(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+) cout

30、yyij ; 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+) for (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-)

31、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 newdundiedai(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

32、+) 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 g(double x,double y,double z); /函數(shù)聲明int main() double z0,x0,y0; /定義變量 couty0z0; int n; coutn; doub

33、le 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+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+

34、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 f(double x,double y,double z) /定義方程組第一個函數(shù) double r; r=-0.01*y-99.99*z; return r; double g(double x,

35、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, int x2, int x3)float a=(yx3-yx2) / (xx3 - xx2);float b = (yx2 - yx1) / (xx2 - xx1);return (a-b)/(xx3-xx

36、1); /求差分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 = 6*(y1