2016-09-12 2 views
1

Я написал короткую программу на C, чтобы выполнить линейную интерполяцию, которая итерации дает корень функции заданному числу десятичных точек. До сих пор для функции f (x):Почему C дает разные значения для Python для этого алгоритма?

long double f(long double x) { 
     return (pow(x, 3) + 2 * x - 8); 
    } 

Программа не сходится к значению 1dp. Программа обновляет переменные a и b, между которыми лежит корень f (x), до тех пор, пока a и b не будут округлены до одного и того же числа с заданной точностью. Использование длинных дублей и выше функции, отладчик показывает в течение первых 2-х итераций:

a = 1.5555555555555556 
    a = 1.6444444444444444 

, хотя это должно было быть:

a = 1.5555555555555556 
    a = 1.653104925053533 

Программа не может изменять значения после этого. Уравнение для линейной интерполяции, которое я использую, представляет собой перестроенную версию математической заданной here, а код, который я использую, - это C-версия программы-питона, которую я написал. Почему реализация C получает разные значения, несмотря на тот же алгоритм, и как я могу это исправить?

ОК я все еще получаю повесить, но, надеюсь, ниже меня Minimal, полный и проверяемый пример:

#include <stdlib.h> 
#include <stdio.h> 
#include <math.h> 

long double a; long double b; long double c; // The values for interpolation 
long double fofa; long double fofb; long double fofc; // The values f(a), f(b) and f(c) 
const int dp = 1; // The number of decimal places to be accurate to 

long double f(long double x) { 
    return (pow(x, 3) + 2 * x - 8); 
} 

int main(void) { 
    a = 1; b = 2; 
    while(roundf(a * pow(10, dp))/pow(10, dp) != roundf(b * pow(10, dp))/pow(10, dp)) { // While a and b don't round to the same number... 

     fofa = f(a); fofb = f(b); // Resolve the functions 
     printf("So f(a) = %g, f(b) = %g\n", (double)fofa, (double)fofb); // Print the results 

     c = (b * abs(fofa) + a * abs(fofb))/(abs(fofb) + abs(fofa)); // Linear Interpolation 

     fofc = f(c); 

     if(fofc < 0) { 
      a = c; 
     } 
     else if(fofc == 0) { 
      a = c; 
      break; 
     } 
     else { 
      b = c; 
     } 

    } 

    printf("The root is %g, to %d decimal places, after %d iterations.\n", (double)a, dp, i); 
} 
+0

не могли бы вы использовать некоторые отладки, чтобы увидеть, какие значения вы получаете для 'fofa' (который, как я полагаю,' f (a) '),' fofb' и все остальные промежуточные значения? – Teepeemm

+0

