Введение в вычислительную математику

Предмет вычислительной математики. Обусловленность задачи, устойчивость алгоритма, погрешности вычислений. Задача численного дифференцирования

Показывать лекцию целиком

Предмет вычислительной математики. Обусловленность задачи, устойчивость алгоритма, погрешности вычислений. Задача численного дифференцирования

Первое применение вычислительных методов принадлежит древним египтянам, которые умели вычислять диагональ квадрата за конечное количество действий. Они также могли находить квадратный корень из 2, скорее всего, с помощью алгоритма, в дальнейшем получившего название формулы Герона, а еще позднее — метода Ньютона:

$$$ u_{k + 1} = \frac{1}{2}\left(u_k + \frac{2}{u_k}\right),\quad u_0 = a. $$$

С именем среднеазиатского врача, философа и математика Аль-Хорезми связано понятие алгоритма. Разработкой вычислительных методов занимались Л.Эйлер, которому принадлежит, по-видимому, первый численный метод для решения обыкновенных дифференциальных уравнений, И.Ньютон, О.Л.Коши, Ж.Л.Лагранж, А.М.Лежандр, П.С.Лаплас, А.Пуанкаре, П.Л.Чебышёв и многие другие известные математики. Решающую роль в развитии вычислительной математики как самостоятельной науки сыграли немецкий математик Карл Рунге и русский математик, механик и кораблестроитель А.Н.Крылов.

В наше время вычислительная математика получила значительный импульс в 1950-е годы, что было связано с развитием ядерной физики, механики полета, аэродинамики спускаемых космических аппаратов. В дальнейшем решались задачи, связанные не только с расчетами действия ядерного взрыва и обтеканием боеголовок стратегических ракет. Численные методы нашли свое применение в таких областях как динамика атмосферы, термогидрография, физика плазмы, механика горных пород и ледников, синергетика, биомеханика, теория оптимизации, математическая экономика и др. Наиболее наукоемки и требуют максимальных вычислительных ресурсов задачи физики, механики и электродинамики сплошных сред. К ним относятся системы уравнений в частных производных Эйлера, Лагранжа, Максвелла и др., кинетической теории газов (уравнения Власова, Берда), а также задачи многомерной оптимизации. Развитие многих вычислительных методов — заслуга ученых, неразрывно связанных с МФТИ. Среди них следует упомянуть академиков А.А. Дородницына, О.М.Белоцерковского, А.А.Самарского, членов-корреспондентов РАН А.С.Холодова, Б.Н.Четверушкина, Ю.П.Попова, С.П.Курдюмова, профессоров В.С.Рябенького, Э.Э.Шноля, Р.П.Федоренко, Л.А.Чудова, В.Ф.Дьяченко, И.М.Гельфанда, Г.А.Тирского, А.П.Фаворского.

