計算數值實驗報告(太原理工大學)

1、本科實驗報告課程名稱： 計算數值方法 實驗地點： 綜合樓五層506室 專業(yè)班級：計科1002 學號： 2010001414 學生姓名： xxx 指導教師： 王崢 2012 年 6 月 20 太原理工大學學生實驗報告學院名稱計算機科學與技術專業(yè)班級計科 1002學號2010001414學生姓名 xxx實驗日期2012.6.5成績課程名稱計算數值方法實驗題目實驗一 方程求根一、課題名稱方程求根：熟悉使用、迭代法、牛頓法、割線法等方法對給定的方程進行根的求解。選擇上述方法中的兩種方法求方程：二分法f(x)=x3+4x2-10=0在1,2內的一個實根，且要求滿足精度|x*-xn|0.510-5迭代法：

2、用迭代公式x=f(x)進行迭代計算，直到滿足|x*-xn|0.510-5 為止 。二分法：設f（x）在a,b上連續(xù)，且f（a1）*f（x1）0,記（a2，b2）=（x1，b1）帶入計算式進行計算 直到 |x*-xn|0.510-5 為止 。二、目的和意義（1）了解非線性方程求根的常見方法，如二分法、迭代法、牛頓法、割線法。（2）加深對方程求根方法的認識，掌握算法。會進行誤差分析，并能對不同方法進行比較。三、計算公式（1）迭代法 1）.首先對給定的計算公式進行變形使其能夠迭代或者找出相應迭代速度較快的式子。 2）.帶入求好的式子到循環(huán)中去比如：（2）二分法：f(x)在區(qū)間（x，y）上連續(xù) 1）.

3、先找到a、b屬于區(qū)間（x，y），使f(a)，f(b)異號，說明在區(qū)間(a,b)內一定有零點，然后求f(a+b)/2, 2）.如果f(a+b)/2=0，該點就是零點， 如果f(a+b)/20,則在區(qū)間（(a+b)/2，b)內有零點，反之在（a，(a+b)/2）內有零點 帶入1）中繼續(xù)。四、主要儀器設備Vc+ 9.0 C-free CodeBlocks五、結構程序設計迭代法: #include #include main() int i; double xn15,y,x1,x2,m ; printf(請輸入x1,x2的值：n ); scanf(%lf%lf,&x1,&x2); printf(請輸入

4、精度要求：n ); scanf(%lf,&m); printf( n xnn);i=0; do xn0=(x1+x2)/2 ; xni+1= sqrt(10/(4+xni); /迭代printf(%5d %5lfn,i,xni); y= fabs(xni+1-xni) ; i+; if(ym)break;while(1); 二分法：#include #include main() int m,n,o,p; double a,b,l; printf(請輸入x3, x2, x的系數和常數p：n);scanf(%d%d%d%d,&m,&n,&o,&p); /1 4 0 -10 printf(請輸入x

5、1,x2:n);scanf(%lf%lf,&a,&b); /1 2 printf(請輸入精度要求:n);scanf(%lf,&l); /0.5x105 printf( n an bn xn f(xn)n); double x,fx; int i=1; do x=(b+a)/2; fx=m*x*x*x+n*x*x+o*x+p; printf(%5d %5f %5f %5f %5fn,i,a,b,x,fx); i+; if(fx=0) break; if(fx0) b=x; else if(fx0) a=x; if(b-a)l) break; /進行計算并返值 while(1);六、結果討論和分析

6、 二分法： 迭代法： 分析討論：使用不同的方法，可以不同程度的求得方程的解，不同的方法速度不同，求得的結果也稍有區(qū)別，當然和要求精度也有關系。剛開始的時候用數組對二分法進行求解，發(fā)現 循環(huán)到第二次 就無法實現值的傳遞，于是換了另外一種方法代替了數組。流程圖：實驗地點綜合樓五層506室指導教師王崢太原理工大學學生實驗報告學院名稱計算機科學與技術專業(yè)班級計科 1002學號2010001414學生姓名 xxx實驗日期2012.6.5成績課程名稱計算數值方法實驗題目實驗二 線性方程組的直接解法一、課題名稱線性方程組的直接解法合理利用Gauss消元法、LU分解法、追趕法求解下列方程組： （n=5,10,