Что вы используете как значение для 'dp' (1, 10, 15, что-то еще)? Почему вы используете 'float'-precision с' roundf() ', когда у вас есть' long double' для 'a' и' fofa' и т. Д.? (Вы должны использовать 'roundl()' и 'powl()' (а не 'pow()'), если вы не используете '', но вы должны упомянуть, что если вы и не хотите используйте 'roundf()', но просто 'round()' ... –

+1

Пожалуйста, предоставьте [mcve], а не игру-головоломку. Также см. [ask] – Olaf

ответ

7

Функция abs() (от <stdlib.h>) имеет подпись int abs(int); - вы получая целые числа из ваших вычислений.

Вы должны использовать long double fabsl(long double); от <math.h>.

Вы также должны использовать powl() вместо pow() (long double против double) и roundl() вместо roundf() (long double против float).

Убедитесь, что вы используете правильные типы, другими словами.


Когда вы исправили проблемы типа, у вас все еще есть проблема с конвергенцией. Это помогло бы, если бы вы создали MCVE (Minimal, Complete, Verifiable Example). Тем не менее, это MCVE, что я могу вывести из вашего вопроса:

#include <math.h> 
#include <stdio.h> 

static inline long double f(long double x) 
{ 
    return(powl(x, 3) + 2 * x - 8); 
} 

int main(void) 
{ 
    long double a = 1.0L; 
    long double b = 2.0L; 
    int dp = 6; 

    while (roundl(a * powl(10, dp))/powl(10, dp) != roundl(b * powl(10, dp))/powl(10, dp)) 
    { 
     long double fofa = f(a); 
     long double fofb = f(b); 
     long double c = (b * fabsl(fofa) + a * fabsl(fofb))/(fabsl(fofb) + fabsl(fofa)); 
     long double fofc = f(c); 
     printf("a = %+.10Lf, f(a) = %+.10Lf\n", a, fofa); 
     printf("b = %+.10Lf, f(b) = %+.10Lf\n", b, fofb); 
     printf("c = %+.10Lf, f(c) = %+.10Lf\n", c, fofc); 
     putchar('\n'); 

     if (fofc < 0.0L) 
     { 
      a = c; 
     } 
     else if (fofc == 0.0L) 
     { 
      a = c; 
      break; 
     } 
     else 
     { 
      b = c; 
     } 
    } 
    printf("Result: a = %Lg\n", a); 
    return 0; 
} 

Выход я получаю от него:

a = +1.0000000000, f(a) = -5.0000000000 
b = +2.0000000000, f(b) = +4.0000000000 
c = +1.5555555556, f(c) = -1.1248285322 

a = +1.5555555556, f(a) = -1.1248285322 
b = +2.0000000000, f(b) = +4.0000000000 
c = +1.6531049251, f(c) = -0.1762579238 

a = +1.6531049251, f(a) = -0.1762579238 
b = +2.0000000000, f(b) = +4.0000000000 
c = +1.6677455452, f(c) = -0.0258828049 

a = +1.6677455452, f(a) = -0.0258828049 
b = +2.0000000000, f(b) = +4.0000000000 
c = +1.6698816424, f(c) = -0.0037639074 

a = +1.6698816424, f(a) = -0.0037639074 
b = +2.0000000000, f(b) = +4.0000000000 
c = +1.6701919841, f(c) = -0.0005465735 

a = +1.6701919841, f(a) = -0.0005465735 
b = +2.0000000000, f(b) = +4.0000000000 
c = +1.6702370440, f(c) = -0.0000793539 

a = +1.6702370440, f(a) = -0.0000793539 
b = +2.0000000000, f(b) = +4.0000000000 
c = +1.6702435859, f(c) = -0.0000115206 

a = +1.6702435859, f(a) = -0.0000115206 
b = +2.0000000000, f(b) = +4.0000000000 
c = +1.6702445357, f(c) = -0.0000016726 

a = +1.6702445357, f(a) = -0.0000016726 
b = +2.0000000000, f(b) = +4.0000000000 
c = +1.6702446735, f(c) = -0.0000002428 

a = +1.6702446735, f(a) = -0.0000002428 
b = +2.0000000000, f(b) = +4.0000000000 
c = +1.6702446936, f(c) = -0.0000000353 

a = +1.6702446936, f(a) = -0.0000000353 
b = +2.0000000000, f(b) = +4.0000000000 
c = +1.6702446965, f(c) = -0.0000000051 

a = +1.6702446965, f(a) = -0.0000000051 
b = +2.0000000000, f(b) = +4.0000000000 
c = +1.6702446969, f(c) = -0.0000000007 

a = +1.6702446969, f(a) = -0.0000000007 
b = +2.0000000000, f(b) = +4.0000000000 
c = +1.6702446970, f(c) = -0.0000000001 

a = +1.6702446970, f(a) = -0.0000000001 
b = +2.0000000000, f(b) = +4.0000000000 
c = +1.6702446970, f(c) = -0.0000000000 

a = +1.6702446970, f(a) = -0.0000000000 
b = +2.0000000000, f(b) = +4.0000000000 
c = +1.6702446970, f(c) = -0.0000000000 

a = +1.6702446970, f(a) = -0.0000000000 
b = +2.0000000000, f(b) = +4.0000000000 
c = +1.6702446970, f(c) = -0.0000000000 

a = +1.6702446970, f(a) = -0.0000000000 
b = +2.0000000000, f(b) = +4.0000000000 
c = +1.6702446970, f(c) = -0.0000000000 

Причина бесконечного цикла ясно; разница между a и b не является небольшой долей. Вам нужно пересмотреть свой критерий конвергенции. Вероятно, это должно быть сравнение fofc с 0.0 с точностью до заданного количества знаков после запятой - или что-то вдоль этих линий.

+0

Спасибо, очень удивительно, что сейчас он сходится (слишком много dp, но я это увижу!) – AWSM

+0

Следуя подсказкам из [Lutzl] (http : //stackoverflow.com/users/3088138/lutzl) 's [answer] (http://stackoverflow.com/a/39459981), вы можете найти статью Wikipedia articl e по методу [Ложная позиция] (https://en.wikipedia.org/wiki/False_position_method). Я еще не попытался исправить математику. –

3

То, что вы используете, называется regula falsi или метод ложного расположения.

Фактически нет необходимости использовать абсолютные значения, если поддерживается противоположное условие знака.

Существует довольно известная проблема с обычным ванильным регуляционным фальситом, заключающаяся в том, что в тот момент, когда функция является выпуклой над оставшимся интервалом, одна конечная точка больше не будет перемещаться к корню. Существуют легкие модификации, чтобы избежать этого, например, вставляя шаги деления пополам. Еще более простой в применении, но сложнее понять, это модификация Иллинойса . См. Статью в Wikipedia для regula falsi.

Или этот вопрос и ответ: Regula-Falsi Algorithm?


адаптировано из ответа в ссылке:

#include<math.h> 
#include<stdio.h> 

long double f(long double x) { 
    return powl(x, 3) + 2 * x - 8; 
} 

int main(void) { 
    const int dp = 18; 
    long double eps=0.5*powl(10,-dp); 
    int i=0; 
    long double a=1, fofa = f(a); 
    long double b=2, fofb = f(b); 

    printf("\na=%.21Lf b=%.21Lf fofa=%.21Lf fofb=%.21Lf\n------\n",a,b, fofa,fofb); 

    if(signbit(fofb)==signbit(fofa)) { 
     printf("Warning, initial values do not have opposite sign!\n"); 
    } 
    do { 
     long double c=(a*fofb-b*fofa)/(fofb-fofa), fofc = f(c); 

     if(signbit(fofc)!=signbit(fofa)) { 
      b=a; fofb=fofa; 
      a=c; fofa=fofc; 
     } else { 
      a=c; fofa=fofc; 
      fofb *= 0.5; 
     } 
     i++; 
     printf("\na=c=%.21Lf b=%.21Lf fofa=fofc=%.21Lf fofb=%.21Lf",c,b, fofc,fofb); 

    } while(fabsl(b-a)>eps); 

    printf("\ngoal reached after %d iterations\n",i); 

    return 0; 
} 

с результатом

a=1.000000000000000000000 b=2.000000000000000000000 fofa=-5.000000000000000000000 fofb=4.000000000000000000000 
------ 

a=c=1.555555555555555555507 b=2.000000000000000000000 fofa=fofc=-1.124828532235939643549 fofb=2.000000000000000000000 
a=c=1.715539947322212467064 b=1.555555555555555555507 fofa=fofc=0.480046589479470829469 fofb=-1.124828532235939643549 
a=c=1.667685780603345490963 b=1.715539947322212467064 fofa=fofc=-0.026500999000164700194 fofb=0.480046589479470829469 
a=c=1.670189362207942139265 b=1.715539947322212467064 fofa=fofc=-0.000573759143326624515 fofb=0.240023294739735414734 
a=c=1.670297511133477909220 b=1.670189362207942139265 fofa=fofc=0.000547652143260468627 fofb=-0.000573759143326624515 
a=c=1.670244695550498326532 b=1.670297511133477909220 fofa=fofc=-0.000000014643676336194 fofb=0.000547652143260468627 
a=c=1.670244696962696986627 b=1.670297511133477909220 fofa=fofc=-0.000000000000373724489 fofb=0.000273826071630234313 
a=c=1.670244696962769068724 b=1.670244696962696986627 fofa=fofc=0.000000000000373706274 fofb=-0.000000000000373724489 
a=c=1.670244696962733028543 b=1.670244696962769068724 fofa=fofc=-0.000000000000000000434 fofb=0.000000000000373706274 
a=c=1.670244696962733028651 b=1.670244696962733028543 fofa=fofc=0.000000000000000000867 fofb=-0.000000000000000000434 
goal reached after 10 iterations 
Смежные вопросы