數(shù)值分析課程設計(求解線性方程組)_第1頁
已閱讀1頁,還剩24頁未讀, 繼續(xù)免費閱讀

下載本文檔

版權說明:本文檔由用戶提供并上傳,收益歸屬內容提供方,若內容存在侵權,請進行舉報或認領

文檔簡介

1、<p><b>  數(shù)值分析 課程設計</b></p><p><b>  作者姓名:</b></p><p><b>  學號: </b></p><p><b>  一、 問題的提出</b></p><p>  分別用SOR方法和高斯消元的L

2、U分解算法(lii=1, i=1,…,n)求解給定的線性方程組AX=B, 以感受迭代法和直接法的不同特點。</p><p><b>  二、 實驗內容</b></p><p>  自定義函數(shù) SOR(A, B, w, MAXN, TOL),以實現(xiàn)SOR方法求解線性方程組AX=B,其中</p><p><b>  A——系數(shù)矩陣;<

3、;/b></p><p><b>  B——常數(shù)列向量;</b></p><p><b>  w——松弛因子;</b></p><p>  MAXN——迭代的最大次數(shù)</p><p>  TOL——達到的精度上限</p><p>  返回值有以下四種可能:</p&

4、gt;<p>  -2:SOR方法不收斂;(不收斂的依據為的某個分量值超出區(qū)間[-108, 108]。)</p><p>  -1:矩陣有一列全為0;</p><p>  0:算法經過MAXN次迭代還未收斂;</p><p>  k:SOR方法經k次迭代收斂,求得方程組的解向量X記錄下來.</p><p>  自定義函數(shù)Dire

5、ct(A, B),以實現(xiàn)高斯LU分解的方法求解線性方程組AX=B,其中</p><p><b>  A——系數(shù)矩陣;</b></p><p><b>  B——常數(shù)列向量;</b></p><p><b>  返回值有兩種可能:</b></p><p>  “LU decomp

6、sition failed.”:分解過程中U的對角線元素至少一個為0;</p><p><b>  X:分解過程中</b></p><p>  分別使用步驟1中定義的函數(shù)SOR(A, B, w, MAXN, TOL)和步驟2中定義的函數(shù)Direct(A, B)進行測試,記錄返回值及X值(算法收斂或有效的情形, 保留4位小數(shù)):</p><p>

7、<b>  測試1:</b></p><p>  MAXN =1000,TOL =10-9,w分別取1, 1.05, 1.1, 1.2, 1.3, 1.6, 1.95;</p><p><b>  測試2:</b></p><p>  MAXN =1000,TOL =10-9,w=1;</p><p&g

8、t;<b>  測試3:</b></p><p>  MAXN =1000,TOL =10-9,w=1.2;</p><p><b>  測試4:</b></p><p>  MAXN =1000,TOL =10-9,w=1, 1.1, 1.3, 1.8;</p><p>  測試5:: n階Hil

9、bert矩陣定義為</p><p>  取n=3, MAXN =1000,TOL =10-9,w=1, 1.3, 1.6, 1.9;</p><p>  測試6:A為4階Hilbert矩陣,MAXN =10000,TOL =10-6,w=1, 1.3, 1.6, 1.8, 1.9.</p><p><b>  三、實驗結果及分析</b><

10、;/p><p><b> ?。ㄒ唬?SOR方法</b></p><p>  1. SOR法分析:</p><p>  (1)利用高斯SOR法可得迭代公式:</p><p>  X1(k+1)=(1-w)X1(k)-w/4(-X2(k)-X4(k))</p><p>  X2(k+1)=(1-w)X2(

11、k)-w/4(-X1(k+1)-X3(k)-X5(k)-5)</p><p>  X3(k+1)=(1-w)X3(k)-w/4(-X2(k+1)-X6(k))</p><p>  X4(k+1)=(1-w)X4(k)-w/4(-X1(k+1)-X5(k)-6)</p><p>  X5(k+1)=(1-w)X5(k)-w/4(-X2(k+1)-X4(k+1)-X(k

12、+1)+2)</p><p>  X6(k+1)=(1-w)X6(k)-w/4(-X3(k)-X5(k)-6);</p><p>  將松弛系數(shù)w的不同德值代入計算出X的值。</p><p>  利用迭代使|X(k+1)-X(k)|<e(e為精度,e=10^(-9))。</p><p>  (2)由于矩陣出現(xiàn)了,一列為0,所以不能使用迭

