i am trying to find such p,that for given function f(p),we have equality
p=f(p); here is code
#include <iostream>
#include<math.h>
using namespace std;
float fixed(float x){
return (float)(pow(x,3)-4*pow(x,2)-10);
}
int main(){
float p=0.0;
float p0=1.开发者_运维问答5;
float tol=(float).001;
int N=25;
int i=1;
while(i<N){
p=(float)fixed(p0);
if((p-p0)<tol){
cout<<p<<endl;
break;
}
i=i+1;
p0=p;
if(i>N){
cout<<"solution not found ";
break;
}
}
return 0;
}
i have tried different initial point,also different tolerance,but result is very nonsense -16 or -15.something,so what is wrong?is code correct?please help
I think you simply don't have a situation in which the iterative algorithm is applicable. See here for some conditions. Your function doesn't have an attractive fixed point near 1.5, and the algorithm diverges.
But why the numerics: Your function is f(x) = x^3 - 4x - 10
, so solving f(x) = x
amounts to finding the zeros of f(x) - x
, and there is only one real zero near 5.35. Howevever, f'(x)
at that point is very large, so even there the iterative algorithm isn't usable.
A numeric root-finding algorithm may be a more appropriate approach.
I'm not sure what you're trying to achieve, but probably using :
if(fabs(p-p0)<tol){
instead of :
if((p-p0)<tol){
will get you closer to where you want to go.
The reason is that all negative values are smaller than .001
, so if the result of the function is a negative value (like in this case), you'll stop immediately.
Btw : this check :
if(i>N){
will never be true. You probably meant to use ==
or >=
instead of >
.
fabs(p - p0) <= tol * fabs(p);
Is the correct formula for equal values. In case of special ranges you also have to care for NaN and Inf.
#include <limits>
#include <cfloat>
inline bool Equal(const double& x, const double& y)
{
const int classX (_fpclass(x));
if(classX != _fpclass(y)) {
return false;
}
switch(classX) {
case _FPCLASS_SNAN:
case _FPCLASS_QNAN:
case _FPCLASS_NINF:
case _FPCLASS_NZ :
case _FPCLASS_PZ :
case _FPCLASS_PINF:
return true;
break;
case _FPCLASS_NN :
case _FPCLASS_ND :
case _FPCLASS_PD :
case _FPCLASS_PN :
default:
break;
}
return fabs(x - y) <= std::numeric_limits<double>::epsilon() * fabs(x);
}
精彩评论