bisectroot.c
上传用户:bjtelijie
上传日期:2010-01-01
资源大小:87k
文件大小:3k
源码类别:

数学计算

开发平台:

Visual C++

  1. #include <stdio.h>
  2. #include <math.h>
  3. #include <malloc.h>
  4. #include <stdlib.h>
  5. double Func(double);
  6. int BisectRoot(double,double,double,double,double *,int,int *);
  7. void main()
  8. {
  9. int i,n,m;
  10. double a,b,h,eps,*x;
  11. n = 3; /*方程根的个数的预估值*/
  12. x = (double*)calloc(n,sizeof(double)); /*开辟内存空间*/
  13. if(x == NULL)
  14. {
  15. printf("内存分配失败n");
  16. exit(1);
  17. }
  18. a = -3; /*区间起始端点*/
  19. b = 7; /*区间终止端点*/
  20. h = 0.1; /*步长*/
  21. eps = 1.e-8; /*要求达到的精度*/
  22. BisectRoot(a,b,h,eps,x,n,&m); /*调用二分法函数*/
  23. printf("y=sin(x)在范围%2.0f和%2.0f之间的根有%d个根n",a,b,m);
  24. printf("它们分别是:n");
  25. for(i = 0;i<n;i++)
  26. printf("x[%d] = %en",i,x[i]);
  27. free(x); /*释放内存空间*/
  28. }
  29. double Func(double x)
  30. {
  31. return(sin(x));
  32. }
  33. int BisectRoot(a,b,h,eps,x,n,m)
  34. double a; /*实型变量,输入参数,求根区间的起始端点*/
  35. double b; /*实型变量,输入参数,求根区间的终止端点*/
  36. double h; /*利用逐步扫描法确定根位置时的步长*/
  37. double eps; /*实型变量,输入参数,控制精度的参数*/
  38. double *x; /*实型一维数组,输出参数,存放计算得到的数组*/
  39. int n; /*输入参数,区间内方程根的个数的预估值*/
  40. int *m; /*输出参数,实际求得的根的个数*/
  41. {
  42. double z,z0,z1,y,y0,y1;
  43. *m = 0;
  44. z = a;
  45. y = Func(z);
  46. while(1) /*无限循环,直到遇到return或者break语句*/
  47. {/*如果逐步扫描到求根区间的右端点或者得到的根的个数达到预估根的个数*/
  48. if((z>b+h/2)||(*m==n))
  49. return(1);
  50. if(fabs(y)<eps) /*如果当前根z对应的函数f(z)满足精度要求*/
  51. {
  52. *m+=1;
  53. x[*m-1] = z; /*将此时的z值赋值给x数组*/
  54. z+=h/2;
  55. y = Func(z);
  56. continue; /*结束本次循环,即跳过循环体中下面尚未执行
  57.  的语句接着进行下一次是否执行循环的判定*/
  58. }
  59. z1 = z+h; /*逐步扫描中小区间的右端点*/
  60. y1 = Func(z1); /*小区间右端点对应的函数值*/
  61. if(fabs(y1)<eps) /*如果右端点恰好满足根的精度要求*/
  62. {
  63. *m+=1;
  64. x[*m-1] = z1;
  65. z = z1+h/2;
  66. y = Func(z);
  67. continue;
  68. }
  69. if(y*y1>0) /*如果对应根乘积大于零,说明该区间内没有根*/
  70. {
  71. y = y1;
  72. z = z1;
  73. continue;
  74. }
  75. while(1) /*如果本while循环执行,说明逐步扫描小区建z和z1间有根*/
  76. {
  77. if(fabs(z1-z)<eps) /*如果满足精度要求*/
  78. {
  79. *m+=1;
  80. x[*m-1]=(z1+z)/2;
  81. z = z1+h/2;
  82. y = Func(z);
  83. break;
  84. }
  85. z0 = (z1+z)/2; /*二分发求根公式*/
  86. y0 = Func(z0);
  87. if(fabs(y0)<eps)
  88. {
  89. *m = *m+1;
  90. x[*m-1] = z0;
  91. z =z0+h/2;
  92. y = Func(z);
  93. break;
  94. }
  95. if(y*y0<0) /*如果乘积小于零,说明根在z和z0之间*/
  96. {
  97. z1 = z0;
  98. y1 = y0;
  99. }
  100. else /*否则根在z0和z1之间*/
  101. {
  102. z = z0;
  103. y = y0;
  104. }
  105. }
  106. }
  107. }