13、代,在程序中會出現(xiàn)r=-1.</p><p>  (3)X1(k+1)=(1-w)X1(k)-w/3(-X2(k)+3X3(k))</p><p>  X2(k+1)=(1-w)X2(k)-w/6(3X1(k+1)+3X3(k))</p><p>  X3(k+1)=(1-w)X3(k)-w/3(3X1(k+1)+3X2(k+1))</p><p

14、>  利用迭代使|X(k+1)-X(k)|<e(e為精度,e=10^(-9))。</p><p>  (4)將方程組變?yōu)椋?lt;/p><p>  X1+0X2+2X3+0X4+3X5+0X6+4X7=3</p><p>  3X1-1X2+0.5X3+8X4+2.2X5+1.6X6+0X7=8</p><p>  3X1+3X2+0

15、.5X3+12.5X4+5.4X5+3.6X6+X7=10</p><p>  5X1+2X2+5.5X3+8X4+2.2X5+1.6X6+3.3X7=12</p><p>  X1-4X2-1.5X3+9X4+2.2X5+1.6X6+3.3X7=9</p><p>  5.5X1+3.5X2+0.5X3+8X4+3.2X5+1.6X6+0X7=6</p>

16、;<p>  -0.5X1-1.5X2+3X3+2X4+0X5+X6-X7=5</p><p><b>  迭代公式為:</b></p><p>  X1(k+1)=(1-w)X1(k)-w(2X3(k)+3X5(k)+4X7(k)-3)</p><p>  X2(k+1)=(1-w)X2(k)+w(3X1(k+1)+0.5X3(

17、k)+8X4(k)+2.2X5(k)+1.6X6(k)+X7(k)-10) </p><p>  X3(k+1)=(1-w)X3(k)+w/0.5(3X1(k+1)+3X2(k+1)+12.5X4(k)+5.4X5(k)</p><p>  +3.6X6(k)+X7(k)-10)</p><p>  X4(k+1)=(1-w)X4(k)+w/8(5X1(k+1)+

18、2X2(k+1)+5.5X3(k+1)+2.2X5(k)</p><p>  +1.6X(k)+3.3X7(k)-12)</p><p>  X5(k+1)=(1-w)X5(k)+w/2.2(X1(k+1)-4X2(k+1)-1.5X3(k+1)+9X5(k+1)</p><p>  +1.6X6(k)+3.3X7(k)-9)</p><p>

19、;  X6(k+1)=(1-w)X6(k)+w/1.6(5.5X1(k+1)+3.5X2(k+1)+0.5X3(k+1)+8X4(k+1)+3.5X5(k+1)-6)</p><p>  X7(k+1)=(1-w)X7(k)-w(-0.5X1(k+1)-1.5X2(k+1)+3X3(k+1)+2X4(k+1)</p><p><b>  +X6(k)-6)</b>&l

20、t;/p><p>  將松弛系數(shù)w的不同德值代入計算出X的值。</p><p>  利用迭代使|X(k+1)-X(k)|<e(e為精度,e=10^(-9))。</p><p> ?。?) 將方程組變?yōu)椋?lt;/p><p>  X1(k+1)=(1-w)X1(k)-2w(1/3X2(k)+1/4X3(k)-1)</p><p

