數(shù)值計算課程設(shè)計實驗合并

1、本科生實驗報告實驗課程 數(shù)值計算方法 學院名稱 信息科學與技術(shù)學院 專業(yè)名稱 計算機科學與技術(shù) 學生姓名 周瑞豪 學生學號 202113030317 指導教師 馬桂媛 實驗地點 6A502 實驗成績 二 15 年 4 月 二 15 年 5 月填寫說明1、 適用于本科生所有的實驗報告印制實驗報告冊除外；2、 專業(yè)填寫為專業(yè)全稱，有專業(yè)方向的用小括號標明；3、 格式要求： 用A4紙雙面打印封面雙面打印或在A4大小紙上用藍黑色水筆書寫。 打印排版：正文用宋體小四號，1.5倍行距，頁邊距采取默認形式上下2.54cm，左右2.54cm，頁眉1.5cm，頁腳1.75cm。字符間距為默認值縮放100%，間距

2、：標準；頁碼用小五號字底端居中。 具體要求：題目二號黑體居中；摘要“摘要二字用小二號黑體居中，隔行書寫摘要的文字局部，小4號宋體；關(guān)鍵詞隔行頂格書寫“關(guān)鍵詞三字，提煉3-5個關(guān)鍵詞，用分號隔開，小4號黑體)； 正文局部采用三級標題；第1章 ××(小二號黑體居中，段前0.5行)1.1 ×××××小三號黑體×××××段前、段后0.5行小四號黑體段前、段后0.5行參考文獻黑體小二號居中，段前0.5行，參考文獻用五號宋體，參照?參考文獻著錄規(guī)那么GB/T 77142005?。實驗一

3、非線性方程求根第1章 問題描述實驗目的：掌握非線性方程求根的根本步驟及方法，。實驗內(nèi)容：試分別用二分法、簡單迭代法、Newton迭代法、弦截法(割線法、雙點弦法)，求 x5-3x3+x-1= 0 在區(qū)間 -8,8上的全部實根，誤差限為10-6。要求：討論求解的全過程，對所用算法的局部收斂性，優(yōu)缺點等作分析及比擬， 第2章 算法思想 算法分析此局部摘抄自課本P23頁 設(shè)函數(shù)有根區(qū)間表示為。將有根區(qū)間用中點分成兩 半，計算函數(shù)值。如果，就得到方程組的實根，否那么檢查與是否同號，假設(shè)同號，那么說明所求的根在的右側(cè)，這時令；否那么，根在的左側(cè)，這時令，這樣新的有根區(qū)間的長度為之半。對壓縮了的有根區(qū)間又

4、施以同樣的方法如此反復二分下去，即可得出一系列有根區(qū)間其中，每個區(qū)間都是前一個區(qū)間的一半，因此二分k次后的有根區(qū)間的長度為可見，如果二分過程無限地下去，這些有根區(qū)間最終必收縮于一點，該點顯然就是所求的根。取有根區(qū)間的中點作為根的近似值，此時的誤差假設(shè)事先給定的誤差要求為，那么只需便可以停止二分計算。 void divide(V a1,V a2) for(;fabs(fabs(a1)-fabs(a2)>=1e-6;) if(f(a1)*f(a1+a2)/2)<0) a2=(a1+a2)/2;else a1=(a1+a2)/2; cout<<endl<<&quo

5、t;經(jīng)過二分法求得的結(jié)果是："<<(a1+a2)/2<<endl<<endl; 2.2 簡單迭代法 算法分析此局部摘抄自肯課本P41面Newton下山法是擴大初值范圍的修正Newton法。將Newton迭代法的計算結(jié)果進行適當?shù)募訖?quán)平均作為新的改良值，即。從而化簡可得簡單法迭代公式 void newton(V a1) for(;fabs(fabs(a1)-fabs(a1-fnewton(a1)>=1e-6;) a1-=fnewton(a1); cout<<"經(jīng)過簡單迭代法得的結(jié)果是："<<a1<