7、100）二、目的和意義（1）了解線性方程組常見的直接解法，如Guass消元法、LU分解法、追趕法。（2）加深對線性方程組求解方法的認識，掌握算法。（3）會進行誤差分析，并能對不同方法進行比較。三、計算公式 高斯分解法：將原方程組化為三角形方陣的方程組：lik=aik/akk aij= aij- lik* akj k=1,2,n-1 i=k+1,k+2, ,n j=k+1,k+2, ,n+1由回代過程求得原方程組的解： xn= ann+1/ ann xk=( akn+1-akj xj)/ akk (k=n-1,n-2, ,2,1) LU分解法：將系數矩陣A轉化為A=L*U, L為單位下三角矩陣，

8、U為普通上三角矩陣，然后通過解方程組l*y=b,u*x=y,來求解x.追趕法:用來求對角方程組；將系數矩陣A轉化為A=L*U, L為普通下n-1對角矩陣，U為單位上n-1對角矩陣，然后通過解方程組l*y=b,u*x=y,來求解x.四、主要儀器設備Vc+ 9.0 C-free CodeBlocks五、結構程序設計 Gauss消元法： #include #include using namespace std;int main()int n,i,j,k;double a100100,b100,o;cout輸入未知數個數:n;cout輸入數列:endl; for (i=1;i=n;i+)for (j

9、=1;jaij;for (i=1;i=n;i+)for (j=i+1;j1e-7)o=aii/aji;for (k=i;k0;i-)bi=ain+1/aii;for (j=i-1;j0;j-)ajn+1=ajn+1-bi*aji;cout解得:endl;for (i=1;i=n;i+) coutbiendl;/system(pause);return 0;列主元素消元法： #include #include #define N 20 using namespace std; void load(); float aNN; int m; int main() int i,j; int c,k,n