21、>  X2(k+1)=(1-w)X2(k)-4w(1/3X1(k+1)+1/5X3(k)-1))</p><p>  X3(k+1)=(1-w)X3(k)-6w(1/4X1(k+1)+1/5X2(k+1)-1)</p><p>  將松弛系數(shù)w的不同德值代入計算出X的值。</p><p>  利用迭代使|X(k+1)-X(k)|<e(e為精度,e=10^

22、(-9))。</p><p>  (6)將方程組變?yōu)椋?lt;/p><p>  X1(k+1)=(1-w)X1(k)-2w(1/3X2(k)+1/4X3(k)+1/5X4(k)-1)</p><p>  X2(k+1)=(1-w)X2(k)-4w(1/3X1(k+1)+1/5X3(k)+1/6X4(k)-1)</p><p>  X3(k+1)=

23、(1-w)X3(k)-6w(1/4X1(k+1)+1/5X2(k+1)+1/6X4(k)-1)</p><p>  X4(k+1)=(1-w)X4(k)-8w(1/5X1(k+1)+1/6X2(k+1)+1/7X3(k+1)-1)</p><p>  將松弛系數(shù)w的不同德值代入計算出X的值。</p><p>  利用迭代使|X(k+1)-X(k)|<e(e為精

24、度,e=10^(-6))。</p><p>  2. 測試SOR方法</p><p><b>  測試1</b></p><p> ?。?)W=1的時候:</p><p> ?。?)W=1.05的時候:</p><p>  (3)W=1.1的時候:</p><p> ?。?

25、)W=1.2的時候</p><p>  (5)W=1.3的時候</p><p> ?。?)W=1.6的時候</p><p>  (7)W=1.95的時候</p><p><b>  測試2</b></p><p><b>  測試3</b></p><p&

26、gt;<b>  測試4</b></p><p><b>  (1)W=1的時候</b></p><p> ?。?)W=1.1的時候</p><p> ?。?)W=1.3的時候</p><p> ?。?)W=1.8的時候</p><p><b>  測試5</

27、b></p><p><b>  (1)W=1的時候</b></p><p> ?。?)W=1.3的時候</p><p> ?。?)W=1.6的時候</p><p> ?。?)W=1.9的時候</p><p><b>  測試6</b></p><p

28、><b> ?。?)W=1的時候</b></p><p>  (2)W=1.3的時候</p><p> ?。?)W=1.6的時候</p><p> ?。?)W=1.8的時候</p><p> ?。?)W=1.9的時候</p><p> ?。ǘ?高斯LU方法</p><p

29、>  1. 高斯LU分析:</p><p><b> ?。?)L矩陣為: </b></p><p>  利用LD=B,算出D為:</p><p><b>  U矩陣為:</b></p><p>  利用UX=D,求出X的值。</p><p><b> ?。?)

30、L矩陣為:</b></p><p><b>  U矩陣:</b></p><p>  由于U矩陣對角線出現(xiàn)了0,所以出現(xiàn)了“LU decompsition failed.”</p><p><b> ?。?)L矩陣:</b></p><p><b>  U矩陣:</b&g

31、t;</p><p>  由于U矩陣對角線出現(xiàn)了0,所以出現(xiàn)了“LU decompsition failed.”</p><p> ?。?)將矩陣A變?yōu)椋?lt;/p><p>  而在將矩陣變?yōu)長與U時候出現(xiàn)了異常,在L與U矩陣中有異常值,且在U矩陣對角線上出現(xiàn)了0值,所以出現(xiàn)了“LU decompsition failed.”</p><p>

32、<b> ?。?)L矩陣:</b></p><p>  利用LD=B,可得到D為:</p><p><b>  U矩陣為:</b></p><p>  利用UX=D,求出X的值出來。 </p><p><b> ?。?)L矩陣:</b></p><p>

33、;<b> ?。ㄒ娤乱豁摚?lt;/b></p><p>  利用LD=B,可以求出D;</p><p><b>  U矩陣為:</b></p><p>  利用UX=D,可以求出X的值出來。</p><p>  2. 測試高斯LU方法</p><p><b>  測試1

34、</b></p><p><b>  測試2</b></p><p><b>  測試3</b></p><p><b>  測試4</b></p><p><b>  測試5</b></p><p><b>

