Табулирование функции на отрезке с дробным шагом (накапливается погрешность) , идеально чистый СИ, visual studio 2010

Всем хай! Без раскачки сразу переходим к делу!

Есть функция y = 2x^2 + 1, на отрезке [2; 3] надо протабулировать. Кол-во разбиний отрезка (n - натуральное) задается с клавиатуры и может быть от 1 до infinity.

Вот код всей прожки:

#include <stdio.h>  
#include <conio.h>      
#include <locale.h>    
#include <math.h>      
 
#define X_LEFT_BORDER 2.0
#define X_RIGHT_BORDER 3.0
#define EPS 0.0001
 
// заданная функция y = 2x^2 + 1
double f(const double px)
{
    return (2 * pow(px, 2) + 1);
}
// главная функция программы
int main(void)
{
    int n;
    float dx;
    float x;
    float y;
    int number;
 
    setlocale(LC_ALL, "rus");
    
    do
    {
        printf("Введите количество разбиений n(натуральное число) отрезка [%0.0f .. %0.0f]: ", X_LEFT_BORDER, X_RIGHT_BORDER);
        scanf("%d", &n);
    }
    while(n < 1);
 
    dx = (X_RIGHT_BORDER - X_LEFT_BORDER) / (float)n;
 
    printf("\nТабулирование функции y = 2x^2 + 1 на отрезке [%0.0f; %0.0f] с шагом dx = %0.4f:\n", X_LEFT_BORDER, X_RIGHT_BORDER, dx);
    printf("\t------------------------------\n");
    printf("\t| № | Аргумент х | Функция y |\n");
    printf("\t------------------------------\n");
    number = 0;
    for(x = X_LEFT_BORDER; x <= (X_RIGHT_BORDER + EPS); x += dx)   // <------------ ПРОБЛЕМА ЗДЕСЬ!
    {
        number++;
        y = f(x);
        printf("\t|%3d| %10.4f | %9.4f |\n", number, x, y);
    }
    printf("\t------------------------------");
 
    getch();
    return 0;
}

схватил лютую проблему в этой строке:

for(x = X_LEFT_BORDER; x <= (X_RIGHT_BORDER + EPS); x += dx)

вообще, по классике должно быть так:

for(x = X_LEFT_BORDER; x <= X_RIGHT_BORDER; x += dx)

результаты (см.рис.).

Прикреплённая картинка

При n = 5 последняя итерация не учитывается. Отладчик показал, что x накапливается с погрешностью какой-то (ПРИЧЕМ В БОЛЬШУЮ СТОРОНУ, а когда n = 10, ТО В МЕНЬШУЮ, я этого вообще не понимаю):
x = 2.000000, x + dx = 2.200000, x + dx = 2.4000001, x + dx = 2.6000001, x + dx = 2.8000002
я кое-что знаю про машинное епсилон. Но мне от этого НЕ легче!

Пробовал победить проблему ну, наверное, 5-ю различными способами:
1. убрал константы в переменные float/double - НЕ ПОМОГЛО
2. добавил "хитрое" условие окончания цикла:

x <= X_RIGHT_BORDER + dx/2.0

- ПОМОГАЕТ на некоторых n, но при БОЛЬШИХ n, где-то 10000 - та же проблема
3. в условии добавлял типа fabs(X_RIGHT_BORDER - x) >= EPS - очевидно, что работать НЕ будет для любых n, т к при некоторых n длина отрезка разбиения может стать вообще меньше EPS и в некоторых случаях приводит к бесконечному циклу
4. приводил явно все к double/float - разумеется, НЕ ПОМОГЛО)
5. что-то еще

Посмотрел несколько прожек в инете, там все легко и просто, пишут аля

for(x = X_LEFT_BORDER; x <= X_RIGHT_BORDER; x += dx)

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

А как нужно сделать условие выхода из цикла, чтобы прожка РАБОТАЛА НА ЛЮБЫХ РАЗУМНЫХ n правильно?? Подскажите как быть-то?

спс.за внимание

Добавлено 05.12.19, 19:56
P.S. забыл дописать, что нужно именно через for)

