西南交通大學(xué)數(shù)值分析上機(jī)實(shí)習(xí)

PAGEPAGE10序言數(shù)值分析是計算數(shù)學(xué)的一個主要部分，它不僅要研究各種數(shù)學(xué)問題的數(shù)值解法，同時也要分析所用的數(shù)值解法在理論上的合理性，如解法所產(chǎn)生的誤差能否滿足精度要求：解法是否穩(wěn)定、是否收斂及熟練的速度等。由于求解的迭代過程很復(fù)雜，計算量很大，所以需要借助其他輔助手段來求解。本文中采用了軟件C++來編寫求解程序，C++語言簡潔緊湊，使用方便靈活，其程序書寫形式自由，壓縮了一切不必要的編程工作并且控制語句非常豐富（如for循環(huán)，while循環(huán)，break語句和if語句），對編寫程序提供了很大的便利。在本文中使用C++編寫了牛頓法、牛頓-Steffensen法方程求解的程序和雅格比法、高斯－賽德爾迭代法求解方程組的程序及Runge-Kutta4階算法，并通過實(shí)例求解驗(yàn)證了其可行性，比較了求解同一種問題時不同方法之間的優(yōu)缺性，其中包含解的精確度和解的收斂速度兩個重要指標(biāo)。選用該軟件編寫主要原因是大學(xué)時期學(xué)習(xí)了C++程序設(shè)計，因此對于改軟件相對而言要熟悉些。目錄第一章牛頓法和牛頓-Steffensen法迭代求解的比較 31.1計算結(jié)果 31.2結(jié)果分析 4第二章Jacobi和Causs-Seidel迭代法比較 42.1計算結(jié)果 42.2結(jié)果分析 4第三章觀察Runge-Kutta4階算法穩(wěn)定區(qū)間的作用 53.1計算結(jié)果 63.2結(jié)果分析 6總結(jié) 6附錄 7第一章牛頓法和牛頓-Steffensen法迭代求解的比較分別用牛頓法，及基于牛頓算法下的Steffensen加速法求ln(x+sinx)=0的根。初值x0分別取0.1,1,1.5,2,4進(jìn)行計算。求sinx=0的根。初值x0分別取1，1.4，1.6,1.8，3進(jìn)行計算。分析其中遇到的現(xiàn)象與問題。1.1計算結(jié)果(1)方程可變形為x+注：x0為初始值，x為收斂解，N為迭代次數(shù)牛頓法x00.111.524x0.5109730.5109730.5109730.5109730.510973N445640牛頓—Steffensen加速法x0.5109730.510973無收斂解0.5109730.510973N3348(2)牛頓法x011.41.61.83x03.1415931.41596.283193.14159N57843牛頓—Steffensen加速法x無收斂解-3.1415925.13276.28319無收斂解N4631.2結(jié)果分析從結(jié)果對比我們可發(fā)現(xiàn)牛頓—Steffensen加速法比牛頓法要收斂的快，牛頓法對于初值的選取特別重要，比如第（1）問中的初值為4的情況，迭代次數(shù)算了40次，遠(yuǎn)大于其余初值的情況；在第（2）問中的初值為1.6的情況，收斂解得31.4159，分析其原因應(yīng)該是，x0=1.6,；在牛頓—Steffensen加速法第(1)問中x=1.5、第二問x=0和3的情況，無收斂解。第二章Jacobi和Causs-Seidel迭代法比較用雅格比法與高斯－賽德爾迭代法解下列方程組Ax=b，研究其收斂性，上機(jī)驗(yàn)證理論分析是否正確，比較它們的收斂速度，觀察右端項(xiàng)對迭代收斂有無影響。(1)A行分別為A1=[6,2,-1],A2=[1,4,-2],A3=[-3,1,4]；b1=[-3,2,4]T,b2=[100,-200,345]T,(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,-10]T,(3)A行分別為A1=[1,3],A2=[-7,1]；b=[4,6]T,2.1計算結(jié)果（1）x0=[0,0,0]T為初始值,N為迭代次數(shù)，0.00001為誤差精度,X為收斂解Jacobi迭代Causs-Seidel迭代Ax=b1X[-0.727271,0.808083,0.252524][-0.727275,0.808082,0.252523]N2012Ax=b2X[36.3636,-2.0707,114.04][36.3636,-2.07071,114.04]N2617（2）Jacobi迭代Causs-Seidel迭代Ax=b1X無收斂解[5.76922,0.76922,-4.23075]N37Ax=b2X無收斂解[36.6923,7.69229,-42.3077]N44（3）Jacobi迭代Causs-Seidel迭代Ax=bX無收斂解無收斂解N2.2結(jié)果分析（1）第一小題的雅可比迭代矩陣為0-0.33330.1667-0.250000.50000.7500-0.25000其特征值分別為-0.5427，0.2713+0.3708i，0.2713-0.3708i，經(jīng)計算譜半徑為小于1,故方程組雅可比迭代收斂。而高斯－賽德爾迭代矩陣為0-0.33330.166700.08330.45830-0.27080.0104其特征值為0，0.0469+0.3504i，0.0469-0.3504i經(jīng)計算譜半徑為0.3535小于1,故原方程組高斯－賽德爾迭代矩陣收斂。第二小題的雅可比迭代矩陣為0-0.8000-0.8000-0.80000-0.8000-0.8000-0.80000征值分別為-1.6，0.8，0.8，經(jīng)計算譜半徑為1.6大于1,故方程組雅可比迭代發(fā)散。而高斯－賽德爾迭代矩陣為0-0.8000-0.800000.6400-0.160000.12800.7680其特征值為0，0.7040+0.1280i，0.704280i經(jīng)計算譜半徑為0.7155小于1,故原方程組高斯－賽德爾迭代矩陣收斂。第三小題的雅可比迭代矩陣為0-370其特征值分別為4.5826i，－4.5826i，譜半徑為4.5826大于1,故方程組雅可比迭代發(fā)散。而高斯－賽德爾迭代矩陣為0-30-21其特征值為0，－21，譜半徑為21大于1,故原方程組高斯－賽德爾迭代矩陣發(fā)散。（2）從2.1中的結(jié)果列表中可以看到，Causs-Seidel迭代法比Jacobi迭代法收斂速度要快，并且方程組右端項(xiàng)對迭代收斂是無影響的。第三章觀察Runge-Kutta4階算法穩(wěn)定區(qū)間的作用用Runge-Kutta4階算法對初值問題y/=-20*y，y(0)=1按不同步長求解，用于觀察穩(wěn)定區(qū)間的作用，推薦兩種步長h=0.1,0.2。注：此方程的精確解為：3.1計算結(jié)果=1.0h=0.1h=0.21.6935131251.693331253.2結(jié)果分析,h=0.2時，h=-4,而Runge-Kutta4階算法的絕對穩(wěn)定區(qū)間是[-2.78,0],故h=0.2時計算不穩(wěn)定；而h=0.1時，h=-2，在絕對穩(wěn)定區(qū)間內(nèi)，計算穩(wěn)定，結(jié)果可靠?？偨Y(jié)通過這次上機(jī)練習(xí)，讓我對數(shù)值分析所介紹的迭代求解方法及其理論有了更深層次的理解，了解了各種方法之間的優(yōu)缺點(diǎn)，并且認(rèn)識到了自己在以前的學(xué)習(xí)中所存在的問題，及時的修補(bǔ)了自己知識上的漏洞。同時也提高了我在編寫程序上的熟練程度，所以，我認(rèn)為這次上機(jī)實(shí)習(xí)是非常有收獲的，給予我學(xué)習(xí)數(shù)值分析的幫助也是非常大的。附錄1.1第(1)問牛頓法#include<iostream>#include<cmath>usingnamespacestd;intmain(){ doublex,x0; inti=0; cout<<"請輸入初值:"; cin>>x0; do{x=x0; x0=x0-(x0+sin(x0)-1)/(1+cos(x0)); i++;}while(fabs(x-x0)>=0.000001);cout<<"x="<<x0<<endl;cout<<"N="<<i<<endl;return0;}1.2第(1)問牛頓—Steffensen加速法#include<iostream>#include<cmath>usingnamespacestd;intmain(){ doublex,x0,x1,x2; inti=0; cout<<"請輸入初值:"; cin>>x0; do{ x=x0;x1=x0-(sin(x0)+x0-1)/(cos(x0)+1); x2=x1-(sin(x1)+x1-1)/(cos(x1)+1); x0=x0-(x1-x0)*(x1-x0)/(x2-2*x1+x0); i++;}while(fabs(x-x0)>=0.000001);cout<<"x="<<x0<<endl;cout<<"N="<<i<<endl;return0;}1.3第(2)問牛頓法#include<iostream>#include<cmath>usingnamespacestd;intmain(){ doublex,x0; inti=0; cout<<"請輸入初值:"; cin>>x0; do{x=x0; x0=x0-sin(x0)/cos(x0); i++;}while(fabs(x-x0)>=0.000001);cout<<"x="<<x0<<endl;cout<<"N="<<i<<endl;return0;}1.4第(2)問牛頓—Steffensen加速法#include<iostream>#include<cmath>usingnamespacestd;intmain(){ doublex,x0,x1,x2; inti=0; cout<<"請輸入初值:"; cin>>x0; do{ x=x0;x1=x0-sin(x0)/cos(x0); x2=x1-sin(x1)/cos(x1); x0=x0-(x1-x0)*(x1-x0)/(x2-2*x1+x0); i++;}while(fabs(x-x0)>=0.000001);cout<<"x="<<x0<<endl;cout<<"N="<<i<<endl;return0;}2.1雅格比法#include<iostream>#include<cmath>usingnamespacestd;intmain(){ doublea[10][10],b[10],sum=0.0,x[10],y[10],c,s=0.0; inti,k=0,n,j; cout<<"請輸入維數(shù):"; cin>>n; cout<<"請輸入數(shù)組A:"; for(i=0;i<n;i++) for(j=0;j<n;j++) cin>>a[i][j];cout<<"請輸入數(shù)組B:"; for(i=0;i<n;i++) cin>>b[i];cout<<"請輸入初始X:"; for(i=0;i<n;i++) cin>>x[i]; do{for(i=0;i<n;i++) {for(j=0;j<i;j++) sum=sum+a[i][j]*x[j];for(j=i+1;j<n;j++) sum=sum+a[i][j]*x[j]; y[i]=(b[i]-sum)/a[i][i]; sum=0; } c=0; for(i=0;i<n;i++) if(c<fabs(x[i]-y[i])) c=fabs(x[i]-y[i]); if(c>0.00001) c=0; elses=1;for(i=0;i<n;i++) x[i]=y[i];k++;}while(s!=1);for(i=0;i<n;i++) cout<<"解為："<<x[i]<<endl;cout<<"N="<<k<<endl;return0;}2.2高斯－賽德爾迭代法#include<iostream>#include<cmath>usingnamespacestd;intmain(){ doublea[10][10],b[10],sum=0.0,x[10],y[10],c,s=0.0,z[10]; inti,k=0,n,j; cout<<"請輸入維數(shù):"; cin>>n; cout<<"請輸入數(shù)組A:"; for(i=0;i<n;i++) for(j=0;j<n;j++) cin>>a[i][j];cout<<"請輸入數(shù)組B:"; for(i=0;i<n;i++) cin>>b[i];cout<<"請輸入初始X:"; for(i=0;i<n;i++) cin>>x[i]; do{for(i=1;i<n;i++) z[i]=x[i]; for(i=0;i<n;i++) {for(j=0;j<i;j++) sum=sum+a[i][j]*x[j];for(j=i+1;j<n;j++) sum=sum+a[i][j]*x[j]; y[i]=(b[i]-sum)/a[i][i]; sum=0; x[i]=y[i]; } c=0; for(i=0;i<n;i++) if(c<fabs(z[i]-y