35、  測試6</b></p><p>  四、 關于本設計的體會</p><p>  通過對高斯LU法和SOR法設計分析和測試以后,我發(fā)現(xiàn)兩種方法各有優(yōu)劣。高斯LU法得出的結果精度比較高,但是卻不適用于所有的方程,使用范圍相對較窄;而在使用SOR法時,雖然精度會稍微差一點,但是通過調整松弛度w,卻可能適用于更多的方程,適用范圍相對較寬。</p><p> 

36、 在做課程設計這個過程中,我發(fā)現(xiàn)自己還有很多很多知識沒有學好,在參考別人的例子的時候好像很簡單,但自己一上機操作寫程序的時候就出現(xiàn)問題。調試的時候系統(tǒng)總有報錯,還有很多警告,每修改一個變量,往往都要調試很久,有時候僅僅只是少了一個大括號,卻地花上近半個小時才能找到問題的瓶頸所在。此外,通過本次的課程設計,我重溫了許多C語言的知識。同時也發(fā)現(xiàn)了自己對C語言的掌握程度有所下降。其實,為了更加方便我今后的學習,我還是有必要對MATLAB進行學

37、習,不斷擴充我的知識。</p><p>  最后,雖然得到的結果未如理想,但我會繼續(xù)努力!謝謝李老師的指導!</p><p><b>  五、參考文獻</b></p><p>  【1】 《標準C語言基礎教程》(第四版) 電子工業(yè)出版社</p><p>  【2】 《數(shù)值分析》(第三版)

38、 北京理工大學出版社</p><p><b>  六、 附錄</b></p><p>  1. SOR方法(源代碼):</p><p>  #include <stdio.h></p><p>  #include<math.h></p><p>  #defi