Сообщение отредактировано: FasterHarder - 05.12.19, 19:56

FasterHarder, может имеет смысл в цикле задавать просто перебор шагов?
Ну а значение x вычислять как dx*номер_шага. Таким образом у тебя на любом
из шагов будет одинаковая погрешность, равная погрешности операции умножения.

Не ты первый, не ты последний.
Правильнее – заранее рассчитать количество итераций как целое число, использовать его для организации цикла, а очередную итерацию рассчитывать линейной интерполяцией между минимумом и максимумом по текущему номеру точки из общего их количества.

Добавлено 05.12.19, 20:09

Цитата JoeUser @ 05.12.19, 20:02

FasterHarder, может имеет смысл в цикле задавать просто перебор шагов?

Вот типа этого, да.

Цитата JoeUser @ 05.12.19, 20:02

FasterHarder, может имеет смысл в цикле задавать просто перебор шагов?
Ну а значение x вычислять как dx*номер_шага. Таким образом у тебя на любом
из шагов будет одинаковая погрешность, равная погрешности операции умножения.

я понял вроде, попробовал, вот так, да?

    //for(x = X_LEFT_BORDER; x <= X_RIGHT_BORDER; x += dx)
    for(i = 0; i <= n; i++)
    {
        //number++;
        //y = f(x);
        x = X_LEFT_BORDER + dx * i;
        y = f(x);
        printf("\t|%3d| %10.4f | %9.4f |\n", (i + 1), x, y);
    }
    printf("\t------------------------------");

сейчас, да, работает чОтко! спс.

но, такой код немного не соот-ет духу табулирования функции) не-не, я оставлю такой вариант, если нельзя сделать так:

for(x = X_LEFT_BORDER; x <= X_RIGHT_BORDER; x += dx)

просто на форумах в сети все показывают именно вариант через for, в котором "скачет" х и все ликуют, мол, ох как хорошо, что все работает исправно. У меня вот не получилось...

а в твоем варианте, как я понимаю, происходит "обнуление" погрешности на каждой итерации, поэтому она НЕ накапливается к концу вычислений. ну, ок)

Добавлено 05.12.19, 20:28
JoeUser, Qraizer

этот вопрос чисто из любопытства (ну, хочется понять эту тонкость)

При n = 5 последняя итерация не учитывается. Отладчик показал, что x накапливается с погрешностью какой-то (ПРИЧЕМ В БОЛЬШУЮ СТОРОНУ, а когда n = 10, ТО В МЕНЬШУЮ, я этого вообще не понимаю)

т е, когда dx = 0.1 (n = 10), x = 2,0000, x + dx = 2.0999999, т е В МЕНЬШУЮ СТОРОНУ (не добирает)
когда dx = 0.2 (n = 5), x = 2.00000. x + dx = 2.200000001, т е в БОЛЬШУЮ СТОРОНУ (перебирает)

почему так происходит??

Сообщение отредактировано: FasterHarder - 05.12.19, 20:28

Цитата FasterHarder @ 05.12.19, 20:21

почему так происходит??

Согласно информатике:

Цитата

Сложение чисел с плавающей точкой

Необходимым условием выполнения этой операции является соответствие разрядов мантисс складываемых чисел, что достигается выравниванием их порядков путем сдвига мантиссы меньшего числа вправо на количество разрядов, равное АР= РЛ – Р1 После каждого сдвига порядок увеличивается на единицу. После выравнивания порядков производится алгебраическое сложение мантисс по правилам для чисел с фиксированной точкой. Порядок результата принимается равным порядку большего слагаемого. В процессе суммирования чисел с разными знаками может нарушиться нормализация результирующей мантиссы. В этом случае полученный результат нормализуется путем сдвига мантиссы результата влево. После каждого сдвига влево порядок результата уменьшается на единицу.

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

Заметь, в тексте есть условия, когда младший разряд теряется (читай уменьшение) и когда увеличивается. Чтобы понять, бери и пробуй произвести операции побитно, согласно реализации сложений вещественных чисел. Уверен, найдешь эти самые "моменты".

