數(shù)值代數(shù)上機報告

1、把矩陣A分解成一個下三角陣與一個上三角陣的乘積，即A=LR，其中L為下三角陣，R為上三角陣，這樣原線性方程組就可以化為Ly b的求解問題，方Rx y便求解。二、算法：1)輸入系數(shù)矩陣2)11 ij rij 和（l ji1行，計算得出矩陣L和R；利用公式rijaijajiljkki)/rii 交錯進k 1Ly b3)回帶到Ly b中得出原線性方程組的解。Rx y三、源程序：#i nclude <stdio.h>#i nclude <stri ng.h>#in clude <math.h>#i nclude <stdlib.h>#define N 1

2、00mai n()int i,j,k,s,n;printf("請輸入系數(shù)矩陣A的階數(shù)：n=");scan f("%d",&n);floataNN=0, LNN=0,RNN=0,sigma1,sigma2,bN,yN,xN;/*為L主對角線元素賦1*/for(i=0;i <n ;i+)Lii=1;printf("請輸入系數(shù)矩陣A: n");/*輸入系數(shù)矩陣A*/ for(i=0;i<n;i+)for(j=0;j<n;j+)scanf("%f",&aij);for(k=0;k<n

3、;k+)for(j=k;j<n;j+)/* 計算矩陣 R*/sigma1=0;for(s=0;s<=k-1;s+)sigma1+=Lks*Rsj;Rkj=akj-sigma1;for(i=k;i<n;i+)/* 計算矩陣 L*/sigma2=0;for(s=0;s<=k-1;s+)sigma2+=Lis*Rsk;Lik=(aik-sigma2)/Rkk;printf("n A 矩陣為： n");/* 輸出矩陣 L、 R*/ for(i=0;i<n;i+)for(j=0;j<n;j+)printf("%5.1f ",ai

4、j);printf("n");printf("n L 矩陣為： n");for(i=0;i<n;i+)for(j=0;j<n;j+)printf("%5.1f ",Lij);printf("n");printf("n R 矩陣為： n");for(i=0;i<n;i+)for(j=0;j<n;j+)printf("%5.1f ",Rij);printf("n");printf("請輸入 b 矩陣n",i+1);f

5、or(i=0;i <n ;i+)scan f("%f",&bi);for(i=0;i<n;i+)/*回代法求解方程組Ly=b*/ sigma1=0;for(k=0;k<=i-1;k+) sigma1+=Lik*yk; yi=bi-sigma1;for(i=n-1;i>=0;i-)sigma2=0;for(k=i+1;k <n; k+) sigma2+=Rik*xk;xi=(yi-sigma2)/Rii;printf("解得 x 為：n");for(i=0;i <n ;i+)prin tf("%5.1f

6、 ",xi); prin tf("n");A矩陣為：2.01.0l.B1.03.02.01.02.02.0四、計算結(jié)果與分析:0 0 0 0 0 1 0 9 & » « 0 1 U 0 _b ra * - 10 6H矩陣為h2.0l.,06-92.50.0 01.61_0 1.0 Fi*e©s an u kel.Sto continue分析：運行結(jié)果與預(yù)想的結(jié)果相近，誤差對結(jié)果的影響不是很大，比較理想五、參考文獻：1 刑志棟矩陣數(shù)值分析陜西：陜西科學(xué)技術(shù)出版社，20052 譚浩強.C語言程序設(shè)計北京：清華大學(xué)出版社，2005伊

7、報告2cholesky 分解目的意義:對稱正定矩陣是實踐中經(jīng)常遇到的一種特殊矩陣類型矩陣，由于矩陣本身的 流量好興致，使得cholesky分解在存儲和運算量上較一般消去法節(jié)省一半左右， 且解的精度高，cholesky分解方法是目前計算機求解該類問題最有效的方法之、算法:1)輸入系數(shù)矩陣A ；Sijaij1Sik 1 jk12)3)lijdiikSij / d jjii 1aiisk l ikk 1交錯進行，計算得出矩陣L和D；回帶到J?bDTy中得出原線性方程組的解三、源程序：#in clude <stdio.h>#in clude <stdlib.h>#in clud

