二分法求解方程例子的说明!!
2008-09-13 10:53
393 查看
正在学习C,今天看的例子是二分法求解方程,源代码如下:
#include <stdio.h>
#include <math.h>
#include <malloc.h>
#include <stdlib.h>
double func(double);
int bisectroot(double,double,double,double,double *,int ,int *);
void main()
{
int i,n,m;
double a,b,h,eps,*x;
n=3;
x=(double *)calloc(n,sizeof(double));
if(x==NULL)
{
printf("memory error!!\n");
exit(1);
}
a=-1;
b=7;
h=0.1;
eps=1.e-8;
bisectroot(a,b,h,eps,x,n,&n);
printf("the root of y=sin(x) between %2.0f and %2.0f have %d\n",a,b,n);
printf("they are:\n");
for(i=0;i<n;i++)
{
printf("x[%d]=%e\n",i,x[i]);
}
free(x);
}
double func(double x)
{
return(sin(x));
}
int bisectroot(double a,double b,double h,double eps,double *x,int n,int *m)
{
double z,z0,z1,y,y0,y1;
*m=0;
z=a;
y=func(z);
while(1)
{
if((z>b+h/2)||(*m==n))
return(1);
if(fabs(y)<eps)
{
*m+=1;
x[*m-1]=z;
z+=h/2;
y=func(z);
continue;
}
z1=z+h;
y1=func(z1);
if(fabs(y1)<eps)
{
*m+=1;
x[*m-1]=z1;
z=z1+h/2;
y=func(z);
continue;
}
if(y*y1>0)
{
y=y1;
z=z1;
continue;
}
while(1)
{
if(fabs(z1-z)<eps)
{
*m+=1;
x[*m-1]=(z1+z)/2;
z=z1+h/2;
y=func(z);
break;
}
z0=(z1+z)/2;
y0=func(z0);
if(fabs(y0)<eps)
{
*m=*m+1;
x[*m-1]=z0;
z=z0+h/2;
y=func(z);
break;
}
if(y*y0<0)
{
z1=z0;
y1=y0;
}
else
{
z=z0;
y=y0;
}
}
}
}
可能有的初学者看了这个程序,会以为这个程序使用于任何函数的求解,其实不是这样.
1.如果有fabs(func(x))和fabs(func(x+h/2))相对于同一点都<eps的情况出现的话,可能会有问题. 2.如果有fabs(func(x))和fabs(func(x+h))在两点的左边且都<eps,再分别加上h/2,形成func(x+h/2)和func(x+h+h/2)在第二个点的两端,就形成了func(x+h/2)*func(x+h+h/2)<0.这样就会重复求得第二个解. 所以,还是要知道要求的函数的大致情况,调整相应的参数.
#include <stdio.h>
#include <math.h>
#include <malloc.h>
#include <stdlib.h>
double func(double);
int bisectroot(double,double,double,double,double *,int ,int *);
void main()
{
int i,n,m;
double a,b,h,eps,*x;
n=3;
x=(double *)calloc(n,sizeof(double));
if(x==NULL)
{
printf("memory error!!\n");
exit(1);
}
a=-1;
b=7;
h=0.1;
eps=1.e-8;
bisectroot(a,b,h,eps,x,n,&n);
printf("the root of y=sin(x) between %2.0f and %2.0f have %d\n",a,b,n);
printf("they are:\n");
for(i=0;i<n;i++)
{
printf("x[%d]=%e\n",i,x[i]);
}
free(x);
}
double func(double x)
{
return(sin(x));
}
int bisectroot(double a,double b,double h,double eps,double *x,int n,int *m)
{
double z,z0,z1,y,y0,y1;
*m=0;
z=a;
y=func(z);
while(1)
{
if((z>b+h/2)||(*m==n))
return(1);
if(fabs(y)<eps)
{
*m+=1;
x[*m-1]=z;
z+=h/2;
y=func(z);
continue;
}
z1=z+h;
y1=func(z1);
if(fabs(y1)<eps)
{
*m+=1;
x[*m-1]=z1;
z=z1+h/2;
y=func(z);
continue;
}
if(y*y1>0)
{
y=y1;
z=z1;
continue;
}
while(1)
{
if(fabs(z1-z)<eps)
{
*m+=1;
x[*m-1]=(z1+z)/2;
z=z1+h/2;
y=func(z);
break;
}
z0=(z1+z)/2;
y0=func(z0);
if(fabs(y0)<eps)
{
*m=*m+1;
x[*m-1]=z0;
z=z0+h/2;
y=func(z);
break;
}
if(y*y0<0)
{
z1=z0;
y1=y0;
}
else
{
z=z0;
y=y0;
}
}
}
}
可能有的初学者看了这个程序,会以为这个程序使用于任何函数的求解,其实不是这样.
1.如果有fabs(func(x))和fabs(func(x+h/2))相对于同一点都<eps的情况出现的话,可能会有问题. 2.如果有fabs(func(x))和fabs(func(x+h))在两点的左边且都<eps,再分别加上h/2,形成func(x+h/2)和func(x+h+h/2)在第二个点的两端,就形成了func(x+h/2)*func(x+h+h/2)<0.这样就会重复求得第二个解. 所以,还是要知道要求的函数的大致情况,调整相应的参数.
相关文章推荐
- TF随笔-14-二分法求解一元方程
- 二分法求解方程根
- 二分法求解方程(有错误,请高手指点)。
- 数值计算方法:二分法求解方程的根(伪代码 python c/c++)
- 使用二分法求解一元N次方程的近似值
- 二分法求解方程的近似解(sicily 1017)
- 二分法求解方程的解
- 关于使用牛顿迭代法和二分法解方程的算法说明
- 金融库QuantLib中最优化求解方程Simplex用法--例子
- 非线性方程的数值解法——二分法求解
- 二分法 简单迭代法 Newton法 弦截法 求解非线性方程的根
- 第三周练习——二分法2 方程求解
- 求解一元多次方程的两种方法:牛顿迭代法和二分法
- C++实现二分法求解方程
- 1、编写程序,分别用二分法和牛顿迭代法求解方程x3 – 3x – 1 = 0在x = 2附近的实根,要求计算精确到小数点后七位数字为止,并将求出的近似结果与理论值2cos20 相比较,二分法的初始迭代
- /* 编程用二分法求解方程x3+4x2-10=0的解。 */
- 关于使用SVD分解方法求解AX=0方程的一点说明
- PCA、SVD、协方差矩阵求解的关系和对比(例子说明)
- 二分法,matlab中利用二分法求解一个多项式方程的近似值。
- 求解方程(二分法)