JoeUser, спс, почитаю. Мдяяя, там не так все просто)
вообще, всякие представления вещ.чисел в памяти всегда какие-то заумные слишком), а всякие стандарты, аля IEEE-754 и т.п какие-то совсем лютые...

всем спс. за помощь

Скрытый текст

раз в 3 месяца тестирую себя по ЕГЭ по инфе (чтобы не забывать основы, так сказать) и набираю стабильно 99-100 баллов из 100, но это явно др.информатика

Цитата FasterHarder @ 05.12.19, 20:21

почему так происходит??

Потому что в компьютере используется двоичное представление чисел, а в нём точно представимы только числа вида n/2^m, то есть только те дроби, знаменатель которых является степенью 2 (есть ещё ограничение на величину n). Ты же в программе шаг задаёшь десятичной дробью (как и все), или задаёшь число шагов, обычно не дающее в результате такой дроби. Поэтому реальный шаг получается или чуть меньше, чем нужно, или чуть больше. В одном случае всё выглядит именно так, как ожидается (именно выглядит, на самом деле всё чуть хуже), в другом у тебя в конце появляется или не хватает одного значения.
Бороться с этим можно несколькими способами.
Один из них - отказаться от вещественного шага в пользу целочисленного счётчика, а нужное значение параметра цикла вычислять интерполяцией. Умножение на предварительно вычисленный шаг при этом не самый лучший способ.
Другой - немного поменять условие завершения, чтобы оно учитывало неточность представления шага - вместо (x = startvalue ; x <= endvalue; x += dx) писать (x = startvalue ; x <= endvalue + TOL; x += dx), выбрав TOL так, чтобы он не сильно превышал накопленную ошибку.

В чём плохость использования вещественного шага? Предположим тебе нужно пройтись по значениям от 0 до 10 с шагом 0.1. Ты добился, чтобы последнее значение было в районе 10, а не 9.9. Но на самом деле ни одно целое из этого интервала, кроме начального 0.0 в значения параметра не попадёт. Будут значения очень близкие к 1, 2, 3..., но не равные им. В то же время цикл for (i = 0; i <= 100; ++i, x = i/10.0) {...} пройдёт по всем нужным x максимально близко.

Цитата JoeUser @ 05.12.19, 20:49

Согласно информатике:

А вот это, на самом деле, никакого отношения к описанной проблеме не имеет.

Цитата amk @ 05.12.19, 21:54

Другой - немного поменять условие завершения, чтобы оно учитывало неточность представления шага - вместо (x = startvalue ; x <= endvalue; x += dx) писать (x = startvalue ; x <= endvalue + TOL; x += dx), выбрав TOL так, чтобы он не сильно превышал накопленную ошибку.

я знаком с этим подходом и в 1-й мессаге о нем упомянул, НО, как понимаю, подобрать этот ТОЛ не представляется возможным, т к кол-во разбиений n задается динамически и может быть, как очень малым, например, 1-2, так и большим, например, 10^10 или около того...но я не уверен в сказанном)

Цитата amk @ 05.12.19, 21:54

числа вида n/2m

т е это рациональные числа, у которых знаменатель может быть представлен через степень 2-ки

amk, а машинное епсилон тут разве не участвует? или ты про него как бы косвенно написал

Добавлено 05.12.19, 22:09

Цитата amk @ 05.12.19, 21:54

В то же время цикл for (i = 0; i <= 100; ++i, x = i/10.0) {...} пройдёт по всем нужным x максимально близко.

интересное примечание, хм...

т е вообще не будет НИ капли погрешности, если, например, в программе оперируем числами формы n/2^m? т е там можно что-угодно складывать, умножать и пр. и погрешность БУДЕТ ОТСУТСТВОВАТЬ ПОЛНОСТЬЮ, да?

Сообщение отредактировано: FasterHarder - 05.12.19, 22:09

Цитата amk @ 05.12.19, 21:54

А вот это, на самом деле, никакого отношения к описанной проблеме не имеет.

Т.е. описание особенностей сложения числе с плавающей никакого отношения к авторскому "x += dx" никакого отношения не имеет? Ну-ну ...