8、e <math.h>#define EPS 1.0e-8#define N 20 double aNN, bN, xN;int n;int zhuyuan(int row);/* 選主元 */void hangjiaohuan(int row1, int row2); /* 行交換 */ void xiaoyuan(int row);/* 消元 */void huidai(); /* 回代 */ void main()printf(" 請輸入方程的維數(shù) n: n = ");scanf("%d", &n);getchar();print

9、f(" 輸入 %d 行 %d 列系數(shù)矩陣 A：n", n, n); for (int i=0; i<n; i+)for (int j=0; j<n; j+)scanf("%lf", &aij);getchar();printf("n 輸入線性方程組右端項 b%d:n", n);for (i=0; i<n; i+)scanf("%lf", &bi);getchar();for (i=0; i<n-1; i+)double rowmark = zhuyuan(i);if (row

10、mark = -1)printf(" 多解 !"); system("pause");return;if (rowmark != i)hangjiaohuan(i, rowmark);xiaoyuan(i);huidai();printf("n 線性方程組的解為 :n ");for (i=0; i<n; i+)printf("x%d=%lf ", i+1, xi);printf("n"); printf("n"); system("pause");i

11、nt zhuyuan(int row)double elem = arowrow; int rowmark = row; for (int i=row+1; i<n; i+)if (elem<airow)elem = airow; rowmark = i; if(fabs(elem) <= EPS) return -1;return rowmark;void hangjiaohuan(int row1, int row2)double tmp; tmp = brow1; brow1 = brow2; brow2 = tmp; for (int j=0; j<n; j+)

12、tmp = arow1j; arow1j = arow2j; arow2j = tmp;void xiaoyuan(int row)for (int i=row+1; i<n; i+)int j=row; double tmp = aij/arowrow;bi -= tmp*brow;for (aij+ = 0; j< n; j+)aij -= tmp*arowj;void huidai()x n-1 = bn-1/a n-1 n-1;for (int i=n_2; i>=0; i-)double sum = 0.0;for (int j=i+1; j< n; j+)

13、sum -= aij*xj;xi = (bi+sum)/aii;四、計算結(jié)果與分析:分析：運行結(jié)果與預(yù)想的結(jié)果相近，誤差對結(jié)果的影響不是很大，比較理想五、參考文獻：1刑志棟矩陣數(shù)值分析陜西：陜西科學(xué)技術(shù)出版社，20052譚浩強.C語言程序設(shè)計北京：清華大學(xué)出版社，2005目的意義:Jacobi迭代法通過有限次的算術(shù)運算得到線性方程組的近似值、算法：1）輸入系數(shù)矩陣A ；2）輸入迭代的初始向量;k 13）利用公式人binkaijXjj1jiai，計算得出原線性方程組的近似解三、源程序:#in clude<stdio.h>#in clude<c oni o.h>#i ncl

14、ude<malloc.h>#in clude<math.h>#define EPS 1e-6#defi ne MAX 100#defi ne N 100float *Jacobi(float aNN,int n)float *x,*y,s;double epsil on;int i,j,k=0;x=(float *)malloc( n*sizeof(float); y=(float *)malloc( n*sizeof(float); printf("請輸入迭代的初始向量：n"); for(i=0;i< n;i+) sca nf("%

15、f", &xi);while(1)epsilon=0;k+;for(i=0;i<n;i+)s=0;for(j=0;j<n;j+)if(j=i) continue; s+=aij*xj;yi=(ain-s)/aii;epsilon += fabs(yi-xi);if(epsilon < EPS)printf("n 迭代次數(shù)為： %d 次 n",k); return y;if(k>=MAX)printf("The Method is disconvergent"); return y;for(i=0;i<n;i

