Краткое объяснение проблемы: я использую алгоритм Ньютона-Рафсона для поиска корней в полиномах и в некоторых случаях не работает. Зачем?
Я взял из «числовых рецептов в c ++» гибридный алгоритм Ньютона-Рафсона, который делит пополам, если New-Raph не сходится должным образом (с низким значением производной или если скорость сходимости не быстрая).
Я проверил алгоритм с несколькими полиномами, и он работал. Сейчас я тестирую в программном обеспечении, которое у меня есть, и всегда получаю ошибку с определенным полиномом. Моя проблема в том, что я не знаю, почему этот многочлен просто не достигает результата, когда так много других. Поскольку я хочу улучшить алгоритм для любого полинома, нужно знать, какой из них является причиной отсутствия сходимости, чтобы я мог правильно его обработать.
Далее я опубликую всю информацию, которую я могу предоставить об алгоритме и полиноме, в котором у меня есть ошибка.
Полином:
f(t)= t^4 + 0,557257315256597*t^3 - 3,68254086033178*t^2 +
+ 0,139389107255627*t + 1,75823776590795
Это первая производная:
f'(t)= 4*t^3 + 1.671771945769790*t^2 - 7.365081720663563*t + 0.139389107255627
Участок:
Корни (автор Matlab):
-2.133112008595826 1.371976341295347 0.883715461977390
-0.679837109933505
Алгоритм:
double rtsafe(double* coeffs, int degree, double x1, double x2,double xacc,double xacc2)
{
int j;
double df,dx,dxold,f,fh,fl;
double temp,xh,xl,rts;
double* dcoeffs=dvector(0,degree);
for(int i=0;i<=degree;i++)
dcoeffs[i]=0.0;
PolyDeriv(coeffs,dcoeffs,degree);
evalPoly(x1,coeffs,degree,&fl);
evalPoly(x2,coeffs,degree,&fh);
evalPoly(x2,dcoeffs,degree-1,&df);
if ((fl > 0.0 && fh > 0.0) || (fl < 0.0 && fh < 0.0))
nrerror("Root must be bracketed in rtsafe");
if (fl == 0.0) return x1;
if (fh == 0.0) return x2;
if (fl < 0.0) { // Orient the search so that f(xl) < 0.
xl=x1;
xh=x2;
} else {
xh=x1;
xl=x2;
}
rts=0.5*(x1+x2); //Initialize the guess for root,
dxold=fabs(x2-x1); //the "stepsize before last,"dx=dxold; //and the last step
evalPoly(rts,coeffs,degree,&f);
evalPoly(rts,dcoeffs,degree-1,&dx);
for (j=1;j<=MAXIT;j++) { //Loop over allowed iterations
if ((((rts-xh)*df-f)*((rts-xl)*df-f) > 0.0) //Bisect if Newton out of range,
|| (fabs(2.0*f) > fabs(dxold*df))) { //or not decreasing fast enough.
dxold=dx;
dx=0.5*(xh-xl);
rts=xl+dx;
if (xl == rts)
return rts; //Change in root is negligible.
} else {// Newton step acceptable. Take it.
dxold=dx;
dx=f/df;
temp=rts;
rts -= dx;
if (temp == rts)
return rts;
}
if (fabs(dx) < xacc)
return rts;// Convergence criterion
evalPoly(rts,coeffs,degree,&f);
evalPoly(rts,dcoeffs,degree-1,&dx);
//The one new function evaluation per iteration.
if (f < 0.0) //Maintain the bracket on the root.
xl=rts;
else
xh=rts;
}
//As the Accuracy asked to the algorithm is really high (but usually easily reached)
//the results precission is checked again, but with a less exigent result
dx=f/df;
if(fabs(dx)<xacc2)
return rts;
nrerror("Maximum number of iterations exceeded in rtsafe");
return 0.0;// Never get here.
}
Алгоритм вызывается со следующими переменными:
x1=0.019
x2=1.05
xacc=1e-10
xacc2=0.1
degree=4
MAXIT=1000
coeffs[0]=1.75823776590795;
coeffs[1]=0.139389107255627;
coeffs[2]=-3.68254086033178;
coeffs[3]=0.557257315256597;
coeffs[4]=1.0;
проблема в том, что алгоритм превышает максимальные итерации и существует корень ошибки приблизительно 0.15
,
Поэтому мой прямой и краткий вопрос: почему этот многочлен не достигает точной ошибки, когда многие (не менее 1000) других очень похожих многочленов делают (с точностью 1е-10 точности и нескольких итераций!)
Я знаю, что вопрос сложный, и, возможно, у него нет действительно прямого ответа, но я застрял с этим на несколько дней и не знаю, как его решить. Большое спасибо за то, что нашли время прочитать мой вопрос.
Я не совсем уверен, почему, но проверка того, достаточно ли быстро уменьшается функция, в данном случае не работает.
Это работает, если я делаю это так:
double old_f = f;
.
.
.
if ((((rts-xh)*df-f)*((rts-xl)*df-f) > 0.0) //Bisect if Newton out of range,
|| (fabs(2.0*f) > old_f)) { //or not decreasing fast enough.
.
.
.
if (fabs(dx) < xacc)
return rts;// Convergence criterion
old_f = f;
ОБНОВИТЬ
Похоже, в вашем коде есть проблема:
evalPoly(rts,dcoeffs,degree-1,&dx);
должно быть
evalPoly(rts,dcoeffs,degree-1,&df);
Не запуская ваш код, я предполагаю, что вы сравниваете значения с плавающей запятой на равенство, чтобы определить, сходится ли ваше решение.
if (xl == rts)
return rts; //Change in root is negligible.
Может быть, вы должны рассчитать это как соотношение:
diff = fabs(xl - rts);
if (diff/xl <= 1.0e-8) // pick your own accuracy value here
return rts;