数值分析课程设计(求解线性方程组)_第1页
已阅读1页,还剩24页未读 继续免费阅读

下载本文档

版权说明:本文档由用户提供并上传,收益归属内容提供方,若内容存在侵权,请进行举报或认领

文档简介

1、<p><b>  数值分析 课程设计</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>  自定义函数 SOR(A, B, w, MAXN, TOL),以实现SOR方法求解线性方程组AX=B,其中</p><p><b>  A——系数矩阵;<

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

4、gt;<p>  -2:SOR方法不收敛;(不收敛的依据为的某个分量值超出区间[-108, 108]。)</p><p>  -1:矩阵有一列全为0;</p><p>  0:算法经过MAXN次迭代还未收敛;</p><p>  k:SOR方法经k次迭代收敛,求得方程组的解向量X记录下来.</p><p>  自定义函数Dire

5、ct(A, B),以实现高斯LU分解的方法求解线性方程组AX=B,其中</p><p><b>  A——系数矩阵;</b></p><p><b>  B——常数列向量;</b></p><p><b>  返回值有两种可能:</b></p><p>  “LU decomp

6、sition failed.”:分解过程中U的对角线元素至少一个为0;</p><p><b>  X:分解过程中</b></p><p>  分别使用步骤1中定义的函数SOR(A, B, w, MAXN, TOL)和步骤2中定义的函数Direct(A, B)进行测试,记录返回值及X值(算法收敛或有效的情形, 保留4位小数):</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>  将松弛系数w的不同德值代入计算出X的值。</p><p>  利用迭代使|X(k+1)-X(k)|<e(e为精度,e=10^(-9))。</p><p>  (2)由于矩阵出现了,一列为0,所以不能使用迭

13、代,在程序中会出现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)将方程组变为:</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>  将松弛系数w的不同德值代入计算出X的值。</p><p>  利用迭代使|X(k+1)-X(k)|<e(e为精度,e=10^(-9))。</p><p>  (5) 将方程组变为:</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>  将松弛系数w的不同德值代入计算出X的值。</p><p>  利用迭代使|X(k+1)-X(k)|<e(e为精度,e=10^

22、(-9))。</p><p>  (6)将方程组变为:</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>  将松弛系数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>  (1)W=1的时候:</p><p>  (2)W=1.05的时候:</p><p>  (3)W=1.1的时候:</p><p>  (4

25、)W=1.2的时候</p><p>  (5)W=1.3的时候</p><p>  (6)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>  (2)W=1.1的时候</p><p>  (3)W=1.3的时候</p><p>  (4)W=1.8的时候</p><p><b>  测试5</

27、b></p><p><b>  (1)W=1的时候</b></p><p>  (2)W=1.3的时候</p><p>  (3)W=1.6的时候</p><p>  (4)W=1.9的时候</p><p><b>  测试6</b></p><p

28、><b>  (1)W=1的时候</b></p><p>  (2)W=1.3的时候</p><p>  (3)W=1.6的时候</p><p>  (4)W=1.8的时候</p><p>  (5)W=1.9的时候</p><p>  (二) 高斯LU方法</p><p

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

30、L矩阵为:</b></p><p><b>  U矩阵:</b></p><p>  由于U矩阵对角线出现了0,所以出现了“LU decompsition failed.”</p><p><b>  (3)L矩阵:</b></p><p><b>  U矩阵:</b&g

31、t;</p><p>  由于U矩阵对角线出现了0,所以出现了“LU decompsition failed.”</p><p>  (4)将矩阵A变为:</p><p>  而在将矩阵变为L与U时候出现了异常,在L与U矩阵中有异常值,且在U矩阵对角线上出现了0值,所以出现了“LU decompsition failed.”</p><p>

32、<b>  (5)L矩阵:</b></p><p>  利用LD=B,可得到D为:</p><p><b>  U矩阵为:</b></p><p>  利用UX=D,求出X的值出来。 </p><p><b>  (6)L矩阵:</b></p><p>

33、;<b>  (见下一页)</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法设计分析和测试以后,我发现两种方法各有优劣。高斯LU法得出的结果精度比较高,但是却不适用于所有的方程,使用范围相对较窄;而在使用SOR法时,虽然精度会稍微差一点,但是通过调整松弛度w,却可能适用于更多的方程,适用范围相对较宽。</p><p> 

36、 在做课程设计这个过程中,我发现自己还有很多很多知识没有学好,在参考别人的例子的时候好像很简单,但自己一上机操作写程序的时候就出现问题。调试的时候系统总有报错,还有很多警告,每修改一个变量,往往都要调试很久,有时候仅仅只是少了一个大括号,却地花上近半个小时才能找到问题的瓶颈所在。此外,通过本次的课程设计,我重温了许多C语言的知识。同时也发现了自己对C语言的掌握程度有所下降。其实,为了更加方便我今后的学习,我还是有必要对MATLAB进行学

37、习,不断扩充我的知识。</p><p>  最后,虽然得到的结果未如理想,但我会继续努力!谢谢李老师的指导!</p><p><b>  五、参考文献</b></p><p>  【1】 《标准C语言基础教程》(第四版) 电子工业出版社</p><p>  【2】 《数值分析》(第三版)

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("请输入系数矩阵的阶数:"

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("请输入最大迭代次数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("请输入松弛系数w=");</p><p>  scanf("%f",&w);</p><p>  printf ("请输入系数矩阵:\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 ("请输入右端项数组: \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("算法经过最大迭代次数还没有收敛,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迭代次数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("请输入系数矩阵的大小:");</p><p>  scanf("%d", &n);</p><p>  printf("请输入系数矩阵值:\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("请输入右端项数组:\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>  /*回代方式计算数组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数组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. 本站的文档不包含任何第三方提供的附件图纸等,如果需要附件,请联系上传者。文件的所有权益归上传用户所有。
  • 3. 本站RAR压缩包中若带图纸,网页内容里面会有图纸预览,若没有图纸预览就没有图纸。
  • 4. 未经权益所有人同意不得将文件中的内容挪作商业或盈利用途。
  • 5. 众赏文库仅提供信息存储空间,仅对用户上传内容的表现方式做保护处理,对用户上传分享的文档内容本身不做任何修改或编辑,并不能对任何下载内容负责。
  • 6. 下载文件中如有侵权或不适当内容,请与我们联系,我们立即纠正。
  • 7. 本站不保证下载资源的准确性、安全性和完整性, 同时也不承担用户因使用这些下载资源对自己和他人造成任何形式的伤害或损失。

评论

0/150

提交评论