16、+)xi=yi;main()int i,j,n;float aNN;float *x;printf(" 請輸入系數(shù)矩陣的階數(shù)： n="); scanf("%d",&n);printf(" 請按行輸入線性方程組的增廣矩陣： n"); for(i=0;i<n;i+)for(j=0;j<n+1;j+)scanf("%f",&aij);x=(float *)malloc(3*sizeof(float); x=Jacobi(a,n);printf("n 原線性方程組的解為： n"

17、;);for(i=0;i< n;i+) prin tf("x%d=%fn",i,xi);getch();四、計算結(jié)果與分析：幕I篙賀蠶黠螫滿卜矩陣，12-2111132 2 15請輸入遞代的初始向量匕回E) 0迭代次數(shù)為|嵌原線性方程組的解為：x0J=l.000000X11>1.000000Fr*uw百 dny key toinueH分析：運行結(jié)果與預(yù)想的結(jié)果相近，誤差對結(jié)果的影響不是很大，比較理想五、參考文獻：1 刑志棟矩陣數(shù)值分析陜西：陜西科學(xué)技術(shù)出版社，20052 譚浩強.C語言程序設(shè)計北京：清華大學(xué)出版社，2005目的意義:階數(shù)超過四次的多項式的根一般不

18、能用有限次運算求得，因而矩陣特征值的計算方法本質(zhì)上總是采用迭代格式，而采用乘幕法可以計算最大的一個或按模最 大的幾個特征值和特征相應(yīng)特征向量。、算法:1)輸入系數(shù)矩陣A ；2）取出事向量為zo 1,1,1* AZk i；3) 利用乘幕法的迭代格式J mk |yj ；計算得出原矩陣按模最大特征值和厶 yk /mk相應(yīng)特征的向量。三、源程序：#i nclude<stdio.h>#in clude<math.h>#define N 100#define M 100#define epsilon 0.00001int mai n()int n;int i,j,k;double

19、xmax,oxmax;static double aNN;static double xN, nxN;printf("輸入矩陣維數(shù)n:");sca nf("%d",&n);if(n >N)prin tf("the in pur n is larger tha n MAX_N,please redefi ne the N.n"); return 1;if(n<=0)prin tf("please in put a nu mber betwee n 1 and %d. n",N); return 1

20、;輸入A矩陣printf("請輸入矩陣 A(%d,%d):n",n,n);for(i=0;i< n;i+)for(j=0;j <n ;j+)sca nf("%lf",&aij);for(i=0;i< n;i+)xi=1;oxmax=0;for(i=0;i<M;i+)for(j=0;j <n ;j+)nxj=0;for(k=0;k< n;k+)n xj+=ajk*xk;xmax=0.0;for(j=0;j <n ;j+)if(fabs( nxj)>xmax) xmax=fabs( nxj);for(j

21、=0;j <n ;j+)n xj/=xmax;for(j=0;j <n ;j+)xj=n xj;if(fabs(xmax-oxmax)<epsilo n)printf("按模最大特征值=%.4fn",xmax); printf("相應(yīng)的特征向量為:n"); for(i=0;i <n ;i+)prin tf("%.4fn", nxi);break;/return 0;oxmax=xmax;return 0;四、計算結(jié)果與分析：2 3 210 3 43 6 1相應(yīng)的特彳正向量為;1 0060Pi'ess a

22、ny Ice to cant in u.e分析：運行結(jié)果與預(yù)想的結(jié)果相近，誤差對結(jié)果的影響不是很大，比較理想五、參考文獻：1 刑志棟矩陣數(shù)值分析陜西：陜西科學(xué)技術(shù)出版社，20052 譚浩強.C語言程序設(shè)計北京：清華大學(xué)出版社，2005 報告5逆幕法、目的意義：逆幕法用于計算非奇異矩陣A的按模最小特征值和相應(yīng)的特征向量、算法：1）輸入系數(shù)矩陣A和近似特征值；2）將乘幕法用于A 1便可求出A的按模最小的特征值和相應(yīng)的特征向量。三、源程序：#in clude<stdio.h>#in clude<math.h>#define e 0.0001#define N 100void

