1626435587-51311eae4652e8ad616b5bdef025cbb3 (844239), страница 9
Текст из файла (страница 9)
п. При решениифизических задач следует стремиться к соблюдению баланса междупогрешностями различных типов. Так, например, если математическаямодель описывает физическую систему с точностью поряка 1 %, врядли имеет смысл бороться за восемь значащих цифр в ответе.Наконец, для удобства представим в виде таблицы основные характеристики типов данных с одинарной, двойной и четверной точностью(float, double и long double20 в языке Си, или binary 32, 64 и 128в соответствии со стандартом IEEE 754). Для использования в программном коде перечисленных в табл.
1 и других аналогичных им параметров, предусмотрены константы FLT_EPSILON, FLT_MAX, FLT_MIN,DBL_EPSILON и т. п. (см. файл <float.h>).Таблица 1Размер в битахКол-во бит мантиссы, Кол-во бит показателя, = 2− (ULP)max показатель maxmax число 2max (2 − 2)min норм. число 21−maxbinary 32(float)322381,19 · 10−71273,40 · 10381,18 · 10−38binary 64(double)6452112,22 · 10−1610231,80 · 103082,23 · 10−308binary 128128112151,93 · 10−34163831,19 · 1049323,36 · 10−4932В качестве дополнительной литературы к данному разделу можнорекомендовать вторую главу замечательной монографии [A8]. Краткоерассмотрение вопросов точности вычислений можно найти во введенииучебного пособия [3].
Классическое рассмотрение ряда алгоритмов инетривиальных вопросов, связанных с вычислениями на ЭВМ, выполнено в [A9]. Для владеющих английским языком небезынтересно будетзнакомство с [A1, гл. 3 и 14].Упражнения1. Напишите на языке Си программу для нахождения значения машинного эпсилон, последовательно прибавляя к единице 2−1 , 2−2 ,20Разрядность longdoubleможет варьироваться на разных платформах.402−3 и т. д. Определите число разрядов в мантиссе, максимальныйи минимальный двоичный порядок чисел при вычислениях с одинарной и двойной точностью.2. Сравните друг с другом четыре машинных числа: 1, 1 + 2 , 1 + и1++ 2 , объясните результат.
То же для чисел 1++ 2 и 1+ 2 +.3. Какое наименьшее число с плавающей точкой вида 2− можноприбавить к числу = 1000, чтобы сумма была больше ?4. Напишите программу, умножающую float x = 1 на 2 до тех пор,пока результат не выйдет за пределы разрядной сетки. Сколькошагов цикла потребовалось для этого? Многократно деля x на2, найдите общее количество степеней двойки, представимых воfloat. Является ли полученное число «круглым» в двоичной системе счисления? Почему?5. Сколько всего битов используется для записи чисел типа float?Выведите на экран значение sizeof(float).
Вычтите знаковыйбит и биты для записи мантиссы (см. задачу 1). Сравните остатокс двоичным логарифмом от количества степеней двойки, полученных в упражнении 4, объясните результат.6. Напишите на языке Си функцию для сравнения друг с другомдвух чисел с плавающей точкой на приближённое равенство сточностью до машинного эпсилон.7.
Напишите на языке Си программу, содержащую следующие двестроки кода, объясните результат её работы:float x = 0.1;if(x == 0.1) puts("equal"); else puts("not equal");8. Аналогично для следующих строк кода:double x = exp(777), y = x - x;puts((y == 0) ? "equal" : "not equal");puts((y == y) ? "equal" : "not equal");9. Скомпилируйте и выполните программу на языке Си, содержащую следующий фрагмент кода, объясните результат:41double L = 1, dx = 0.1;int i = 0;while(L > 0) ++i, L -= dx;printf("A section of length L=1 was divided into %d ""sub-sections of length %f\n", i, dx);10. Что напечатает программа, содержащая следующие строки:float x0 = FLT_EPSILON, x1 = x0/2, x2 = x1*2;printf("x2-x0 = %e\n", x2-x0);11.
Что будет, если в предыдущем примере заменить FLT_EPSILON наFLT_MIN? В чем смысл указанных констант?12. Пусть () := (100, ), где (, ) — решение задачи Коши = − ,для обыкновенного дифференциального уравнения 0 ≤ ≤ 100, (0) = . Используя общее решение ОДУ () =1 + + · , вычислите значения (1) и (tg( 4 )).
Объясните полученные результаты.13. Определите в gnuplot функцию (, ) = sin( + 2 · 10 ) и постройте график (, = const) при фиксированном целочисленном значении , выполнив приведённые ниже две строки кода.Объясните результат:f(x,y)=sin(x+2*pi*10**y) #a**b означает "a в степени b"plot [0:2*pi] f(x,0), f(x,3), f(x,5), f(x,15)14. Проверьте, что интерполяционный полином второй степени () =− 12 2 + (1 + 32 ) − проходит через точки (1, 1), (2, 2) и (3, 3 − ).Решая квадратное уравнение () = 0 на калькуляторе илиЭВМ, вычислите ближайший к нулю корень при = 10−4 , 10−6 ,1.3 · 10−8 . Сколько значащих цифр содержит ответ? Как следуетпереписать выражение для повышения точности?15. Опишите алгоритм, позволяющий определить систему счисления,используемую калькулятором для выполнения вычислений (двоичная или десятичная).
Известно, что калькулятор работает счислами с плавающей точкой; ввод данных и отображение результатов производится в десятичной системе счисления.16. Вычислите значение B(100, 200), где B(, ) = Γ()Γ()/Γ( + ) —бета-функция Эйлера, Γ( + 1) = ! — гамма-функция Эйлера.4217. Пусть ≤ — два машинных числа с плавающей точкой. Всегдали проверка условий ≤ 12 ( + ) и 12 ( + ) ≤ на ЭВМ будетдавать логическую истину? (H.
Björk [A9]).18. Продемонстрируйте нарушение неравенства Коши для чисел сплавающей точкой — найдите числа и , для которых при вычислении на ЭВМ справедливо ( + )2 > 2(2 + 2 ) [A9].19. Каких чисел с плавающей точкой больше — положительных илиотрицательных? Бо́льших единицы или меньших единицы?3. Решение конечных уравненийРешение большинства физических задач приводит к возникновениюконечных уравнений, т. е. уравнений, которые могут быть представлены в форме () = 0. Если () является полиномом, такие уравненияназывают алгебраическими, в противном случае — трансцендентными.Зачастую точное решение конечных уравнений не может быть выражено через элементарные функции.
Примером может служить хрестоматийная задача из квантовой механики о поиске уровня энергии основного состояния в прямоугольной яме конечной глубины:{︂~2 ′′−0 , || ≤ , + ( () − ) = 0, () =−0, || > .2Волновая функция основного состояния чётная и имеет вид{︂22 cos , || < () =, 2 = 2 (0 + ), 2 = − 2 .− , || > ~~43Условие непрерывности решения () и его первой производной в точке = приводит к трансцендентному уравнению, точное решениекоторого не может быть представлено в элементарных функциях:√︂√︂−22 01(1 − ) =− 1, =ctg.(4)2~0Неразрешимые в элементарных функциях трансцендентные уравнения часто возникают при решении спектральных задач для дифференциальных операторов второго порядка в частных производных, при решении задачи Коши для обыкновенных дифференциальных уравненийи в ряде других случаев.
В данной главе будут рассмотрены нескольконаиболее употребительных методов численного решения таких уравнений.3.1. Метод деления отрезка пополамПожалуй, самым простым в реализации методом поиска корнейфункции одной вещественной переменной является,или деления отрезка пополам. Пусть на отрезке [, ] задана непрерывная функция : [, ] → ℛ, причём () · () ≤ 0.
Из курса анализа известно, что непрерывная функция принимает на отрезке все промежуточные значения и, следовательно, существует точка* ∈ [, ] : (* ) = 0. Наша задача состоит в нахождении такойточки (корня ()).Для решения поставленной задачи вычислим значение функции всредней точке отрезка, ( 12 ( + )).
Если () · ( 12 ( + )) ≤ 0, положим1 = , 1 = 21 ( + ), в противном случае возьмем 1 = 12 ( + ), 1 = .Шаг дихотомии выполнен. Мы свели задачу о поиске корня на отрезке [, ] к аналогичной задаче для отрезка [1 , 1 ], ширина которого вдва раза меньше. Повторяя процесс многократно, получим последовательность отрезков [ , ], причём − → 0 при → ∞. Предлагаем читателю самостоятельно доказать сходимость процесса к корню* , используя свойства непрерывности () и принцип математическойиндукции, а здесь ограничимся лишь несколькими замечаниями, касающимися использования данного метода на ЭВМ.На каждом шаге метода дихотомии можно было бы выполнятьпроверку равенства нулю значения функции в средней точке отрезка ( 12 (+)) с остановом итерационного процесса при достижении корня.Ввиду конечности множества машинных чисел вероятность «точного»попадания в ноль функции () хотя и не равна нулю, но тем не менее достаточно мала на большинстве шагов дихотомии.
Как следствие,44метод дихотомиивведение дополнительной проверки не приведет к повышению эффективности программы, однако может сделать программный код болеечитаемым и понятным для некоторых пользователей.В эффективных реализациях алгоритма функцию следует вычислять один раз на каждом шаге, сохраняя найденные ранее значения () и () в памяти ЭВМ. «Лишние» вычисления значений функции могут привести к существенному снижению скорости работы программы в случаях, когда вычисление функции требует большого количества операций.При проверке знака функции в середине отрезка, ( 12 (+)), следуетиметь в виду возможность выхода за пределы разрядной сетки, а такжепринципиальную возможность «точного» попадания в искомый кореньв процессе итераций.
Условие выхода из цикла итераций должно бытьосновано на сравнении ширины отрезка и требуемой точности (не превосходящей точность, обеспечиваемую используемым типом данных сплавающей точкой), в противном случае может произойти зацикливание программы.Максимальная абсолютная погрешность определения корня * равна 21 ( − ) и на каждой итерации уменьшается вдвое. Следовательно,если искомое значение * известно по порядку величины, каждый шагдихотомии даёт один дополнительный бит мантиссы. Это позволяетоценить необходимое число шагов для достижения требуемой точности:log(0 /),(5)≈log 2где 0 и — точность начального приближения и требуемая точностьнахождения корня соответственно.Метод деления отрезка пополам обладает рядом ограничений: онпозволяет находить только простые корни (а также корни, кратностькоторых нечётна) и не обобщается на функции нескольких переменных.Ниже будет дан обзор других методов, которые при определённых условиях могут обладать более высокой скоростью сходимости, что можетбыть существенно в случае, когда вычисление каждого значения ()занимает достаточно много времени, либо когда требуется решениебольшого количества уравнений (например, на каждом шаге численного интегрирования системы дифференциальных уравнений).
Вместес тем дихотомия обладает важным преимуществом: гарантированнойсходимостью к корню непрерывной функции () : () · () < 0.453.2. Метод простых итерацийПерепишем уравнение () = 0 в виде () = , где () := ()+.Для простоты будем считать вначале, что функция () гладкая иудовлетворяет условию |′ ()| ≤ , где < 1 — некоторая постоянная.Пусть нам также известно некоторое начальное приближение 0 к корню * . Организуем итерационный процесс вида+1 = ( ), = 0, 1, 2, . . .(6)аДанный процесс изображён на рис.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