Вычислительная математика отличается от других математических дисциплин и обладает специфическими особенностями.

  • Вычислительная математика имеет дело не только с непрерывными, но и с дискретными объектами. Так, вместо отрезка прямой часто рассматривается система точек $$\left\{t_k\right\}_{k = 0}^K$$, вместо непрерывной функции f(x) — табличная функция $$\left\{f_k\right\}_{k = 0}^K$$, вместо первой производной — ее разностная аппроксимация, например,

    $$\frac{f_{k + 1} - f_k}{x_{k + 1} - x_k},$$ $$k = 0 \div K,\quad x_{k+1} > x_k.$$

    Такие замены, естественно, порождают погрешности метода.

  • В машинных вычислениях присутствуют числа с ограниченным количеством знаков после запятой из-за конечности длины мантиссы при представлении действительного числа в памяти ЭВМ. Другими словами, в вычислениях присутствует машинная погрешность (округления) $$\delta _M.$$ Это приводит к вычислительным эффектам, неизвестным, например, в классической теории обыкновенных дифференциальных уравнений, уравнений математической физики или в математическом анализе.
  • В вычислительной практике большое значение имеет обусловленность задачи, т.е. чувствительность ее решения к малым изменениям входных данных.
  • В отличие от "классической" математики выбор вычислительного алгоритма влияет на результаты вычислений.
  • Существенная черта численного метода — экономичность вычислительного алгоритма, т.е. минимизация числа элементарных операций при выполнении его на ЭВМ.
  • Погрешности при численном решении задач делятся на две категории — неустранимые и устранимые. К первым относят погрешности, связанные с построением математической модели объекта и приближенным заданием входных данных, ко вторым — погрешности метода решения задачи и ошибки округления, которые являются источниками малых возмущений, вносимых в решение задачи.
  • Специфику машинных вычислений можно пояснить на нескольких элементарных примерах.

    1.1. Обусловленность задачи

    Пример 1.1. Вычислить все корни уравнения

    $$x^4 - 4x^3 + 8x^2 - 16x + 15.\underbrace{99999999}_8 = {(x - 2)}^4 - 10^{- 8} = 0.$$

    Точное решение задачи легко найти:

    $$(x - 2)^{2} = \pm 10^{- 4}, \\ x_{1}= 2,01; x_{2}= 1,99; x_{3,4}= 2 \pm 0,01i.$$

    Если компьютер работает при $$\delta _M > 10^{ - 8}$$, то свободный член в исходном уравнении будет округлен до 16,0 и, с точки зрения представления чисел с плавающей точкой, будет решаться уравнение (x-2)4= 0, т.е. x1,2,3,4 = 2, что, очевидно, неверно. В данном случае малые погрешности в задании свободного члена $$\approx 10^{-8}$$ привели, независимо от метода решения, к погрешности в решении $$\approx 10^{-2}.$$

    Пример 1.2. Решается задача Коши для обыкновенного дифференциального уравнения 2-го порядка:

    u''(t) = u(t), u(0) = 1, u'(0) = - 1.

    Общее решение имеет вид

    u(t) = 0,5[u(0) + u'(0)]et  + 0,5[u(0) - u'(0)]e- t.

    При заданных начальных данных точное решение задачи: u(x) = e-t, однако малая погрешность $$\delta$$ в их задании приведет к появлению члена $$\delta e^t$$, который при больших значениях аргумента может существенно исказить решение.

    Пример 1.3. Пусть необходимо найти решение обыкновенного дифференциального уравнения

    $$u^{\prime} = 10u,\quad u = u(t),\\ u(t_0) = u_0,\quad t \in [0,1].$$

    Его решение: $$u(t) = u_0e^{10(t - t_0 )}$$, однако значение u(t0) известно лишь приближенно: $$u(t_0) \approx u_0^*$$, и на самом деле $$u^*(t) = u_0^*e^{10(t - t_0)}.$$

    Соответственно, разность u* - u будет

    $$u^* - u = (u_0^* - u_0)e^{10(t - t_0)}.$$

    Предположим, что необходимо гарантировать некоторую заданную точность вычислений $$\varepsilon > 0$$ всюду на отрезке $$t \in [0,1].$$ Тогда должно выполняться условие

    $$|{u^*(t) - u(t)}| \le \varepsilon.$$

    Очевидно, что$$\max\limits_{t \in [0,1]} |{u^*(t) - u(t)}| = |{u*(1) - u(1)}| = |{u_0^* - u_0}|e^{10(1 - t_0)}.$$

    Отсюда можно получить требования к точности задания начальных данных $$\delta\colon |{u_0^* - u_0}| < \delta , \delta \le \varepsilon e^{ - 10}$$ при t0= 0.

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

    Решение оказывается очень чувствительным к заданию начальных данных. Такого рода задачи называются плохо обусловленными.

    Пример 1.4. Решением системы линейных алгебраических уравнений (СЛАУ)

    $$\left\{ \begin{array}{l} u + 10v = 11 \\ 100u + 1001v = 1101 \\ \end{array} \right.$$

    является пара чисел {1, 1}.

    Изменив правую часть системы на 0,01, получим возмущенную систему

    $$\left\{ \begin{array}{l} u + 10v = 11.01 \\ 100u + 1001v = 1101 \\ \end{array} \right.$$

    с решением {11.01; 0.00}, сильно отличающимся от решения невозмущенной системы. Эта система также плохо обусловлена.

    Пример 1.5. Рассмотрим полином

    (x - 1)(x - 2)...(x - 20)=x20 - 210x19 + ...,

    корни которого x1 = 1, x2 = 2, …, x20 = 20.

    Положим, что коэффициент (-210) при x19 увеличен на $$\approx 10^{-7}.$$ В результате вычислений с 11-ю значащими цифрами получим совершенно иные корни: $$x_{1} = 1,00; x_{2} = 2,00; x_{3} = 3,00; x_{4} = 4,00; x_{5} = 5,00; x_{6} = 6,00; x_{7} = 7,00; x_{8} = 8,01; x_{9} = 8,92; x_{10,11} = 10,1 \pm 0,644i; x_{12,13} = 11,8 \pm 1,65i; x_{14,15} = 14,0 \pm 2,52i; x_{16,17} = 16,7 \pm 2,81i; x_{18,19} = 19,5 \pm 19,4i; x_{20} = 20,8.$$

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

    1.2. Влияние выбора вычислительного алгоритма на результаты вычислений

    Пример 1.6. Пусть необходимо вычислить значение выражения $$$ \left({\frac{{\sqrt{2} - 1}}{{\sqrt{2} + 1}}}\right)^3 . $$$

    Избавившись от знаменателя, получаем $$(\sqrt{2} - 1)^6 = (3 - 2\sqrt{2})^3 = 99 - 70\sqrt{2}.$$

    Полагая

    а)

    $$$ \sqrt{2} \approx \frac{7}{5} = 1,4 , $$$

    в)

    $$$ \sqrt{2} \approx \frac{17}{12} = 1,41(6) $$$

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

    $$\sqrt{2}$$ $${(\sqrt{2} - 1)}^6$$ $${(3 - 2\sqrt{2})}^3$$ $$99 - 70\sqrt{2}$$
    7/5 0,004096 0,008000 1
    17/12 0,005233 0,004630 - 0,1(6)

    Очевидно, что столь значительное различие в результатах вызвано влиянием ошибки округления в задании $$\sqrt{2}.$$

    Пример 1.7. Вычисление функции sin x с помощью ряда Тейлора.

    Из курса математического анализа известно, что функция синус представляется своим рядом Тейлора

    $$$ \sin x = x - \frac{x^3}{3!} + \frac{x^5}{5!} - \frac{x^7}{7!} + \ldots , $$$

    причем радиус сходимости ряда равен бесконечности — ряд сходится при любых значениях x.

    Вычислим значения синуса при двух значениях аргумента. Пусть сначала $$x = 0,5236 (30^\circ).$$ Будем учитывать лишь члены ряда, большие, чем 10- 4. Выполнив вычисления с четырьмя значащими цифрами, получим sin (0.5236) = 0.5000, что соответствует принятой точности.

    Пусть теперь $$x = 25,66 (1470^\circ).$$ Если вычисления по данной формуле проводить с восемью значащими цифрами, то получим абсурдный результат: $$sin (25,66) \approx 24$$ (учитывались члены ряда, большие, чем 10-8 ).

    Разумеется, выходом из создавшейся ситуации может быть использование формул приведения.

    Пример 1.8. Вычисление функции ex с помощью ряда Тейлора.

    Из курса математического анализа известно, что экспонента представляется своим рядом Тейлора

    $$$ e^x = 1 + x + \frac{x^2}{2!} + \frac{x^3}{3!} + \ldots , $$$

    радиус сходимости этого ряда также равен бесконечности.

    Приведем некоторые результаты расчетов ( $$e^x_M$$ — значения экспоненты, вычисленные на компьютере).

    x ex exM
    1 2,718282 2,718282
    20 $$4,8516520 \cdot 10^8$$ $$4,8516531 \cdot 10^8$$
    -1 0,3678795 0,3678794
    -10 $$4,5399930 \cdot10 ^{-5}$$ $$-1,6408609 \cdot 10^{-4}$$
    -20 $$2,0611537 \cdot 10^{-9}$$ 1,202966

    Выходом из этой ситуации может быть использование для отрицательных аргументов экспоненты формулы

    $$$ e^{- x} = \frac{1}{e^x } = \frac{1}{1 + x + \frac{x^2}{2!} + \ldots } . $$$

    Естественно ожидать рост ошибок округления при вычислении рассматриваемой функции при больших значениях аргумента x. В этом случае можно использовать формулу ex = en + a = enea, где n = [x].

    Пример 1.9. Рассмотрим следующий метод вычисления интеграла

    $$I_n = \int\limits_0^1 {x^n}{e^{1-x}dx}, n=0,1,2\ldots$$

    Интегрирование по частям дает

    $$I_n = \int\limits_0^1 {x^n}{d(-e^{1-x})} = - x^{n-1}\cdot e^{1-x} |\limits_0^1 + \int\limits_0^1 e^{1-x}\cdot d(x^n) = -1 + \int\limits_0^1 nx^{n-1}e^{1-x}dx,$$

    откуда следует$$I_0 = \int\limits_0^1 e^{1 - x}dx = e - 1 \approx 1,71828.$$ $$I_n = nI_{n-1}-1, n \ge 1.$$

    Тогда

    $$\begin{gather*} I_1 = 1 \cdot I_0 - 1 \approx 0.71828, \quad I_2 = 2I_1 - 1 \approx 0.43656, \quad I_3 = 3I_2 - 1 \approx 0.30968, \\ I_4 = 4I_3 - 1 \approx 0.23872, \quad I_5 = 5I_4 - 1 \approx 0.1936, \quad I_6 = 6 \cdot I_5 - 1 \approx 0.16160, \\ I_7 = 7I_6 - 1 \approx 0.13120, \quad I_8 = 8I_7 - 1 \approx 0.00496, \quad I_9 = 9I_8 - 1 \approx - 0.55360, \\ {I}_{10} = 10\cdot I_9 - 1 \approx -6.5360 \end{gather*} $$

    Очевидно, что отрицательные значения при n = 9,10 не имеют смысла. Дело в том, что ошибка, сделанная при округлении I0 до 6-ти значащих цифр сохранилась при вычислении I1, умножилась на 2! при вычислении I2, на 3! — при вычислении I3, и так далее, т.е. ошибка растет очень быстро, пропорционально n!.

    Пример 1.10. Рассмотрим методический пример вычислений на модельном компьютере, обеспечивающем точность $$\delta_M = 0,0005.$$ Проанализируем причину происхождения ошибки, например, при вычитании двух чисел, взятых с точностью до третьей цифры после десятичной точки u = 1,001, v = 1,002, разность которых составляет $$\Delta = |v_{M} - u_{M}| = 0,001.$$

    В памяти машины эти же числа представляются в виде

    $$u_M = u(1 + \delta_M^u), v_M = v(1 + \delta_M^v), \mbox{ причем } \mid \delta_M^u\mid \le \delta_M \mbox{ и } \mid \delta_M^v\mid \le \delta_M.$$

    Тогда $$u_M - u \approx u\delta_M^u, v_M - v \approx v\delta_M^v.$$

    Относительная ошибка при вычислении разности uM - vM будет равна

    $$$ \delta = \frac{(u_M - v_M) - (u - v)}{(u - v)} = \frac{(u_M - u) - (v_M - v)}{(u - v)} = \frac{u\delta_M^u - v\delta_M^v}{(u - v)}. $$$

    Очевидно, что $$$ \delta = \left|{\frac{\delta_M^u - \delta_M^v}{\Delta }} \right| \le \frac{2\delta_M}{0,001} \approx 2000\delta_M = 1, $$$ т.е. все значащие цифры могут оказаться неверными.

    Пример 1.11. Рассмотрим рекуррентное соотношение ui+1 = qui, $$i \ge 0$$, u0 = a, q > 0, ui > 0.

    Пусть при выполнении реальных вычислений с конечной длиной мантиссы на i -м шаге возникла погрешность округления, и вычисления проводятся с возмущенным значением $$u_i^M = u_i + \delta_i$$, тогда вместо ui+1 получим $$u_{i + 1}^M = q(u_i + \delta_i) = u_{i + 1} + q\delta_i$$, т.е. $$\delta_{i + 1} = q\delta_i,\quad i = 0,1,\ldots.$$

    Следовательно, если |q| > 1, то в процессе вычислений погрешность, связанная с возникшей ошибкой округления, будет возрастать ( алгоритм неустойчив ). В случае $$\mid q\mid \le 1$$ погрешность не возрастает и численный алгоритм устойчив.

    1.3. Экономичность вычислительного метода

    Пример 1.12. Пусть требуется вычислить сумму S = 1 + x + x2 + … + x1023 при 0 < x < 1. Для последовательного вычисления $$x, x^2 = x \cdot x, \ldots, x^{1023} = x^{1022} \cdot x$$ необходимо проделать 1022 умножения, а затем столько же сложений.

    Однако если заметить, что $$$ S = \frac{1 - x^{1024}}{1 - x} $$$, то количество арифметических действий значительно уменьшается; в частности, для вычисления x1024 требуется всего 10 умножений: $$x^2 = x \cdot x; x^4 = (x)^2 (x)^2, \ldots , x^{1024} = (x)^{512}(x)^{512}.$$

    Пример 1.13. Вычисления значений многочленов. Если вычислять значение многочлена P(x) = a0 + a1x + a2x2 + … + anxn "в лоб", т.е. вычислять значения каждого члена и суммировать, то окажется, что необходимо выполнить (n2 + [n/2]) умножений и n сложений. Кроме того, такой способ вычислений может привести к накоплению ошибок округления при вычислениях с плавающей точкой.

    Его очевидным улучшением является вычисление каждого члена последовательным умножением на x. Такой алгоритм требует (2n - 1) умножение и n сложений.

    Еще более экономичным алгоритмом является хорошо известная в алгебре схема Горнера:

    P(x) = ((...((anx + an - 1)x + an - 2)x + ... + a0),

    требующая n операций сложения и n операций умножения. Этот метод был известен в средние века в Китае под названием Тянь-Юань и был заново открыт в Европе в начале XIX века англичанином Горнером и итальянцем Руффини.

    Пример 1.14. Рассмотрим систему линейных алгебраических уравнений (СЛАУ) вида $${\mathbf{Au}} = {\mathbf{f}}, {\mathbf{u}} = \{ u_1 , \ldots , u_n \}^T, {\mathbf{f}} = \{ f_1 , \ldots ,f_n \}^T$$, с трехдиагональной матрицей

    $$\mathbf{A} = \left( \begin {array}{cccccc} b_1 c_1 0 \ldots 0 \\ a_2 b_2 c_2 \ldots 0 \\ 0 a_3 b_3 c_3 \ldots 0 \\ 0 \ldots 0 a_n b_n \\ \end {array} \right).$$

    Если проводить решение такой системы без учета специфической структуры матрицы (например, с помощью метода Гаусса), то количество арифметических действий будет порядка n3, если же учесть эту структуру, то количество операций можно уменьшить до n.

    1.4. Погрешность метода

    Оценим погрешность при вычислении первой производной при помощи соотношения: $$$ f^{\prime}(x) \approx \frac{f(x + h) - f(x)}{h} $$$:

    $$$ \frac {f(x + h ) - f(x)}{h} = \frac {[f(x) + hf^{\prime}(x) + O(h^2)] - f(x)}{h} = f^{\prime}(x) + O(h), $$$

    где O(h) есть погрешность метода. В данном случае под погрешностью метода понимается абсолютная величина разности $$$ \left|{f^{\prime}(x) - \frac{f(x + h) - f(x)}{h}}\right| $$$, которая составляет O(h) (более точно $$$ \frac{h}{2}f^{\prime\prime}(\xi ) $$$, где $$\xi \in [x,x + h]$$ ).

    Если же взять другой метод вычисления производной $$$ f^{\prime}(x) \approx \frac{f(x + h) - f(x - h)}{2h}$$$, то получим, что его погрешность составляет O(h2), это оказывается существенным при малых h. Однако уменьшать h до бесконечности не имеет смысла, что видно из следующего примера. Реальная погрешность при вычислении первой производной будет

    $$$ \Delta = \frac{h}{2} \max\limits_{\xi \in [x,x + h]} \left|{f^{\prime\prime}(\xi )}\right| + \frac{2\delta_M}{h} = O (h) + O(h^{- 1}), $$$

    поскольку абсолютная погрешность вычисления значения функции за счет машинного округления не превосходит $$$ \frac{2\delta_M}{h} $.$$

    В этом случае можно найти оптимальный шаг h. Будем считать полную погрешность в вычислении производной $$\Delta$$ функцией шага h. Отыщем минимум этой функции. Приравняв производную $$\Delta '_{h}(h)$$ к нулю, получим оптимальный шаг численного дифференцирования

    $$$ h_{\text{опт}} = 2\sqrt{\frac{\delta_M}{\max\limits_{\xi \in [x,x + h]} \left|{f^{\prime\prime}(\xi )}\right|}}. $$$

    Выбирать значение h меньше оптимального не имеет смысла, так как при дальнейшем уменьшении шага суммарная погрешность начинает расти из-за возрастания вклада ошибок округления.

    1.5. Элементы теории погрешностей

    Определение. Пусть u и u* — точное и приближенное значение некоторой величины, соответственно. Тогда абсолютной погрешностью приближения u* называется величина $$\Delta (u^{*})$$, удовлетворяющая неравенству

    $$\left|{u - u^*}\right| \le \Delta (u^*).$$

    Определение. Относительной погрешностью называется величина $$\delta (u^*)$$, удовлетворяющая неравенству

    $$$ \left|{\frac{u - u^*}{u^*}}\right| \le \delta (u^*). $$$

    Обычно используется запись $$u = u^*(1 \pm \delta (u^*)).$$

    Определение. Пусть искомая величина u является функцией параметров $$t_1, \ldots t_n \in \Omega$$, u* — приближенное значение u. Тогда предельной абсолютной погрешностью называется величина

    $$D(u^*) = \sup\limits_{(t_1, \ldots ,t_n) \in \Omega } \left|{u(t_1, \ldots ,t_n) - u^*}\right|,$$

    Предельной относительной погрешностью называется величина D(u*)/| u*|.

    Пусть $$\left|{t_j - t_j^*}\right| \le \Delta (t_j^* ), j = 1 \div n$$ — приближенное значение $$u* = u(t_1^*, \ldots ,t_n^* ).$$ Предполагаем, что u - непрерывно дифференцируемая функция своих аргументов. Тогда, по формуле Лагранжа,

    $$u(t_1, \ldots ,t_n) - u^* = \sum\limits_{j = 1}^n \gamma_j (\alpha )(t_j - t_j^*),$$

    где $$\gamma_j (\alpha ) = u'_{t_j}(t_1^* + \alpha (t_1 - t_1^*), \ldots ,t_n^* + \alpha (t_n - t_n^*)), 0 \le \alpha \le 1.$$

    Отсюда$$\left|{u(t_1, \ldots ,t_n) - u^*}\right| \le D_1 (u^*) = \sum\limits_{j = 1}^n b_j \Delta (t_j^*),$$ где $$b_j = \sup\limits_\Omega \left|{u^{\prime}_{t_j}(t_1, \ldots ,t_n)}\right|.$$

    Можно показать, что при малых $$\rho = \sqrt{{(\Delta (t_1^* ))}^2 + \ldots + {(\Delta (t_n^* ))}^2 }$$ эта оценка не может быть существенно улучшена. На практике иногда пользуются грубой (линейной) оценкой

    $$\left|{u(t_1, \ldots ,t_n) - u^*}\right| \le D_2 (u^*), \mbox{ где }D_2 (u^*) = \sum\limits_{j = 1} \left|{\gamma_j (0)}\right| \Delta (t^*).$$

    Несложно показать, что

    a) $$\Delta ( \pm t_1^* \pm \ldots \pm t_n^*) = \Delta (t_1^* ) + \ldots + \Delta (t_n^* ) $$, предельная погрешность суммы или разности равна сумме предельных погрешностей.

    b) Предельная относительная погрешность произведения или частного приближенного равна сумме предельных относительных погрешностей

    $$\delta (t_1^* \cdots t_m^* \cdot d_1^{* - 1} \cdots d_m^{* - 1} ) = \delta (t_1^* ) + \ldots + \delta (t_m^*) + \delta (d_1^*) + \ldots + \delta (d_n^*).$$

    1.6. Задача численного дифференцирования

    В пункте 1.2 уже была введена простейшая формула численного дифференцирования. Рассмотрим задачу приближенного вычисления значения производной подробнее.

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

    x0 = a, 
    xj = a + jh,, j = 0, 1, ... , N

    Пусть f(xj) = fj — функция, определенная в узлах сетки. Такие функции будут называться табличными, или сеточными функциями. Считаем, кроме того, что рассматриваемая сеточная функция есть проекция (или ограничение) на сетку некоторой гладкой нужное число раз непрерывно дифференцируемой функции f(x). По определению производной

    $$$ f^{\prime}(x) = \lim\limits_{\Delta x \to 0}\frac{f(x + \Delta x) - f(x)}{\Delta x}, $$$

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

    $$$ f^{\prime}(x_j) \approx \frac{f(x_j + h) - f(x_j)}{h}. $$$

    Если параметр h достаточно мал, то можно считать полученное значение производной достаточно точным. Погрешность формулы (1.1) оценена в пункте 1.2. Как показано выше, при уменьшении шага сетки h ошибка будет уменьшаться, но при некотором значении h ошибка может возрасти до бесконечности. При оценке погрешности метода обычно считается, что все вычисления были точными. Но существует ошибка округления. При оценке ее большую роль играет машинный $$\varepsilon$$ — мера относительной погрешности машинного округления, возникающей из-за конечной разрядности мантиссы при работе с числами в формате с плавающей точкой. Напомним, что по определению машинным $$\varepsilon$$ называют наибольшее из чисел, для которых в рамках используемой системы вычислений выполнено $$1 + \varepsilon = 1.$$ Тогда абсолютная погрешность при вычислении значения функции (или представлении табличной функции) есть $$f(x_j) \cdot \varepsilon .$$ Максимальный вклад погрешностей округления при вычислении производной по формуле (1.1) будет$$\frac{2f(x_j) \cdot \varepsilon }{h},$$ тогда, когда члены в знаменателе (1.1) имеют ошибки разных знаков.

    Пусть k = max |f(x)|, максимум ищется на отрезке, на котором вычисляются значения производных. Тогда суммарная ошибка, состоящая из погрешности метода и погрешности округления, есть $$$ \Delta = \frac{2k\varepsilon }{h} + \frac{M_2 h}{2},\quad M_2 = \max \left|{f^{\prime\prime}(x)}\right|.$$$

    Для вычисления оптимального шага численного дифференцирования найдем минимум суммарной ошибки, как функции шага сетки $$$\frac{M_2}{2} - \frac{2k\varepsilon }{h^2} = 0 $$$, откуда

    $$$h_{opt} = 2\sqrt{\frac{k\varepsilon }{M_2}} . $$$

    Если требуется повысить точность вычисления производных, необходимо воспользоваться формулами, имеющими меньшие погрешности метода. Так, из курсов математического анализа известно, что

    $$$ f^{\prime}(x) = \lim\limits_{\Delta x \to 0} \frac{f(x + \Delta x) - f(x - \Delta x)}{2\Delta x}.$$$

    По аналогии напишем конечно-разностную формулу

    $$$ f^{\prime}(x_j) \approx \frac{f(x_j + h) - f(x_j - h)}{2h}. $$$

    (1.2) — формула с центральной разностью. Исследуем ее на аппроксимацию, т.е. оценим погрешность метода. Предположим, что функция, которую спроектировали на сетку, трижды непрерывно дифференцируема, тогда

    $$$ f(x_j + h) = f(x_j) + f^{\prime}(x_j)h + f^{\prime\prime}(x_j) \frac{h^2}{2} + f^{\prime\prime\prime}(x_j + \theta_1 h) \frac{h^3}{6}, \\ f(x_j - h) = f(x_j) - f^{\prime}(x_j)h + f^{\prime\prime}(x_j) \frac{h^2}{2} - f^{\prime\prime\prime}(x_j - \theta_2 h) \frac{h^3}{6}.$$$

    Погрешность метода определяется 3-й производной функции. Введем $$M_3 = \max\limits_{x \in [a,b]} \left|{f^{\prime\prime\prime}(x)}\right|$$, тогда суммарная погрешность при вычислении по формуле с центральной разностью есть

    $$$ \Delta = \frac{M_3 h^2 }{3} + \frac{k\varepsilon}{h},$$$

    для вычисления оптимального шага, находя минимум погрешности, как функции шага сетки, имеем $$$ \frac{2M_3 h}{6} - \frac{k\varepsilon }{h^2} = 0 $$$, откуда $$$ h_{opt} = \sqrt[3]{\frac{3k\varepsilon }{M_3}}$.$$

    Для более точного вычисления производной необходимо использовать разложение более высокого порядка, шаг hopt будет увеличиваться.

    Формула (1.1) — двухточечная, (1.2) — трехточечная: при вычислении производной используются точки (узлы) xj (узел входит с нулевым коэффициентом), xj + h, xj - h — совокупность узлов, участвующих в каждом вычислении производной, в дальнейшем будем иногда называть сеточным шаблоном.

    Введем на рассматриваемом отрезке шаблон из нескольких точек.

    Считаем, что сетка равномерная — шаг сетки постоянный, расстояния между любыми двумя соседними узлами равны. Используем для вычисления значения первой производной следующую приближенную (конечно-разностную) формулу:

    $$$ f^{\prime}(x_j) \approx \frac{1}{h}\sum\limits_{k = - l}^{m} \alpha _k f(x_j + kh), $$$

    шаблон включает l точек слева от рассматриваемой точки xj и m справа. Коэффициенты $$\alpha_k$$ — неопределенные коэффициенты. Формула дифференцирования может быть и односторонней — либо l, либо m могут равняться нулю. В первом случае иногда называют (на наш взгляд, не слишком удачно) такую приближенную формулу формулой дифференцирования вперед, во втором — формулой дифференцирования назад. Потребуем, чтобы (1.3) приближала первую производную с точностью O(hl + m) . Используем разложения в ряд Тейлора в окрестности точки xj. Подставляя их в (1.3), получим

    $$$ \frac{1}{h}\sum\limits_{k = - l}^{m} \alpha_k f(x_j + kh) = \frac{1}{h}f(x_j)\sum \alpha_k + f^{\prime}(x_j)\sum k\alpha_k + f^{\prime\prime}(x_j)\sum \frac{k^2}{2}\alpha_k h + \\ + f^{\prime\prime\prime}(x_j)\sum \frac{k^3}{6}\alpha_k h^2 + \ldots + f^{\left( n\right)}(x_j) \sum \alpha_k \frac{k^n}{n!}h^{n - 1} + \ldots $$$

    Потребуем выполнение условий:

    $$$ \sum \alpha_k = 0, \sum k\alpha_k = 1, \sum \alpha_k \frac{k^2}{2} = 0, \ldots \sum \alpha_k \frac{k^n}{n!} = 0, \ldots $$$

    Получаем систему линейных алгебраических уравнений для неопределенных коэффициентов $$\alpha$$ (1.4). Матрица этой системы есть

    $$\left( \begin{array}{cccc} 1 1 \ldots 1 \\ - l - l + 1 \ldots m \\ l^2 {(l - 1)}^2 \ldots m^2 \\ {(- l)}^3 {(- l + 1)}^3 \ldots m^3 \\ \ldots \ldots \ldots \ldots \\ \end{array} \right) .$$

    Вектор правых частей (0, 1, 0, ..., 0)T.

    Определитель данной матрицы — детерминант Вандермонда. Из курса линейной алгебры следует, что он не равен нулю. Тогда существует единственный набор коэффициентов $$\alpha$$, который позволяет найти на шаблоне из (1 + l + m) точек значение первой производной с точностью O(hl + m).

    Для нахождения второй производной можно использовать ту же самую формулу (1.3) с небольшой модификацией

    $$$ f^{\prime}(x_j) \approx \frac{1}{h^2}\sum\limits_{k = - l}^{m} \alpha_k f(x_j + kh) , $$$

    только теперь $$$ \sum k\alpha_k = 0, \sum \alpha_k k^2 = \frac{1}{2} $.$$

    Очевидно, что и данная система уравнений для нахождения неопределенных коэффициентов имеет единственное решение. Для получения с той же точностью приближенных значений производных до порядка l + m включительно с точностью O(hl + m) модификации формулы (1.3) и условий (1.4) очевидны, набор неопределенных коэффициентов находится единственным образом.

    Таким образом, доказано следующее утверждение. На сеточном шаблоне, включающем в себя N + 1 точку, с помощью метода неопределенных коэффициентов всегда можно построить единственную формулу для вычисления производной от первого до n порядка включительно с точностью O(hN).

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

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

    Читателям предлагается оценить значение оптимального шага при вычислениях по формулам типа (1.3) самостоятельно.

    1.7. Задачи

  • Найти абсолютную предельную погрешность, погрешность по производной, линейную погрешность для функции u = t10, если заданы точка приближения t* = 1, значение функции u* в этой точке и погрешность $$\Delta t^{*} = 10^{ -1}.$$

    Решение : обозначим$$b=\sup\limits_{\mid{t - 1} \mid \le 0,1} \left|{u^{\prime}_t (t)}\right| = \sup\limits_{\mid{t - 1}\mid \le 0.1} 10 \cdot t^9 = 10 \cdot {(1,1)}^9 \approx \\ \approx 23, \ldots$$

    Абсолютная предельная погрешность может быть определена как

    $$D(u^*) = \sup\limits_{\mid {t - 1}\mid \le 0,1} \left|{t^{10} - 1}\right| = {(1.1)}^{10} - 1 \approx 1,5 \ldots$$

    Оценка погрешности u при вычислении значения функции по максимуму производной и линейная оценка соответственно будут $$D_{1}(u^{*}) = b\Delta (t^{*}) = 2,3 ...$$; $$D_{2}(u^{*}) = |\gamma (0)|\Delta (t^{*}) = 1.$$

  • Дать линейную оценку погрешности при вычислении неявной функции $$\varphi (u,t_1,t_2, \ldots ,t_n) = 0$$, если известны точка приближения $$\{ t_1^*, \ldots ,t_n^* \}$$, значение функции в точке приближения u* и погрешность в определении аргументов $$\Delta t_1^*, \ldots ,\Delta t_{n}^* .$$

    Решение. Дифференцируя по tj, получим

    $$$ \frac{\partial \varphi }{\partial u} \frac{\partial u} {\partial t_j} + \frac{\partial \varphi }{\partial t_j} = 0, $$$

    откуда

    $$$ \frac{\partial u}{\partial t_j} = - (\frac{\partial \varphi }{\partial t_j})(\frac{\partial \varphi }{\partial u})^{ - 1}. $$$

    При заданных $$\{ t_1^*, \ldots ,t_n^* \}$$, можно найти u* как корень уравнения $$\varphi (u,t_1 ,t_2 , \ldots ,t_{n} ) = 0 $$, а затем — значения

    $$$ b_j(0) = \left. { - (\frac{\partial \varphi }{\partial t_j })(\frac{\partial \varphi }{\partial u})^{ - 1} } \right|_{(u^*,t_1^* , \ldots ,t_n^* )} , $$$

    откуда можно получить линейную оценку погрешности функции D2(u*).

  • Вычислить относительную погрешность в определении значения функции$$u = xy^{2}z^{3}, \ если \ x^{*} = 37,1, y^{*} = 9,87, z^{*} = 6,052, \Delta x^{*} = 0,3, \Delta y^{*} = 0,11, \\ \Delta z^{*} = 0,016.$$

    Решение:

    $$$ \delta_x = \frac{0.3}{37.1} \approx 0.81 \cdot 10^{ - 2}, \delta_y = \frac{0,11}{9,87} \approx 1,12 \cdot 10^{ - 2}, \delta_z = \frac{0,016}{6,052} \approx 0,26 \cdot 10^{ - 2}, \\ \delta (u) = \delta (x^*) + 2\delta (y^*) + 3\delta (z^*) = 3.8 \cdot 10^{ - 2} . $$$
  • Оценить погрешность в определении корней квадратного уравнения $$\varphi (u,t_1,t_2) = u^2 + t_1 u + t_2 = 0$$, если заданы приближения $$t_1^*,t_2^*,\Delta (t_1^*),\Delta (t_2^*).$$

    Пусть u* — решение уравнения

    $$u*^2 + t_1^*u^* + t_2^* = 0.$$

    Из формулы

    $$$ b_j(0) = \left. { - (\frac{d\varphi }{dt_j})(\frac{d\varphi }{du})^{ - 1}}\right|_{(u^*,t_1^* , \ldots ,t_{n}^*)}$$$

    получим

    $$$ b_1 (0) = \left. {\frac{du}{dt_1 }} \right|_{(t_1^*,t_2^* )} = - \frac{u^*}{2u^* + t_1^*}, \\ b_2 (0) = \left. {\frac{du}{dt_2 }} \right|_{(t_1^*,t_2^* )} = - \frac{1}{2u^* + t_1^* }.$$$

    Следовательно, линейная оценка будет

    $$$ D_2 (u^*) = \frac{\left|{u^*}\right| \cdot \Delta (t_1^*) + \Delta (t_2^* )}{\left|{2u^* + t_1^*}\right|}.$$$
  • 1.8. Задачи для самостоятельного решения

  • Найти абсолютную предельную погрешность, погрешность по производной и линейную оценку погрешности для функций u = sin t,$$u = \frac{1}{t^2 - 5t + 6}.$$

    Заданы точка приближения t = t* и погрешность $$\Delta t.$$

  • Определить шаг $$\tau$$, при котором погрешность вычисления производной u'(t), приближенно вычисляемой в соответствии с формулами

    $$u^{\prime}(t) \approx \frac{f(t + \tau ) - f(t)}{\tau }, \\ u^{\prime}(t) \approx \frac{f(t + \tau ) - f(t - \tau )}{2\tau } ,$$

    не превосходит 103. Известно, что $$\left|{u^{\prime\prime}(t)}\right| \le 1,\quad \left|{u^{\prime\prime\prime}(t)}\right| \le 1$$ для любых t.

  • Пусть для вычисления функции u = f(t) используется частичная сумма ряда Маклорена

    $$$ u(t) \approx u(0) + \frac{u^{\prime}(0)}{1!}t + \ldots + \frac{u^{(n)}(0)}{n!}t^n $$$,

    причем аргумент задан с погрешностью $$\Delta t = 10^{ -3}.$$

    Найти n такое, чтобы погрешность в определении функции u(y) по данной формуле не превышала $$\Delta t.$$ Рассмотреть отрезки $$t \in [0,1],\quad t \in [10,11] .$$

    Предложить более совершенный алгоритм для вычисления функций u(t) = sin t, u(t) = et на отрезке $$t \in [10,11] .$$

  • Определить оптимальный шаг численного дифференцирования $$\tau$$ при использовании для вычисления производной приближенной формулы

    $$u^{\prime}(t) \approx \frac{u(t - 2\tau ) - 8(t - \tau ) + 8(t + \tau ) - u(t + 2\tau )}{12t},$$

    имеющей четвертый порядок точности, если известно, что $$\left|{u^{(5)}(t)}\right| \le M_5$$, а значения функций вычисляются с точностью $$\varepsilon.$$

  • Вычислить относительную погрешность в определении значения функции u(x,y,z) = x2y2/z4, если заданы$$x^{*} = 37,1, y^{*} = 9,87, z^{*} = 6,052, \Delta (x^{*}) = 0,1; \\ \Delta (y^{*}) = 0,05; \Delta (z^{*}) = 0,02.$$
  • Вернуться к учебному плану