23、mai n()int i,j,k, n,p,s,kk=1,flag=1;double q,m, lm,L m,m0;double ANN,BNN,LNN,RNN,lNN,rNN;double z0N,z1N,y0N,y1N,yN;printf("請輸入矩陣的階數(shù)：n=");scan f("%d",&n);printf("n 請輸入矩陣 A:n");for(i=1;i<=n ;i+)for(j=1;j<=n ;j+)sca nf("%lf",&Aij);printf("n請輸入近

24、似特征值：");scanf("%lf",&lm);/* 令 B=A-lm*I , 用 Doolittle 分解 B=LR */ for(i=1;i<=n;i+) for(j=1;j<=n;j+)if(j=i)Bij=Aij-lm;elseBij=Aij;for(k=1;k<=n;k+)for(i=k;i<=n;i+)for(q=0,p=1;p<k;p+) q=q+Lkp*Rpi;Rki=Bki-q;for(i=k+1;i<=n;i+)for(q=0,p=1;p<k;p+) q=q+Lip*Rpk;Lik=(Bik-

25、q)/Rkk;for(i=1;i<=n;i+)for(j=i;j<=n;j+)if(i=j)Lij=1;elseLij=0;for(i=2;i<=n;i+)for(j=1;j<i;j+)Rij=0;/* 第一步迭代，用 Doolittle 計算 R*y1=y0 中的 y1*/ for(i=1;i<=n;i+)y0i=1;for(k=1;k<=n;k+)for(i=k;i<=n;i+)for(q=0,p=1;p<k;p+)q=q+lkp*rpi; rki=Rki-q;for(i=k+1;i<=n;i+)for(q=0,p=1;p<k;p

26、+) q=q+lip*rpk;lik=(Rik-q)/rkk;for(i=1;i<=n;i+)for(q=0,k=1;k<i;k+) q=q+lik*yk;yi=y0i-q;for(i=n;i>=1;i-)for(q=0,k=i+1;k<=n;k+)q=q+rik*y1k; y1i=(yi-q)/rii;/*計算m,z1(存放在下列代碼的z0中)*/ for(m=y11,i=2;i<=n;i+)if(fabs(y1i)>fabs(y1i-1) m=y1i;for(i=1;i<=n;i+) z0i=y1i/m;/*第n (>1)次迭代*/ kk=2

27、;while(flag=1)/*用 Doolittle 解方程 LR*y2=z1 中的 y2,其中 B=LR */ for(k=1;k<=n;k+)for(i=k;i<=n;i+)for(q=0,p=1;p<k;p+)q=q+lkp*rpi;rki=Bki-q;for(i=k+1;i<=n;i+)for(q=0,p=1;p<k;p+)q=q+lip*rpk;lik=(Bik-q)/rkk;for(i=1;i<=n;i+) for(q=0,k=1;k<i;k+)q=q+lik*yk;yi=z0i-q;for(i=n;i>=1;i-) for(q=0

28、,k=i+1;k<=n;k+)q=q+rik*y1k;y1i=(yi-q)/rii;/*規(guī)范化 */for(m=y11,i=2;i<=n;i+) if(fabs(y1i)>fabs(y1i-1)m=y1i;for(i=1;i<=n;i+)；否則得解。 */z1i=y1i/m;/*判斷：若|Zk-Zk-1|>=e，則進行下一輪迭代for(s=0,i=1;i<=n;i+)if(fabs(z0i)-fabs(z1i)<e)s+; if(s=n)flag=0;else kk+;for(i=1;i<=n;i+) y0i=y1i;z0i=z1i;m0=1/m

29、;Lm=lm+mO;printf("最接近的特征值為：n"); prin tf("%fnn",Lm);printf("相應(yīng)的特征向量為：n"); for(i=1;i<=n ;i+)prin tf("%lf,",z1i);prin tf("nn");四、計算結(jié)果與分析：分析：運行結(jié)果與預(yù)想的結(jié)果相近，誤差對結(jié)果的影響不是很大，比較理想 五、參考文獻：1 刑志棟矩陣數(shù)值分析陜西：陜西科學(xué)技術(shù)出版社，20052 譚浩強.C語言程序設(shè)計北京：清華大學(xué)出版社，2005報告6古典Jacobi算法目的意

30、義:尋找相似變換，是對稱矩陣 A經(jīng)過變換之后所得矩陣的非對角線元素的平 方和減少，對角線元素的平方和增大，且保持對稱性不變，不斷的施行這種正交變換，最終是非對角元素的平方和任意的接近與零， 對角線元素平方的取極大值。、算法:1）輸入系數(shù)矩陣A ；2）取AA，形成一個相似矩陣序列 人k3）選取Rk，使得元素apq =0，這時 應(yīng)滿足4）令 y appaqq， x sgn（aPp ）RkAk iRk， k=1,2,tan2 點敢a ppaqq 1)*2apVk 1aqq因此可以得到：cos21/2（1 譏 _ ） ，sin 1/2* x/、, _ * 1/cos三、源程序：#i nclude &q

31、uot;stdio.h"#i nclude"math.h"#define P 0.1void mai n()double a202020,b202020,r202020,l202020;double max,u,y,x,t,g,f,h,v,w;int i,j,k, n,c,d;printf("請輸入矩陣的維數(shù)：n=");sca nf("%d",&n);printf("請輸入矩陣A:n");for(i=1;i<=n ;i+)for(j=1;j<=n ;j+)sca nf("%l

32、f",&aij1);printf("n所有的特征值為：");for(i=1;i<=n ;i+)for(j=1;j<=n ;j+)if(aij1<0)bij1=-aij1;else bij1=aij1;for(k=1;k<=n ;k+)for( max=0,i=1;i<=n ;i+)for(j=1;j<=n ;j+)if(i<j&&bijk>max+0.000001)max=bijk;c=i;d=j;u=bcck-bddk;if(u>=0) y=u;else y=-u;if(u>=0

33、) x=2*bcdk;else x=-2*bcdk;t=sqrt(x*x+y*y);g=sqrt(0.5*1+0.5*y/t);f=0.5*x/t*(1/g);for(i=1;i<=n;i+)for(j=1;j<=n;j+)if(i!=j&&i!=c,j!=d&&i!=d,j!=c)rijk=0;if(i=c&&j=c) rcck=g;if(i=d&&j=d) rddk=g;if(i=c&&j=d) rcdk=f;if(i=d&&j=c) rdck=-f; if(i=j&&

34、;i!=c&&j!=d) rijk=1;printf("n");for(i=1;i<=n;i+) for(j=1;j<=n;j+) aick+1=bick*g+bidk*f;acik+1=bick*g+bidk*f;aidk+1=-bick*f+bidk*g;adik+1=-bick*f+bidk*g; if(i!=c&&j!=d&&i!=d&&j!=c) aijk+1=bijk;h=(bcck-bddk)*f*g;v=bcdk*(g*g-f*f);acdk+1=h+v;adck+1=h+v;acck+1=bcck*g*g+2*bcdk*f*g+bddk*f*f;addk+1=bcck*f*f-2*bcdk*f*g+bddk*g*g;for(i=1;i<=n;i+)for(j=1;j<=n ;j+) if(aij1<0)bijk+1=-aijk+1; elsebijk+1=aijk+1;w=0;for(i=1;i<=n ;i+)for(j=1;j<=n ;j+)if(i!=j)w=w+bijk+1*bijk+1; elsecon ti nue;prin tf("