10、,p,r; float xN,lNN,s,d; coutm; coutendl; cout請按順序輸入增廣矩陣a:endl; load(); for(i=0;im;i+) for(j=i;jfabs(aii)?j:i; /*找列最大元素*/ for(n=0;nm+1;n+) s=ain; ain=acn; acn=s; /*將列最大數防在對角線上*/ for(p=0;pm+1;p+) coutaipt; coutendl; for(k=i+1;km;k+) lki=aki/aii; for(r=i;r=0;i-) d=0; for(j=i+1;jm;j+) d=d+aij*xj; xi=(ai

11、m-d)/aii; /*求解*/ cout該方程組的解為:endl; for(i=0;im;i+) coutxi=xit; return 0; void load() int i,j; for(i=0;im;i+) for(j=0;jaij; LU分解法：#include void solve(float l100,float u100,float b,float x,int n) int i,j; float t,s1,s2; float y100; for(i=1;i=n;i+) /* 第一次回代過程開始 */ s1=0; for(j=1;j=1;i-) /* 第二次回代過程開始*/ s2

12、=0; for(j=n;ji;j-) t=-uij; s2=s2+t*xj; xi=(yi+s2)/uii; void main() float a100100,l100100,u100100,x100,b100; int i,j,n,r,k; float s1,s2; for(i=1;i=99;i+)/*將所有的數組置零，同時將L矩陣的對角值設為1*/ for(j=1;j=99;j+) lij=0,uij=0; if(j=i) lij=1; printf (輸入方程組的個數 n:n);/*輸入方程組的個數*/ scanf(%d,&n); printf (讀取原矩陣 A（x的系數）:n);/*

13、讀取原矩陣A*/ for(i=1;i=n;i+) for(j=1;j=n;j+) scanf(%f,&aij); printf (讀取列矩陣 B（y的值）:n);/*讀取列矩陣B*/ for(i=1;i=n;i+) scanf(%f,&bi); for(r=1;r=n;r+)/*求解矩陣L和U*/ for(i=r;i=n;i+) s1=0; for(k=1;k=r-1;k+) s1=s1+lrk*uki; uri=ari-s1; for(i=r+1;i=n;i+) s2=0; for(k=1;k=r-1;k+) s2=s2+lik*ukr; lir=(air-s2)/urr; printf(輸

14、出矩陣 L:n);/輸出矩陣Lfor(i=1;i=n;i+) for(j=1;j=n;j+) printf(%7.3f ,lij); printf(n); printf(輸出矩陣 U:n);/輸出矩陣U for(i=1;i=n;i+) for(j=1;j=n;j+) printf(%7.3f ,uij); printf(n); solve(l,u,b,x,n); printf(解為:n); for(i=1;i=n;i+) printf(x%d=%fn,i,xi); 追趕法： #include #define N 3 main() double A33,b3; printf(請按順序輸入x的系數

15、:n); int a,c; for(a=0;a3;a+) for(c=0;c3;c+) scanf(%lf,&Aac); printf(請按順序輸入y的值:n); int k; for(k=0;k3;k+)scanf(%lf,&bk); int i; A01=A01/A00; for(i=1;i2;i+) Aii+1=Aii+1/(Aii-Aii-1*Ai-1i); for(i=1;i3;i+) Aii=Aii-Aii-1*Ai-1i; b0=b0/A00; for(i=1;i=0;i-) bi=bi-Aii+1*bi+1; for(i=0;i3;i+) printf(x%d=%.6lfn,i

16、,di); 六、結果討論和分析Gauss消元法： 列主元素消元法： LU分解法： 追趕法： 分析討論從消元過程可以看出，對于n階線性方程組，只要各步主元素不為零，經過n-1步消元，就可以得到一個等價的系數矩陣為上三角形陣的方程組，然后再利用回代過程可求得原方程組的解. 由于列主元素法相似且優(yōu)于完全主元素法 所以省略了后者。消元過程相當于分解 A為單位下三角陣L與上三角陣U的乘積，解方程組Ly=b回代過程就是解方程組Ux=y。其中的L為n階單位下三角陣、U為上三角陣. 在 A 的LU 分解中, L取下三角陣, U 取單位上三角陣,這樣求解方程組Ax=d 的方法稱為追趕法。另外是追趕法和其他方法求

17、同一方程結果不一樣，我多次修改源程序，也不知道原因。再就是追趕法有很大的局限性 還待改良。流程圖：實驗地點綜合樓五層506室指導教師王崢太原理工大學學生實驗報告學院名稱計算機科學與技術專業(yè)班級計科 1002學號2010001414學生姓名 xxx實驗日期2012.6.5成績課程名稱計算數值方法實驗題目實驗三 線性方程組的迭代解法一、課題名稱線性方程組的迭代解法使用雅可比迭代法或高斯-賽德爾迭代法對下列方程組進行求解。 二、目的和意義學習使用雅可比迭代法或高斯-賽德爾迭代法三、計算公式雅克比迭代法：設線性方程組Ax=b的系數矩陣A可逆且主對角元素a11,a22,ann均不為零，令D=diag(a

18、11,a22,ann)并將A分解成A=(A-D)+D從而線性方程組可寫成Dx=(D-A)x+b則有迭代公式x(k+1)=B1x(k)+f1其中，B1=I-D-1A,f1=D-1b。四、主要儀器設備 Vc+ 9.0 C-free CodeBlocks五、結構程序設計雅克比迭代法：#include #include main() int i; double x120 ,x220,x320; double x10, x20, x30; printf(請輸入x1,x2,x3的初值：n);scanf(%lf%lf%lf,&x10,&x20, &x30); printf( n x1n x2n x3n n)

19、; for(i=0;i18;i+) x10=x10; x20=x20; x30=x30; x1i+1=0.1*x2i+0.2*x3i+0.72;x2i+1=0.1*x1i+0.2*x3i+0.83;x3i+1=0.2*x1i+0.2*x2i+0.84;printf(%5d %5lf %5lf %5lfn,i,x1i,x2i,x3i); 六、實驗結果與分析：雅克比迭代法：分析討論： 其實，這兩個迭代法是之前迭代法的升級，多了幾個迭代式子而已，而且兩者相差不大比較簡單，所以選擇了雅克比迭代法進行求解，但是沒有與另一種方法 高斯賽德爾迭代法進行實質性的比較。流程圖：實驗地點綜合樓五層506室指導教師

20、王崢太原理工大學學生實驗報告學院名稱計算機科學與技術專業(yè)班級計科 1002學號2010001414學生姓名 xxx實驗日期2012.6.10成績課程名稱計算數值方法實驗題目實驗四 矩陣特征值與特征向量問題一、課題名稱使用冪法求A模為最大的特征值及其相應的特征向量。二、目的和意義（1）了解矩陣特征值與特征向量問題解法，掌握冪法。（2）加深對矩陣特征值與特征向量問題求解方法的認識，掌握算法。三、計算公式冪法：由已知的非零向量x0和矩陣A的乘冪構造向量序列xn以計算矩陣A的按模最大特征值及其特征向量的方法，稱為冪法。迭代公式：結果可取 四、主要儀器設備Vc+ 9.0 C-free CodeBlock

21、s五、結構程序設計源代碼： #include#include#define N 3#define eps 1e-6#define KM 30float MaxValue(float x,int n) float Max=x0; int i; for (i=1;ifabs(Max)Max=xi; return Max; void PowerMethod(float *A) float UN,VN,r1,r2,temp; int i,j,k=0; while(kKM) k+; for(i=0;iN;i+) temp=0; for(j=0;jN;j+)temp+=*(A+i*N+j)*Uj; Vi=

22、temp; for(i=0;iN;i+)Ui=Vi/MaxValue(V,N); if(k=1)r1=MaxValue(V,N); else r2=MaxValue(V,N); if(fabs(r2-r1)eps)break; r1=r2; printf(r=%fn,r2); for(i=0;iN;i+)printf(y%d=%fn,i+1,Ui);void main() float ANN=2,-1,0,-1,2,-1,0,-1,2 ;float UN; /A的值 U0=1; U1=1; U2=1;/x0的值 PowerMethod(A0);六、結果討論和分析分析討論由于該程序將A矩陣和x0

23、的值編寫在程序中，所以要想修改成其他矩陣時比較麻煩，所以也有一定的局限性。 冪法是一種求任意矩陣A的按模最大特征值及其對應特征向量的迭代算法。該方法的最大優(yōu)點是計算簡單，容易在計算機上實現，對稀疏矩陣較為適合，但有時收斂速度很慢。流程圖：實驗地點綜合樓五層506室指導教師王崢太原理工大學學生實驗報告學院名稱計算機科學與技術專業(yè)班級計科 1002學號2010001414學生姓名 xxx實驗日期2012.6.10成績課程名稱計算數值方法實驗題目 實驗五 代數插值一、課題名稱（使用拉格朗日插值法或牛頓插值法求解：已知f(x)在6個點的函數值如下表所示，運用插值方法，求f(0.596)的近似值。x0.

24、400.550.650.800.901.05f(x)0.410750.578150.696750.888111.026521.25386二、目的和意義學習使用拉格朗日插值法或牛頓插值法求解三、計算公式設函數在區(qū)間a,b上n+1互異節(jié)點x0,x1,xn上的函數值分別為y0,y1,yn,求n次插值多項式Pn(x),滿足條件Pn(xj)=yj, j=0,1,n令Ln(x)=y0l0(x)+y1l1(x)+ynln(x)= yili(x)其中l(wèi)0(x),l1(x), ln(x) 為以x0,x1,xn為節(jié)點的n次插值基函數，則Ln(x)是一次數不超過n的多項式，且滿足Ln(xj)=yj, L=0,1,n

25、再由插值多項式的唯一性，得Pn(x)Ln(x)四、主要儀器設備Vc+ 9.0 C-free CodeBlocks五、結構程序設計#include #include main() int I;char L; double M100100; double x100,y100; double X=1,xx=0,w=1,N=0,P,R=1; int n=5; /coutn; /for(int i=0;i=n;i+) / /*cout請輸入xi的值：xi; cout請輸入yi的值：yi; Mi0=xi; Mi1=yi; */用二維保存所有數據 M00=0.40;M01=0.41075;M10=0.55;

26、M11=0.57815;M20=0.65;M21=0.69675;M30=0.80;M31=0.88811;M40=0.90;M41=1.02652;M50=1.05;M51=1.25386;for( int j=2;j=n+1;j+) for(int i=1;i=n;i+) Mij=(Mij-1-Mi-1j-1)/(Mi0-Mi-j+10); for(int i=1;i=n;i+) cout其i階均差為：Mii+1endl; coutxx; for(int i=0;in;i+) X*=xx-Mi0; N+=Mi+1i+2*X; P=M01+N; cout其函數值：y=Pendl; / 六、結

27、果討論和分析分析討論拉格朗日插值的優(yōu)點是插值多項式特別容易建立，缺點是增加節(jié)點是原有多項式不能利用，必須重新建立，即所有基函數都要重新計算，這就造成計算量的浪費。所以該程序 選擇了牛頓法，由于在輸入數據的時候比較麻煩 所以將數據存在了程序里面，只能在程序里修改數據。 另外由于精度問題計算結果有一定的誤差。流程圖：實驗地點綜合樓五層506室指導教師王崢太原理工大學學生實驗報告學院名稱計算機科學與技術專業(yè)班級計科 1002學號2010001414學生姓名 xxx實驗日期2012.6.10成績課程名稱計算數值方法實驗題目 實驗六 最小二乘法擬合多項式一、課題名稱給定數據點（xi ,yi），用最小二乘

28、法擬合數據的多項式，并求平方誤差。xi00.50.60.70.80.91.0yi11.751.962.192.442.713.00二、目的和意義1熟練運用已學計算方法求解方程組2加深對計算方法技巧，選擇正確的計算方法來求解各種方程組3培養(yǎng)使用電子計算機進行科學計算和解決問題的能力三、計算公式建立正規(guī)方程組：（xij+k）ak=xijyi ,j=0,1,n 平方誤差：I=（akxik-yi）2四、主要儀器設備Vc+ 9.0 C-free CodeBlocks五、結構程序設計源代碼： #include#include#define N 15double power(double &a,int n)

29、double b=1;for(int i=0;in;i+)b*=a;return b;void Gauss();double XN,YN,sumXN,sumYN,aNN,bN,lNN,xN;void main()ofstream outdata;ifstream indata;double s;int i,j,k,n,index;/coutn;n=7; coutendl;/cout請輸入X和Y:endl; /輸入給定數據X0=0.0;Y0=1.00;X1=0.5;Y1=1.75;X2=0.6;Y2=1.96;X3=0.7;Y3=2.19;X4=0.8;Y4=2.44;X5=0.9;Y5=2.71;X6=1.0;Y6=3.00; /綁定數據 /可以解綁由下for循環(huán) 輸入任何數據for(i=0;in;i+)/coutXiXi;sumX1+=Xi;/coutYiYi;sumY1+=Yi;/coutendl;coutsumX1=