Цитата JoeUser @ 05.12.19, 22:21

Ну-ну ...

Эта особенность в описанном случае увеличивает погрешность весьма незначительно. А в случае точного представления шага вообще никак не влияет на погрешность. Проблема им енно в приближённости представления вещественных чисел.

Цитата FasterHarder @ 05.12.19, 22:03

НО, как понимаю, подобрать этот ТОЛ не представляется возможным, т к кол-во разбиений n задается динамически и может быть, как очень малым, например, 1-2, так и большим

Задай TOL = dx/8. Или dx/16. Или dx/256. Вовсе не обязательно его жёстко фиксировать, можно вычислять непосредственно пере использованием.
Проблемы возникнут, если отношение endval/dx слишком большое, скажем 10^10.

Цитата amk @ 05.12.19, 23:43

Эта особенность в описанном случае увеличивает погрешность весьма незначительно.

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

Цитата amk @ 05.12.19, 21:54

В то же время цикл for (i = 0; i <= 100; ++i, x = i/10.0) {...} пройдёт по всем нужным x максимально близко.

Частный случай предложенного мною варианта.
Основная проблема в арифметике с плавающей точкой – это накапливание погрешностей. Но прикол в том, что этим накапливанием в общем случае можно пренебречь. Если мы просто имеем сложную формулу с кучей операций, то среднестатистически погрешности её вычисления можно считать белым шумом, суммируемым в идеале к нулю. Идеал, естественно недостижим, однако практически всегда погрешность не будет превышать погрешности единственной операции с неким малым >1 коэффициентом. Чем больше операций, тем медленнее растёт этот >1 коэффициент, см. нормальное распределение.
Однако конкретно тут у нас не сложная формула с кучей операций, у нас формула с единственной операцией, рассчитывающей dx, результат которой, этот самый dx, входит в другую формулу, т.б. +=, применяемой многократно. И это уже не белый шум, это систематическая погрешность, которая с каждой итерацией не взаимонейтрализуется, а накапливается. Выход один: избавиться от систематической погрешности. JoeUser предложил ненакапливаемый вариант, т.е. у него абсолютная погрешность будет продолжать расти, но относительная будет оставаться постоянной. Я и amk предложили вообще отказаться от фиксированного dx, рассчитывая сразу результат по единственной формуле, использующей точное промежуточное значение. В перспективе это может дать постоянную абсолютную погрешность, но это не гарантируется. В примере amk, например, так и будет при min/max диапазона одного порядка, но не будет, если они разных порядков.

Добавлено 06.12.19, 03:25

Цитата JoeUser @ 05.12.19, 22:21

Т.е. описание особенностей сложения числе с плавающей никакого отношения к авторскому "x += dx" никакого отношения не имеет?

Имеет, но является частным случаем несколько иного эффекта – потери точности. Даже единственная операция может дать очень большую относительную погрешность. Вычитаем из e число 2,718282 и внезапно для 32-битного формата представления получаем относительную погрешность ~72%. Вот это реальный бич высококритичного ПО, где никак нельзя не учитывать того факта, что абсолютно верные мат.модели не могут вот так просто быть взятыми и реализованными один в один в коде.

Добавлено 06.12.19, 03:27
P.S. Нынешний наш проект по верификации такого высококритичного ПО очень радуется численным интегрированиям по времени с шагом 8мс процессов, длиной в часы.

Цитата Qraizer @ 06.12.19, 03:11

потери точности.

Так это и есть тот самый тупичок, об который бъются поезда и вагоны нецелочисленной арифметики!

Ну он как бы не единственный. Однако самое западло в том, что появляется он всегда из-за угла громким скриммером. Ещё не осознал, а штаны уже мокрые. А то и пахнут.

Цитата Qraizer @ 06.12.19, 03:45

Однако самое западло в том, что появляется он всегда из-за угла громким скриммером

Вот и говорю, если перефразировать принцип Парето, то "20% незаметных косяков с вещественными - привносят 80% всеобщего гемора с результатами". А все из-за накапливаемых ошибок по причине уменьшения точности от итерации к итерации.