6、;<endl<<endl; 2.3 Newton迭代法 算法分析(此局部摘抄自書上P37面) 對于非線性方程，假設(shè)根的一個近似值在這里是，將在處展成一階泰勒公式，忽略高次項，有。右端是直線方程，用這個直線方程來近似非線性方程。將非線性方程的根代入，即解得這就是Newton迭代公式。在具體的應(yīng)用中，寫為：，那么這樣獲得的即為按Newton迭代法求得的近似解。 void newton2(V a1) double t=1.0; for(;fabs(f(a1)<=fabs(f(a1-t*fnewton(a1);t/=2.0) for(;fabs(fabs(a1)-fabs(a1

7、-t*fnewton(a1)>=1e-6;) a1-=t*fnewton(a1); cout<<"使用Newton迭代法求得的結(jié)果是："<<a1<<endl<<endl; 算法分析此局部摘抄自課本P44面由于Newton迭代法有時收斂速度較慢，而且有時函數(shù)的一階導數(shù)不易求得或較為復雜，因此，改用兩個端點都在變動的弦，即用差商替代Newton迭代公式中的導數(shù)，從而導出。這就是雙點弦截法。雙點弦截法詳細的計算步驟為：1選定初始值，計算，。2按雙點弦法迭代公式計算，并求。3判斷：如果給定精度，那么迭代停止；否那么，用和分別代替

8、重復2和3。第3章 測試結(jié)果及分析 測試結(jié)果，要求按表格形式輸出每次迭代的結(jié)果表格格式如書上例題，并比擬不同方法的收斂速度、優(yōu)劣附錄：源程序清單#include ""#include <iostream> typedef double V;double finish(V x); double fnewton(V x); void divide(V a1,V a2); void newton(V a1); void cord(V a1,V a2); void newton2(V a1); using namespace std;double f(V x) retu

9、rn pow(x,5)-6*pow(x,3)+x-1; double fnewton(V x) if(5*pow(x,4)-9*pow(x,2)+1)!=0) return (pow(x,5)-3*pow(x,3)-1)/(5*pow(x,4)-9*pow(x,2)+1); elsereturn 1;void divide(V a1,V a2) for(;fabs(fabs(a1)-fabs(a2)>=1e-6;) if(f(a1)*f(a1+a2)/2)<0) a2=(a1+a2)/2;else a1=(a1+a2)/2; cout<<endl<<&quo

10、t;經(jīng)過二分法求得的結(jié)果是："<<(a1+a2)/2<<endl<<endl; void newton(V a1) for(;fabs(fabs(a1)-fabs(a1-fnewton(a1)>=1e-6;) a1-=fnewton(a1); cout<<"經(jīng)過Newton法得的結(jié)果是："<<a1<<endl<<endl;void cordline(V a1,V a2) int n=0; for(;fabs(f(a2-f(a2)*(a2-a1)/(f(a2)-f(a1)-0)

11、>=1e-6;n+) if(n=50) break; else if(f(a2-f(a2)*(a2-a1)/(f(a2)-f(a1)*f(a1)>0|(f(a2)-f(a1)!=0) a1=a2-f(a2)*(a2-a1)/(f(a2)-f(a1); else a2=a2-f(a2)*(a2-a1)/(f(a2)-f(a1); if(n=50) cout<<"使用雙點弦法求解失敗！"<<endl<<endl; elsecout<<"經(jīng)過雙點弦法求得的結(jié)果是："<<( a2-f(a2)

12、*(a2-a1)/(f(a2)-f(a1) )<<endl<<endl;void newton2(V a1) double t=1.0; for(;fabs(f(a1)<=fabs(f(a1-t*fnewton(a1);t/=2.0) for(;fabs(fabs(a1)-fabs(a1-t*fnewton(a1)>=1e-6;) a1-=t*fnewton(a1); cout<<"使用Newton下山法求得的結(jié)果是："<<a1<<endl<<endl;int main() double x

13、1,x2; cout<<"方程是：x5-3*x3+x-1= 0n請給出求根區(qū)間的左值和右值：n左值：" cin>>x1; cout<<"右值：" cin>>x2; if(f(x1)*f(x2)>0) cout<<"此區(qū)間上方程無根"<<endl; system("PAUSE"); return 0; divide(x1,x2); newton(x1+x2)/2); cordline(x1,x2); newton2(x1+x2)/2); s

14、ystem("PAUSE"); return 0; 實驗2 線性方程組數(shù)值解-直接法一、實驗目的：掌握Gauss消元法、列主元的Gauss消元法、Gauss-Jordan消元法、LU分解、改良平方根法等求解線性方程組的數(shù)值解的直接求解。二、實驗內(nèi)容：1. 用Gauss消元法、選主元的Gauss消元法求解以下線性方程組的解Ax=b，要求輸出消元過程中系數(shù)矩陣的變化情況2. 用Gauss-Jordan消元法解題1中的線性方程組，并求A的逆矩陣A-13. 對上題中的線性方程組的系數(shù)矩陣A做LUDoolittle分解，要求輸出L、U矩陣，并求線性方程組的解 (1) 高斯消元法算法分

15、析摘抄自課本在線性代數(shù)里就有用Gauss消元法線性方程組的解。通常，對線性方程組的增廣矩陣A進行初等行變換化為上三角矩陣的方法，就成為Gauss消元法。對于任意一矩陣，均可以用Gauss消元法求解對應(yīng)線性方程組的解。當用k表示消元過程的次序時，Gauss消元法的計算步驟為：1Gauss消元過程：設(shè)對計算2回代求解過程：(2) 列主元的Gauss消元法 算法分析列主元的Gauss消元法是在Gauss消元法的根底之上加上了列選主元的步2出絕對值最大的，然后通過行交換將其交換到的位置上。設(shè)主元在第個方程，即。假設(shè)，將和方程互易位置，使新的成為主元，然后繼續(xù)進行，這一步驟稱為列選主元。2.源代碼#in

16、clude<iostream>#include<iomanip>using namespace std;void DisplayMatrix(double *a,int size)cout<<endl<<"="<<endl;for(int i=0;i<size;i+)for(int j=0;j<size+1;j+)cout<<setw(15)<<setprecision(8)<<aij;cout<<endl;cout<<endl<<

17、"="<<endl;void Gauss(double *a,int size)int i,j,k;double tmp;cout<<"方程組對應(yīng)的增廣矩陣為:"<<endl;DisplayMatrix(a,size);/輸出矩陣,需自己實現(xiàn)該函數(shù) /*消元,第k列對角線以下全消成0*/for(k=0;k<size-1;k+)for(i=k+1;i<size;i+)tmp=aik/akk;for(j=k;j<=size;j+)aij-=tmp*akj;cout<<endl<<&

18、quot;第"<<k+1<<"列消元后的矩陣為："<<endl;DisplayMatrix(a,size);/*回代求解*/for(i=size-1;i>=0;i-)for(j=size-1;j>i;j-)aisize-=ajsize*aij;aisize/=aii;/*輸出方程組的解，放在數(shù)組最后一列中*/cout<<endl<<"方程的根為："<<endl;for(i=0;i<size;i+)cout<<"x"<&

19、lt;i+1<<"="<<setiosflags(ios:left)<<setw(15)<<setprecision(8)<<aisize;cout<<endl;int main() int n,i,j; cout<<"input the size of equations!"<<endl; cin>>n; double *p=new double *n; for( i=0;i<n;i+) pi=new double n+1; cout<

20、;<"Input the matrix of Equations!"<<endl; for(i=0;i<n;i+) for(j=0;j<=n;j+) cin>>pij; Gauss(p,n); for ( i=0;i<n;i+) delete pi; delete p; system("pause");截圖：實驗3 線性方程組數(shù)值解-迭代法一、實驗目的：掌握Jaccobi迭代法、Guass-Sidel迭代法、松弛法求解線性方程組的數(shù)值解。二、實驗內(nèi)容：1分別用雅格比法與高斯賽德爾迭代法解以下方程組Ax=b，

21、研究其收斂性，上機驗證理論分析是否正確，比擬它們的收斂速度，觀察右端項對迭代收斂有無影響。(1)A行分別為A1=6,2,-1,A2=1,4,-2,A3=-3,1,4； b1=-3,2,4T, b2=100,-200,345T,(2) A行分別為A1=1,0,8,0.8,A2 0.8,1,0.8,A3=0.8,0.8,1；b1=3,2,1 T, b2=5,0,-10T,(3)A行分別為A1=1,3,A2=-7,1；b=4,6T,2. 選作松弛因子對SOR法收斂速度的影響。用SOR法求解方程組Ax=b,其中要求程序中不存系數(shù)矩陣A,分別對不同的階數(shù)取進行迭代，記錄近似解x(k)到達|x(k)-x(

22、k-1)|<10-6時所用的迭代次數(shù)k，觀察松弛因子對收斂速度的影響，并觀察當w£0或w³2會有什么影響？1Jacobi迭代法摘抄自課本 算法分析 線性方程組Ax=b，系數(shù)矩陣A非奇異，且aii0。即 可以簡寫為從上式中別離出變量，將它寫成據(jù)此便建立了Jacobi迭代公式即可寫為這就是解線性方程組的Jacobi迭代公式。也是Jacobi迭代法求解線性方程組的理論根底。源代碼#include <iostream>#include<cmath>#include<iomanip>using namespace std;void jacco

23、bi(double *a,int size)cout<<"雅可比迭代法"double x100=0;double x_tp100=0;for(int i=0;i<17;i+)for(int i=0;i<size;i+) double tmp=aisize;for(int j=0;j<i;j+)tmp-=aij*xj;for(int j=i+1;j<size;j+)tmp-=aij*xj;x_tpi=tmp/aii;for(int i=0;i<size;i+) xi=x_tpi; cout<<endl<<&qu

24、ot;第"<<i+1<<"迭代：" for(int i=0;i<size;i+) cout<<"x"<<i+1<<"="<<left<<setw(9)<<setfill(' ')<<xi<<"|" cout<<endl; int main()int n,i,j;cout<<"請輸入維數(shù)："cin>>n;doub

25、le *mat=new double *n;for( i=0;i<n;i+) mati=new double n+1; cout<<"請輸入增廣矩陣："<<endl; for(i=0;i<n;i+) for(j=0;j<=n;j+) cin>>matij; jaccobi(mat,n);gs(mat,n);for ( i=0;i<n;i+) delete mati; delete mat; system("PAUSE");return 0;2Guass_Sidel迭代法(摘抄自課本)算法分析在G

26、uass_Sidel迭代時，每次迭代充分利用當前最新的迭代值。在迭代收斂時，因新值比老值更準確些，求出新值后，用代替前一次的迭代值繼續(xù)進行計算，這就充分利用新值建立起Guass_Sidel迭代公式。對于一般形式的方程組，Guass_Sidel迭代公式為源代碼#include <iostream>#include<cmath>#include<iomanip>using namespace std;void gs(double *a,int size)cout<<"高斯賽德爾迭代法"double x100=0;double x_

27、tp100=0;for(int k=0;k<17;k+)for(int i=0;i<size;i+)double tmp=aisize;for(int j=0;j<i;j+)tmp-=aij*xj;for(int j=i+1;j<size;j+)tmp-=aij*xj;xi=tmp/aii;/for(int i=0;i<size;i+)/xi=x_tpi;cout<<"第"<<right<<setw(2)<<i+1<<"迭代："for(int i=0;i<si

28、ze;i+)cout<<"x"<<i+1<<"="<<left<<setw(9)<<setfill(' ')<<xi<<"|"cout<<endl;int main()int n,i,j;cout<<"請輸入維數(shù)："cin>>n;double *mat=new double *n;for( i=0;i<n;i+)mati=new double n+1; cout&

29、lt;<"請輸入增廣矩陣："<<endl;for(i=0;i<n;i+)for(j=0;j<=n;j+)cin>>matij;jaccobi(mat,n);gs(mat,n);for ( i=0;i<n;i+)delete mati;delete mat;system("PAUSE");return 0;雅各布解的迭代過程a:高斯賽德爾迭代過程a：用雅各比經(jīng)過 19 步迭代得到解為，。用高斯-塞德爾算法經(jīng)過 11 步迭代得到的解為雅各布解的迭代過程b高斯賽德爾迭代過程b用雅各布經(jīng)過 24 步迭代得到方程組的

30、解是x1 =36.37, x2 =-2.071, x3 =114.譜半徑 (B) = 0.54用高斯-塞德爾經(jīng)過 16 迭代得到方程組的解是x12高斯賽德爾迭代過程a,b用雅各布算法求的譜半徑(B) = 1.6 > 1,由迭代矩陣的收斂條件知其不收斂。 用高斯-塞德爾經(jīng)過 34 步迭代得到方程組的解是x1 =5.7691, x2 =0.7692, x3 =-4.2307 譜半徑 (B) = 0.72 知用高斯-塞德爾法是收斂的。用雅各布用雅各布算法求的譜半徑(B) = 1.6 > 1,由迭代矩陣的收斂條件知其不收斂。7用高斯-塞德爾經(jīng)過 40 步迭代解方程組的解為x1譜半徑 (B)

31、 = 0.72 得用高斯-塞德爾法是收斂的。(3)用雅各比算法迭代得到(B) = 4.6 >1,故迭代過程不收斂。 用高斯-塞德爾算法迭代得(B) = 21>1,迭代過程不收斂。 (1)在同樣的精度要求和初始向量情況下,高斯-塞德爾的收斂速度要快于雅各比算法。 (2)迭代矩陣的收斂與否只和迭代矩陣B 的譜半徑有關(guān),同初始向量和方程組的右端項無關(guān) 。(3)一般情況下高斯-塞德爾收斂比雅各比快,但是也存在雅各比收斂而高斯-塞德爾不收 斂以及兩種算法均不收斂的情況。實驗4 Lagrange插值與三次樣條插值效果的比擬一、實驗目的：掌握Lagrange插值、三次樣條插值的原理，了解Rung

32、e現(xiàn)象及防止方法二、實驗內(nèi)容：實際問題中兩次過程記錄測得(x,y)數(shù)據(jù)點如下：其中yi=f(xi) 第一次的數(shù)據(jù)為：x=-5.0000 -4.0000 -3.0000 ;y=-0.0000 2 -0.0004 -0.0366 -0.3679 0 0.3679 2第二次的數(shù)據(jù)為： x= ;y =1203 321針對第一次采集的數(shù)據(jù)請分別采用Lagrange插值與三次樣條插值，估算在-5,5區(qū)間上各20等分點處的函數(shù)值，并將計算結(jié)果與第二次采集的結(jié)果進行比對，并分析誤差原因。要求：以表格形式輸出各計算結(jié)果，如：xyLarange(x)S(x)三、算法與步驟：1、拉格朗日插值：1輸入數(shù)據(jù)點總數(shù)為n+

33、1即輸入n的值，節(jié)點Xi和對應(yīng)的Yi的值i從0到n，令Lnx= 0；2當i從0到n時，計算Li(x)= (X-Xj)/(Xi-Xj)其中j從0到n且j不等于i；令s=1，j從0到n：如果j=i，那么s=s；否那么令s=s*(x-xj/(xi-xj);循環(huán)結(jié)束令Lix=s，Ln(x) = (Ln(x)+Li(x)*Yi； 3)將待估側(cè)值代入Lnx求值；2、三次樣條插值假定有n+1個數(shù)據(jù)節(jié)點a. 計算步長 (i = 0, 1, , n-1)b. 將數(shù)據(jù)節(jié)點和指定的首位端點條件帶入矩陣方程c. 解矩陣方程，求得二次微分值。該矩陣為三對角矩陣，具體求法參見我的上篇文章：三對角矩陣的求解。d.

34、 計算樣條曲線的系數(shù)：其中i = 0, 1, , n-1e. 在每個子區(qū)間 中，創(chuàng)立方程 四、實驗截圖：1、拉格朗日插值2、三次樣條插值附錄：源程序清單1、拉格朗日插值源碼#include<iostream>#include<iomanip>#define N 21 /插值節(jié)點數(shù)目using namespace std;void main()float aN; /所求點橫坐標，將插值節(jié)點的橫坐標復制float fx = 0, tmp = 1;int i, j;/差值節(jié)點橫坐標float xN = -5.0000, -4.5000, -4.0000,

35、-3.5000, -3.0000,-2.5000, -2.0000, -1.5000, -1.0000, -0.5000,0, 0.5000, 1.0000, 1.5000, 2.0000, 4.0000, 4.5000,5.0000;/差值節(jié)點縱坐標float yN = -0.0000, -0.0001, -0.0002, -0.0003, -0.0004, -0.0048, -0.0366, , 0, 0.3894, 0.3679, 0.1581, , 0.0004, 0.0003, 0.0002, 0.0001, 0.0000 ;cout << "插值節(jié)點橫坐標:n

36、"for (int i = 0; i < N; i+)cout << xi << " "cout << "n"cout << "插值節(jié)點縱坐標：n"for (int i = 0; i < N; i+)cout << yi << " "cout << "n"cout << "拉格朗日插值的結(jié)果為：" << endl;cout << &qu

37、ot;x"<< setw(20) << " y" << setw(20) << " Larange(x)" << endl;for (int i = 0; i < N; i+)ai = xi;for (int k = 0; k < N; k+)for (i = 0; i < N; i+)tmp = 1; /開始前復原tmpfor (j = 0; j < N; j+)if (i != j)tmp = tmp*(ak-xj) / (xi - xj);fx = fx

38、+ tmp*yi;cout <<xk<<setw(20)<<yk<<setw(20)<< fx << endl;fx = 0;2、三次樣條插值源碼#include<iostream>#include<iomanip>#include<cmath>#define N 21 /插值節(jié)點數(shù)目using namespace std; int main()/差值節(jié)點橫坐標double xN = -5.0000, -4.5000, -4.0000, -3.5000, -3.0000,-2.5000,

39、 -2.0000, -1.5000, -1.0000, -0.5000,0, 0.5000, 1.0000, 1.5000, 2.0000, 4.0000, 4.5000,5.0000;/差值節(jié)點縱坐標double yN = -0.0000, -0.0001, -0.0002, -0.0003, -0.0004, -0.0048, -0.0366, , 0, 0.3894, 0.3679, 0.1581, , 0.0004, 0.0003, 0.0002, 0.0001, 0.0000 ;int choice = 0;double *xx,*a, *b, *a1, *b1, *h, *m;xx = new doubleN;a = new doubleN;b = new doubleN;a1 = new doubleN;b1 = new doubleN;h = new doubleN-1;m = new doubleN+1;for (int j = 0; j < N - 1; j+)hj = xj + 1 - xj;cout << "三次樣條插值結(jié)果為：n