На главную Наши проекты:
Журнал   ·   Discuz!ML   ·   Wiki   ·   DRKB   ·   Помощь проекту
ПРАВИЛА FAQ Помощь Участники Календарь Избранное RSS
msm.ru
Модераторы: Qraizer, Hsilgos
  
> Табулирование функции на отрезке с дробным шагом (накапливается погрешность) , идеально чистый СИ, visual studio 2010
    Всем хай! Без раскачки сразу переходим к делу!

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

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

    ExpandedWrap disabled
      #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;
      }


    схватил лютую проблему в этой строке:
    ExpandedWrap disabled
      for(x = X_LEFT_BORDER; x <= (X_RIGHT_BORDER + EPS); x += dx)


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

    Посмотрел несколько прожек в инете, там все легко и просто, пишут аля
    ExpandedWrap disabled
      for(x = X_LEFT_BORDER; x <= X_RIGHT_BORDER; x += dx)
    , но на рис.видно, что результаты нестабильны в этом случае.

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

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

    Добавлено
    P.S. забыл дописать, что нужно именно через for)
    Сообщение отредактировано: FasterHarder -
      FasterHarder, может имеет смысл в цикле задавать просто перебор шагов?
      Ну а значение x вычислять как dx*номер_шага. Таким образом у тебя на любом
      из шагов будет одинаковая погрешность, равная погрешности операции умножения.
        Не ты первый, не ты последний. ;)
        Правильнее – заранее рассчитать количество итераций как целое число, использовать его для организации цикла, а очередную итерацию рассчитывать линейной интерполяцией между минимумом и максимумом по текущему номеру точки из общего их количества.

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


          я понял вроде, попробовал, вот так, да?
          ExpandedWrap disabled
                //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------------------------------");


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

          но, такой код немного не соот-ет духу табулирования функции) не-не, я оставлю такой вариант, если нельзя сделать так:
          ExpandedWrap disabled
            for(x = X_LEFT_BORDER; x <= X_RIGHT_BORDER; x += dx)


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

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

          Добавлено
          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 -
            Цитата FasterHarder @
            почему так происходит??

            Согласно информатике:
            Цитата
            Сложение чисел с плавающей точкой

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

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

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

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

              Скрытый текст
              раз в 3 месяца тестирую себя по ЕГЭ по инфе (чтобы не забывать основы, так сказать) и набираю стабильно 99-100 баллов из 100, но это явно др.информатика :no:
                Цитата FasterHarder @
                почему так происходит??
                Потому что в компьютере используется двоичное представление чисел, а в нём точно представимы только числа вида n/2m, то есть только те дроби, знаменатель которых является степенью 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 @
                Согласно информатике:
                А вот это, на самом деле, никакого отношения к описанной проблеме не имеет.
                  Цитата amk @
                  Другой - немного поменять условие завершения, чтобы оно учитывало неточность представления шага - вместо (x = startvalue ; x <= endvalue; x += dx) писать (x = startvalue ; x <= endvalue + TOL; x += dx), выбрав TOL так, чтобы он не сильно превышал накопленную ошибку.


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

                  Цитата amk @
                  числа вида n/2m

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

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

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

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

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

                    :blink: Т.е. описание особенностей сложения числе с плавающей никакого отношения к авторскому "x += dx" никакого отношения не имеет? Ну-ну ...
                      Цитата JoeUser @
                      Ну-ну ...
                      Эта особенность в описанном случае увеличивает погрешность весьма незначительно. А в случае точного представления шага вообще никак не влияет на погрешность. Проблема им енно в приближённости представления вещественных чисел.
                      Цитата FasterHarder @
                      НО, как понимаю, подобрать этот ТОЛ не представляется возможным, т к кол-во разбиений n задается динамически и может быть, как очень малым, например, 1-2, так и большим
                      Задай TOL = dx/8. Или dx/16. Или dx/256. Вовсе не обязательно его жёстко фиксировать, можно вычислять непосредственно пере использованием.
                      Проблемы возникнут, если отношение endval/dx слишком большое, скажем 10^10.
                        Цитата amk @
                        Эта особенность в описанном случае увеличивает погрешность весьма незначительно.

                        С ростом итераций "незначительно" плавно превращается в "писец". Именно поэтому я и предложил не манипулировать числами с плавающей путем последовательного сложения, а вычислять на каждой итерации нужное значение путем умножения дельты на номер итерации и его сложения с базой.
                          Цитата amk @
                          В то же время цикл for (i = 0; i <= 100; ++i, x = i/10.0) {...} пройдёт по всем нужным x максимально близко.
                          Частный случай предложенного мною варианта.
                          Основная проблема в арифметике с плавающей точкой – это накапливание погрешностей. Но прикол в том, что этим накапливанием в общем случае можно пренебречь. Если мы просто имеем сложную формулу с кучей операций, то среднестатистически погрешности её вычисления можно считать белым шумом, суммируемым в идеале к нулю. Идеал, естественно недостижим, однако практически всегда погрешность не будет превышать погрешности единственной операции с неким малым >1 коэффициентом. Чем больше операций, тем медленнее растёт этот >1 коэффициент, см. нормальное распределение.
                          Однако конкретно тут у нас не сложная формула с кучей операций, у нас формула с единственной операцией, рассчитывающей dx, результат которой, этот самый dx, входит в другую формулу, т.б. +=, применяемой многократно. И это уже не белый шум, это систематическая погрешность, которая с каждой итерацией не взаимонейтрализуется, а накапливается. Выход один: избавиться от систематической погрешности. JoeUser предложил ненакапливаемый вариант, т.е. у него абсолютная погрешность будет продолжать расти, но относительная будет оставаться постоянной. Я и amk предложили вообще отказаться от фиксированного dx, рассчитывая сразу результат по единственной формуле, использующей точное промежуточное значение. В перспективе это может дать постоянную абсолютную погрешность, но это не гарантируется. В примере amk, например, так и будет при min/max диапазона одного порядка, но не будет, если они разных порядков.

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

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

                            Так это и есть тот самый тупичок, об который бъются поезда и вагоны нецелочисленной арифметики! :)
                              Ну он как бы не единственный. Однако самое западло в том, что появляется он всегда из-за угла громким скриммером. Ещё не осознал, а штаны уже мокрые. А то и пахнут.
                                Цитата Qraizer @
                                Однако самое западло в том, что появляется он всегда из-за угла громким скриммером

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


                                Рейтинг@Mail.ru
                                [ Script execution time: 0,1090 ]   [ 18 queries used ]   [ Generated: 29.03.24, 01:19 GMT ]