39、ne N 20</p><p>  float get(float, float);</p><p>  void main ()</p><p><b>  {</b></p><p>  int n,Q,p=0;</p><p>  printf("請輸入系數(shù)矩陣的階數(shù):"

40、;);</p><p>  scanf("%d",&n);</p><p>  int i, j, k=0,r;</p><p>  float a[N][N], b[N];</p><p>  float x[N], y[N],z[N];</p><p>  float e,t,w;<

41、;/p><p>  printf("請輸入最大迭代次數(shù)MAXN=");</p><p>  scanf("%d",&Q);</p><p><b>  float v;</b></p><p>  printf("請輸入精度TOL=");</p>

42、<p>  scanf("%f",&v);</p><p>  printf("請輸入松弛系數(shù)w=");</p><p>  scanf("%f",&w);</p><p>  printf ("請輸入系數(shù)矩陣:\n ");</p><p

43、>  for (i = 0; i < n; i++)</p><p><b>  {</b></p><p><b>  p=0;</b></p><p>  for (j = 0; j < n; j++)</p><p><b>  {</b></p&

44、gt;<p>  printf("a[%d][%d]:",j+1,i+1);</p><p>  scanf ("%f", &a[j][i]);</p><p>  if(a[j][i]==0)</p><p><b>  p=p+1;</b></p><p>

45、;<b>  }</b></p><p><b>  if(p!=n)</b></p><p><b>  continue;</b></p><p><b>  else</b></p><p><b>  r=-1;</b><

46、;/p><p>  printf("\n矩陣有一列全為0,松弛法不能迭代,r=%d\n",r);</p><p><b>  exit(0);</b></p><p><b>  }</b></p><p>  printf ("請輸入右端項數(shù)組: \n");&l

47、t;/p><p>  for (i = 0; i < n; i++)</p><p><b>  {</b></p><p>  printf("b[%d]:",i+1);</p><p>  scanf ("%f", &b[i]);</p><p&g

48、t;<b>  }</b></p><p>  for (i = 0; i < n; i++)</p><p><b>  x[i]=0;</b></p><p><b>  do</b></p><p><b>  {</b></p>

49、<p>  for (i = 0; i < n; i++)</p><p><b>  {</b></p><p><b>  t =0.0 ;</b></p><p>  for (j = 0; j < n; j++)</p><p><b>  {</b&

50、gt;</p><p>  if (j < i)</p><p>  t += a[i][j] * y[j];</p><p>  else if (j > i)</p><p>  t += a[i][j] * x[j];</p><p><b>  }</b></p>

51、<p>  z[i] =(b[i] - t) / a[i][i];</p><p>  y[i]=(1-w)*x[i]+w*z[i];</p><p><b>  }</b></p><p><b>  e = 0.0;</b></p><p>  for (j = 0; j <

52、n; j++)</p><p><b>  {</b></p><p>  e+= get(x[j],y[j]);</p><p><b>  }</b></p><p>  for (j=0; j < n; j++)</p><p><b>  {</

53、b></p><p>  x[j] = y[j];</p><p><b>  }</b></p><p><b>  k=k+1;</b></p><p><b>  if(k>Q)</b></p><p><b>  {<

54、/b></p><p><b>  r=0;</b></p><p>  printf("算法經過最大迭代次數(shù)還沒有收斂,r=%d\n",r);</p><p><b>  exit(0);</b></p><p><b>  }</b></p&

55、gt;<p>  for(j=0;j<n;j++)</p><p><b>  {</b></p><p>  if(y[j]>100000000||y[j]<-100000000)</p><p><b>  {</b></p><p><b>  r=-

56、2;</b></p><p>  printf("\n迭代方法不收斂,r=%d\n",r);</p><p><b>  exit(0);</b></p><p><b>  }</b></p><p><b>  }</b></p>

57、<p><b>  }</b></p><p>  while (e > v);</p><p>  printf("方程組的解為:\n");</p><p>  for (i = 0; i < n; i++)</p><p>  printf ("x%d=%f\t

58、\n", i+1, x[i]);</p><p>  printf ("\n");</p><p>  printf("\n迭代次數(shù)k=%d\n",k);</p><p><b>  return 0;</b></p><p><b>  }</b>

59、</p><p><b>  float</b></p><p>  get (float x, float y)</p><p><b>  {</b></p><p><b>  float t;</b></p><p>  t=fabs(x-y);

60、</p><p><b>  return t;</b></p><p><b>  }</b></p><p>  2. 高斯LU分解法(源代碼):</p><p>  #include <stdio.h></p><p>  #include <stdl

61、ib.h></p><p>  #define N 100</p><p>  float getmx(float a[N][N], float x[N], int i, int n)</p><p><b>  {</b></p><p>  float mx = 0;</p><p>&

62、lt;b>  int r;</b></p><p>  for(r=i+1; r<n; r++)</p><p><b>  {</b></p><p>  mx += a[i][r] * x[r];</p><p><b>  }</b></p><p&

63、gt;  return mx;</p><p><b>  }</b></p><p>  float getx(float a[N][N], float b[N], float x[N], int i, int n)</p><p><b>  {</b></p><p>  float resu

64、lt;</p><p>  if(i==n-1) </p><p>  result = (float)(b[i]/a[n-1][n-1]);</p><p><b>  else </b></p><p>  result = (float)((b[i]-getmx(a,x,i,n))/a[i][i]);</p&

65、gt;<p>  return result;</p><p><b>  }</b></p><p>  void main()</p><p>  { float l[N][N]={0}; </p><p>  float u[N][N]={0}; </p><p>  floa

66、t y[N]={0};</p><p>  float x[N]={0}; </p><p>  float a[N][N]; </p><p>  float b[N]; </p><p>  float sum=0;</p><p>  int i,j,k;</p><p><b>

67、;  int n;</b></p><p>  int flag=1;</p><p>  printf("請輸入系數(shù)矩陣的大?。?quot;);</p><p>  scanf("%d", &n);</p><p>  printf("請輸入系數(shù)矩陣值:\n");<

68、/p><p>  for(i=0; i<n; i++)</p><p><b>  {</b></p><p>  for(j=0; j<n; j++)</p><p><b>  {</b></p><p>  printf("a[%d][%d]: &qu

69、ot;, i, j);</p><p>  scanf("%f", &a[i][j]);</p><p><b>  }</b></p><p><b>  }</b></p><p>  printf("請輸入右端項數(shù)組:\n");</p>

70、;<p>  for(i=0; i<n; i++)</p><p><b>  {</b></p><p>  printf("b[%d]: ", i);</p><p>  scanf("%f", &b[i]);</p><p><b>  

71、}</b></p><p>  for(i=0; i<n; i++)</p><p><b>  {</b></p><p>  for(j=0; j<n; j++)</p><p><b>  {</b></p><p>  if(i==j) l[i

72、][j] = 1;</p><p><b>  }</b></p><p><b>  }</b></p><p>  for(i=0; i<n; i++)</p><p><b>  {</b></p><p>  u[0][i] = (floa

73、t)(a[0][i]/l[0][0]);</p><p><b>  }</b></p><p>  for(i=0; i<n-1; i++)</p><p><b>  {</b></p><p>  for(j=i+1; j<n; j++)</p><p>&

74、lt;b>  {</b></p><p>  for(k=0,sum=0; k<n; k++)</p><p><b>  {</b></p><p>  if(k != i) sum += l[j][k]*u[k][i];</p><p><b>  }</b></p

75、><p>  l[j][i] = (float)((a[j][i]-sum)/u[i][i]);</p><p><b>  }</b></p><p>  for(j=i+1; j<n; j++)</p><p><b>  {</b></p><p>  for(k=0

76、,sum=0; k<n; k++)</p><p><b>  {</b></p><p>  if(k != i+1) sum += l[i+1][k]*u[k][j];</p><p><b>  }</b></p><p>  u[i+1][j] = (float)((a[i+1][j]

77、-sum));</p><p><b>  }</b></p><p><b>  }</b></p><p><b>  int g, h;</b></p><p><b>  h=0;</b></p><p>  for(g=0

78、;g<n+1;g++)</p><p><b>  {</b></p><p>  if(u[g][g]==0)</p><p><b>  h=h+1;</b></p><p><b>  }</b></p><p><b>  if(

79、h=n)</b></p><p><b>  {</b></p><p>  printf("\nLU decompsition failed.\n\n");</p><p><b>  exit(0);</b></p><p><b>  }</b&

80、gt;</p><p>  /*回代方式計算數(shù)組X*/</p><p>  for(i=n-1; i>=0; i--)</p><p><b>  {</b></p><p>  x[i] = getx(u,y,x,i,n);</p><p><b>  }</b>&l

81、t;/p><p>  printf("\n\n數(shù)組X:\n");</p><p>  for(i=0; i<n; i++)</p><p><b>  {</b></p><p>  printf("x%d = %0.3f\n", i+1,x[i]);</p>&l

溫馨提示

  • 1. 本站所有資源如無特殊說明,都需要本地電腦安裝OFFICE2007和PDF閱讀器。圖紙軟件為CAD,CAXA,PROE,UG,SolidWorks等.壓縮文件請下載最新的WinRAR軟件解壓。
  • 2. 本站的文檔不包含任何第三方提供的附件圖紙等,如果需要附件,請聯(lián)系上傳者。文件的所有權益歸上傳用戶所有。
  • 3. 本站RAR壓縮包中若帶圖紙,網頁內容里面會有圖紙預覽,若沒有圖紙預覽就沒有圖紙。
  • 4. 未經權益所有人同意不得將文件中的內容挪作商業(yè)或盈利用途。
  • 5. 眾賞文庫僅提供信息存儲空間,僅對用戶上傳內容的表現(xiàn)方式做保護處理,對用戶上傳分享的文檔內容本身不做任何修改或編輯,并不能對任何下載內容負責。
  • 6. 下載文件中如有侵權或不適當內容,請與我們聯(lián)系,我們立即糾正。
  • 7. 本站不保證下載資源的準確性、安全性和完整性, 同時也不承擔用戶因使用這些下載資源對自己和他人造成任何形式的傷害或損失。

評論

0/150

提交評論