вычисл матем


МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ РОССИЙСКОЙ ФЕДЕРАЦИИ
федеральное государственное бюджетное образовательное учреждение
высшего профессионального образования
«Тольяттинский государственный университет»
УЧЕБНО-МЕТОДИЧЕСКОЕ ПОСОБИЕ
по проведению практических занятий
по дисциплине
«ВЫЧИСЛИТЕЛЬНАЯ МАТЕМАТИКА»
по направлению
010300.68 Математическое обеспечение и администрирование информационных систем
Тольятти 2015

Содержание
TOC \o "1-3" \h \z \u 1Планы практических занятий PAGEREF _Toc434310026 \h 31.1Занятие № 1. Теория погрешности. PAGEREF _Toc434310027 \h 31.2Занятие № 2. Решения алгебраических уравнений. Отделение корней. PAGEREF _Toc434310028 \h 111.3Занятие № 3. Метод хорд (секущих) и метод деления пополам. PAGEREF _Toc434310029 \h 131.4Занятие № 4. Метод Ньютона (касательных). PAGEREF _Toc434310030 \h 171.5Занятие № 5. Комбинированный метод. PAGEREF _Toc434310031 \h 201.6Занятие № 6. Метод итераций. PAGEREF _Toc434310032 \h 211.7Занятие № 7. Прямые методы решения систем линейных алгебраических уравнений. PAGEREF _Toc434310033 \h 241.8Занятие № 8. Метод итераций для систем уравнений. PAGEREF _Toc434310034 \h 291.9Занятие № 9. Равномерное приближение функций многочленами. Интерполяционный многочлен Лагранжа. PAGEREF _Toc434310035 \h 351.10Занятие № 10. Интерполяционный многочлен Ньютона. Многочлены Чебышева. PAGEREF _Toc434310036 \h 441.11Занятие № 11. Интерполяция с кратными узлами и формула Эрмита. PAGEREF _Toc434310037 \h 651.12Занятие № 12. Понятие сплайна. Виды сплайнов. PAGEREF _Toc434310038 \h 691.13Занятие № 13. Численное дифференцирование. PAGEREF _Toc434310039 \h 761.14Занятие № 14. Численное интегрирование. Квадратурные формулы Ньютона-Котеса. PAGEREF _Toc434310040 \h 791.15Занятие № 15. Квадратурные формулы Гаусса. PAGEREF _Toc434310041 \h 881.16Занятие № 16. Численные методы решения задачи Коши для обыкновенных дифференциальных уравнений. Классификация методов. Метод Пикара. PAGEREF _Toc434310042 \h 931.17Занятие № 17. Методы Эйлера. PAGEREF _Toc434310043 \h 961.18Занятие № 18. Методы Рунге-Кутта. PAGEREF _Toc434310044 \h 1051.19Занятие № 19. Пошаговый контроль точности. Метод Кутты-Мерсона. PAGEREF _Toc434310045 \h 1101.20Занятие № 20. Многошаговые методы Адамса. PAGEREF _Toc434310046 \h 1141.21Занятие № 21. Методы прогноза и коррекции. PAGEREF _Toc434310047 \h 1191.22Занятие № 22. Метод Милна четвертого порядка PAGEREF _Toc434310048 \h 122Приложение 1 PAGEREF _Toc434310049 \h 126Приложение 2. PAGEREF _Toc434310050 \h 128Индивидуальное домашнее задание № 1 PAGEREF _Toc434310051 \h 129Индивидуальное задание №2. PAGEREF _Toc434310052 \h 134Индивидуальное задание №3 PAGEREF _Toc434310053 \h 139Индивидуальное домашнее задание № 4 PAGEREF _Toc434310054 \h 143Литература PAGEREF _Toc434310055 \h 147
Планы практических занятийЗанятие № 1. Теория погрешности.Источники и виды погрешности
На некоторых этапах решения задачи на ЭВМ могут возникать погрешности, искажающие результаты вычислений. Оценка степени достоверности получаемых результатов является важнейшим вопросом при организации вычислительных работ. Это особенно важно при отсутствии опытных или других данных для проверки адекватности модели.
Погрешность решения задачи обуславливается следующими причинами:
Математическое описание задачи является неточным, в частности, неточно заданы исходные данные.
Применяемый для решения метод часто не является точным: получение точного решения возникающей математической задачи требует неограниченного или неприемлемо большого числа арифметических операций, и, поэтому вместо получения точного решения задачи приходится прибегать к приближенному.
При вводе данных в машину, при выполнении арифметических операций и при выводе данных производятся округления. Погрешности, соответствующие этим причинам, называют:
неустранимой погрешностью,
погрешностью метода,
вычислительной погрешностью.
Часто неустранимую погрешность разделяют на две части:
а) неустранимой погрешностью называют лишь погрешность, являющуюся следствием неточности задания числовых данных, входящих в математическое описание задачи;
в) погрешность, являющуюся следствием несоответствия математического описания задачи реальности, называют погрешностью математической модели.
Математическая модель, принятая для описания данного процесса или явления, может внести существенные погрешности, если в ней не учтены какие-либо существенные черты рассматриваемой задачи. В частности, математическая модель может прекрасно работать в одних условиях и быть совершенно неприемлемой в других; поэтому важно правильно учитывать область её применимости.
Численный метод также является источником погрешностей. Это связано, например, с заменой интеграла суммой, усечением рядов при вычислениях значений функций, интерполированием табличных данных и т. п. Как правило, погрешность численного метода регулируема, т. е. она может быть уменьшена до любого разумного значения путем изменения некоторого параметра (например, шага интегрирования, числа членов усеченного ряда и т.п.). Погрешность метода обычно стараются довести до величины, в несколько раз меньшей погрешности исходных данных. Дальнейшее снижение погрешности не приведет к повышению точности результатов, а лишь увеличит стоимость расчетов из-за необоснованного увеличения объема вычислении.
При вычислениях с помощью ЭВМ неизбежны погрешности округлений, связанные с ограниченностью разрядной сетки машины. Обычно после выполнения операции производится не округление результата, а простое отбрасывание лишних разрядов с целью экономии машинного времени. Правда, в современных машинах предусмотрена свобода выбора программистом способа округления; соответствующими средствами располагают и некоторые алгоритмические языки.
Максимальная относительная погрешность при округлении есть δmax=0.5α1-k, где α - основание системы счисления, к—количество разрядов мантиссы числа. При простом отбрасывании лишних разрядов эта погрешность увеличивается вдвое.
В современных машинах с памятью, измеряемой в байтах, принята шестнадцатеричная система счисления, и любое число с плавающей точкой содержит шесть значащих цифр. Следовательно, α =16, к = 6, максимальная погрешность округления αmах= 0.5* 16-5 ≈ 0.5* 10-8.
Несмотря на то, что при решении больших задач выполняются миллиарды операций, это вовсе не означает механического умножения погрешности при одном округлении на число операций, так как при отдельных действиях погрешности могут компенсировать друг друга (например, при сложении чисел разных знаков). Вместе с тем иногда погрешности округлений в сочетании с плохо организованным алгоритмом могут сильно исказить результаты.
Перевод чисел из одной системы счисления в другую также может быть источником погрешности из-за того, что основание одной системы счисления не является степенью основания другой (например, 10 и 2). Это может привести к тому, что в новой системе счисления число становится иррациональным.
Например, число 0.1 при переводе в двоичную систему счисления примет вид 0.1=0.00011001100... Может оказаться, что с шагом 0.1 нужно при вычислениях пройти отрезок [0,1] от х=1 до х=0; десять шагов не дадут точного значения х=0.
Абсолютная и относительная погрешности
Если А - точное значение некоторой величины, a a - известное приближение к нему, то абсолютной погрешностью приближения a числа А называют некоторую величину Δ(a) удовлетворяющую условию:

Относительной погрешностью называют некоторую величину δ(А), для которой выполняется условие:

Относительную погрешность часто выражают в процентах.
Информацию о том, что а является приближенным значением числа А с абсолютной погрешностью Δ(а), принято записывать в виде:

Числа а и Δ(а) записываются с одинаковым количеством знаков после запятой.
Информацию о том, что а является приближенным значением числа А с относительной погрешностью δ(а), записывают в виде:

Число δα, заведомо не меньшее относительной погрешности называют предельной относительной погрешностью, т.е.

Т.к. на практике A≈ α, то приближенно можно принять, что


Погрешность функции
Пусть искомая величина Y является функцией параметров a1, а2, …, an; т.е. Y = Y(a), и известна область G в пространстве переменных a1, а2, …, an, которой принадлежат параметры. Необходимо получить приближение y к Y и оценить его погрешность.
Если у - приближённое значение величины Y, то предельной абсолютной погрешностью А(у) называют наилучшую при имеющейся информации оценку погрешности величины у:

Предельной относительной погрешностью δ(у) называют величину

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

Для оценки предельной относительной погрешности функции имеют место выражения

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

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

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

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

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

Решение: Произведем вычисления

Определим относительные погрешности аргументов




Устойчивость, корректность, сходимость
Устойчивость. Рассмотрим погрешности исходных данных. Поскольку это так называемые неустранимые погрешности и вычислитель но может с ними бороться, то нужно хотя бы иметь представление об их влиянии па точность окончательных результатов. Конечно, мы вправе надеяться на то, что погрешность результатов имеет порядок погрешности исходных данных. Всегда ли это так? К сожалению, нет. Некоторые задачи весьма чувствительны к неточностям в исходных данных. Эта чувствительность характеризуется так называемой устойчивостью.
Пусть в результате решения задачи по исходному значению величины х находится значение искомой величины у. Если исходная величина имеет абсолютную погрешность Δх, то решение имеет погрешность Δу. Если решение непрерывно зависит от входных данных, т.е. всегда ||Δу||→0 при ||Δх||→0 , то задача называется устойчивой по входным данным.
Отсутствие устойчивости обычно означает, что даже незначительные погрешности в исходных данных приводят к большим погрешностям в решении или вовсе к неверному результату. О подобных устойчивых задачах также говорят, что они чувствительны к погрешностям исходных данных.
Примером такой задачи является отыскание действительных корней уравнения вида: tg(x) -х = 1.Изменение аргумента х на малую величину Δх может привести к существенному изменению функции.
Корректность. Задача называется поставленной корректно, если для любых значений исходных данных из некоторого класса ее решение существует, единственно и устойчиво по исходным данным.
Рассмотренные выше примеры неустойчивой задачи являются некорректно поставленными. Применять для решения таких задач численные методы, как правило, нецелесообразно, поскольку возникающие в расчетах погрешности округлений будут сильно возрастать в ходе вычислений, что приведет к значительному искажению результатов.
Вместе с тем отметим, что в настоящее время развиты методы решения некоторых некорректных задач. Это в основном так называемые методы регуляризации. Они основываются на замене исходной задачи корректно поставленной задачей. Последняя содержит некоторый параметр, при стремлении которого к нулю решение этой задачи переходит в решение исходной задачи.
Неустойчивость методов. Иногда, даже если задача устойчива, то численный алгоритм может быть неустойчивым. Например, если производные заменяются разностями, то приходится вычитать близкие числа и сильно теряется точность. Эти неточные промежуточные результаты используются в дальнейших вычислениях, и ошибки могут сильно нарастать.
Понятие сходимости. При анализе точности вычислительного процесса одним из важнейших критериев является сходимость численного метода. Она означает близость получаемого численного решения задачи к истинному решению. Строгие определения разных оценок близости могут быть даны лишь с привлечением аппарата функционального анализа. Здесь мы ограничимся некоторыми понятиями сходимости, необходимыми для понимания последующего материала.
Рассмотрим понятие сходимости итерационного процесса. Этот процесс состоит в том, что для решения некоторой задачи и нахождения искомого значения определяемого параметра (например, корня нелинейного уравнения) строится метод последовательных приближений. В результате многократного повторения этого процесса (или итераций) получаем последовательность значении х1, х2, …, хп ,... эта последовательность сходится к точному решению х=а, если при неограниченном возрастании числа итераций предел этой последовательности существует и равен a: limxn = а. В этом случае имеем сходящийся численный метод.
Другой подход к понятию сходимости используется в методах дискретизации. Эти методы заключаются в замене задачи с непрерывными параметрами на задачу, в которой значения функций вычисляются в фиксированных точках. Это относится, в частности, к численному интегрированию, решению дифференциальных уравнений и т. п. Здесь под сходимостью метода понимается стремление значений решения дискретной модели задачи к соответствующим значениям решения исходной задачи при стремлении к нулю параметра дискретизации (например, шага интегрирования).
Таким образом, для получения решения задачи с необходимой точностью, ее постановка должна быть корректной, а используемый численный метод должен обладать устойчивостью и сходимостью.
Задание. Выполнить задание 1 ИДЗ№1.
Занятие № 2. Решения алгебраических уравнений. Отделение корней.Уравнение f(x)=0 называется алгебраическим уравнением n-й степени, если f(x) является полиномом,
Всякое неалгебраическое уравнение называется трансцендентным. Любое значение A, обращающее функцию f(x) в нуль называется корнем уравнения f(x)=0.
Для алгебраических уравнений, степень которых выше четырех, не существует формул, которые выражали бы величины корней через коэффициенты уравнений. Сравнительно редко удается найти точное значение корней и трансцендентных уравнений. Поэтому, важное значение приобретают методы приближенного нахождения корней уравнения и оценка степени их точности.
Сам процесс вычисления корней состоит из двух операций:
отделение корней, т.е. установление возможно тесных промежутков [a,b], в которых содержится один и только один корень уравнения;
уточнение приближенных корней, т.е. доведение их до заданной степени точности.
Для отделения корней часто используют теорему Больцано-Коши.
Теорема: Если непрерывная функция f(x) принимает значение разных знаков на концах отрезка [a,b] то внутри этого отрезка содержится по крайней мере один корень уравнения f(x)=0. Корень заведомо будет единственным, если производная f'(x) существует и сохраняет знак внутри интервала [a,b].
Отделение корней уравнения f(x)=0 можно выполнить графически, построив график функции f(x), по которому можно судить о том, в каких промежутках находятся точки пересечения его с осью 0х. В некоторых случаях целесообразно представить исходное уравнение в эквивалентном виде: с таким расчетом, чтобы графики функций y=f1(x) и y=f2(x) строились проще, чем график y=f(x). Корень уравнения f(x)=0 представляет собой абсциссу точки пересечения графиков y=f1(x) и y=f2(x).
Пример. Отделить корни уравнения
Решение. 1-способ. В данном случае f(x)=x 3+2 х -1, . Поскольку f'(x >0 при всех x, то функция f(x) возрастает в промежутке Корень считается отделенным, если указан конечный промежуток [a,b], на котором он находится.Методом проб находим отрезок [a,b], для которого на концах отрезка функция f(x) принимает значения разных знаков. Для этого вычислим значения функции при некоторых значениях аргумента
f(-1)=-4<0, f(0)=-1<0, f(1)=2>0.
Согласно теореме Больцано-Коши корень находится на отрезке [0,1], причем он единственный, т.к. производная f'(x)>0 сохраняет знак внутри интервала [0,1].

2-способ. Корень данного уравнения можно отделить и графически. Придадим уравнению вид х3=-2х+1, т. е.
вид f1(x)=f2(x), и построим графикифункций у=х3и y=-2 х +1.
Эти графики пересекаются в точке М, абсцисса которой принадлежит интервалу (0,1).
Задания. Выполнить задание 2.1 ИДЗ№1.
Занятие № 3. Метод хорд (секущих) и метод деления пополам.Метод деления пополам (метод дихотомии).
Наиболее надёжным алгоритмом нахождения корня уравнения f(x)=0, особенно когда о поведении функции f(x) мало что известно, является метод половинного деления.
Пусть f(x) - функция действительной переменной и пусть известен интервал [a,b], на котором функция меняет знак, следовательно, между а и b существует точка, в которой функция обращается в нуль. Если разделить интервал пополам и определить положительна или отрицательна функция в точке деления, то тем самым найдём подынтервал, в котором функция меняет знак.
В принципе, повторным применением этого приема (деление интервала пополам) можно сколь угодно близко "подойти к корню". Так как каждый шаг делит интервал, в котором лежит корень, пополам, то 10 шагов уменьшают интервал в 210=1024 раз. При заданной абсолютной точности ε алгоритм метода деление пополам состоит из следующих шагов:
Вычислить f(a) и f(b).
Положить с = (а + b) / 2. Вычислить f(c).
Если sign(f(c))=sign(f(a)), то заменить а на с; в противном случае заменить b на с (функция sign(f(c)) означает знак в точке с).
Если b - a >ε, то перейти к шагу 2; в противном случае прекратить вычисления, поскольку мы достигли требуемой точности. Любой из концов отрезка или их полусумма может быть использован в качестве корня уравнения f(x) = 0.
Алгоритм деления пополам довольно медлителен, но зато абсолютно застрахован от неудач.
Пример. Пользуясь методом половинного деления вычислить с точностью до ε=0,1 корень уравнения x3+3x2-3=0, заключенный в отрезке [-3,-2].
Решение. Описанный выше процесс решения удобно оформить в виде
вычислительного бланка:
n an bn f(an) f(bn) cn f(cn) (bn-an)/2
0 -3 -2 -3 1 -2,5 0,125 0.5
1 -3 -2,5 -3 0,125 -2,75 -1,11 0.25
2 -2,75 -2,5 -1,11 0,125 -2,625 -0,42 0.125
3 -2,625 -2,5 -0,42 0,125 -2,5625 -0,129 0.0625
Из таблицы следует, что A=-2,5625±0,0625. Или округлив результат, получаем А=-2,6±0,1.
Метод пропорциональных частей (метод хорд).
Мы будем предполагать выполнение следующих условий:
функция f(x) в промежутке [a, b] непрерывна вместе со своими производными
значения f(a) и f(b) функции на концах промежутка имеют разные знаки: т.е. f(a)*f(b)<0;
обе производные f' (x) и f'' (x) сохраняют каждая определенный знак во всем промежутке [а, b].
Тогда возможны следующие четыре случая, изображенные на рисунке 1.
204470-702310

1. Рассмотрим случай f(a)f"(x)<0 , которому соответствуют графики, изображенные на рис.1б или 1в. Заменим дугу ММ' кривой (рис.1а)- хордой ММ'. Уравнение хорды может быть написано, например, в виде
(1)
Наше правило, по существу, сводится к тому, что вместо точки А - пересечения кривой f(x) с осью 0x, определяется точка D - пересечение с осью 0x соответственно хорды MM'. Полагая у=0 в формуле (1), получим значение точки a1:
(2)
Через точки (a1, f(a1)) и (b, f(b)) новую хорду получим точку a2. Продолжая этот процесс, получим последовательность точек сходящихся к корню А. Формулы для расчета последовательности an следующие
(3)
2. В случае f(a)f"(x)>0 (см. рис.1а или 1г) приближение к корню А происходит со стороны точки b, а последовательность точек определяется по рекуррентной формуле
(4)
Для оценки погрешности метода хорд будем предполагать, что производная функции f(x) непрерывна и сохраняет знак на отрезке [a,b]. Тогда, если А -точный корень, а xn - приближение к корню, то по формуле конечных приращений (теорема Лагранжа), имеем
(5)
где точка с лежит между А и xn.
Из формулы (2) получаем
(6)
где m наименьшее значение модуля производной функции f(x) на отрезке [a,b], т.е.
Другая формула для оценки погрешности метода хорд
(7)
где М и m соответственно наибольшее и наименьшее значение модуля производной функции f(x).
Пример. Найти корень уравнения x3 -2x2 -4x-7=0 интервале [3,4], с точность до 0,01.
Решение. Обозначим левую часть уравнения через f(x). Подсчитаем
f(3)=-10<0и< f(4) = 9>0. Нетрудно показать, что обе производные f'(x)=3x2- 4x-4, f''(x ) = 6x- 4 в промежутке[3,4] положительны.
Действительно, для f"(x) это очевидно, а т.к. f"(x)>0, то f'(x) монотонно возрастает. Так как f'(3) = 11>0, то в других точках отрезка [3,4] значения производной будут заведомо не меньше 11. Попутно мы показали, что наименьшее значение первой производной равно m= f'(3) = 11.
Т.к. f(3)f"(x)<0, то для дальнейших расчетов следует воспользоваться формулами (3) и (6). Обозначим
, тогда получим, что

Описанный выше процесс решения удобно оформить в виде вычислительного бланка.
Процесс прервался после третьего шага, т. к. мы достигли нужную точность: εn=0,004<0,01. Приближенное значение корня А=3,630±0,01.
Задания. Выполнить задание 2.2 ИДЗ№1.
Занятие № 4. Метод Ньютона (касательных).Цель - ознакомить студентов с методом Ньютона решения алгебраических уравнений.
Метод касательных (Метод Ньютона)
Пусть корень уравнения f(x)=0 отделен на отрезке [a,b], причем f'(x) и f''(x) непрерывны и сохраняют свои знаки на этом отрезке. Тогда, как было сказано выше, возможны четыре случая, графически изображенные на рис.1 а)-г).
Геометрически метод Ньютона эквивалентен замене дуги кривой ММ' касательной, проведенной к некоторой точке данной кривой.

1. Рассмотрим случай f(a)f"(x)<0 , которому соответствуют графики изображенные на рис.1 б) или 1 в). В этом случае, приближение к корню А происходит со стороны точки b (см. рис.2). Касательные проводятся в точках M0, M1, M2, ... .Пересечением касательных с осью 0х является последовательность точек b0,b1,b2, которая сходится к точному корню А. Уравнение касательной к функции f(x) в точке Mn с координатами (bn ,f(bn)) имеет вид:
(8)
Полагая y=0 , x=bn+1 получаем
(9)
Заметим, что если бы мы в данном случае провели касательную в точке (a,f(a)), то точка пересечения этой касательной с осью 0х лежала бы вне отрезка [a,b], и требуемой сходимости к корню не наблюдалось бы.
2. В случае f(a)f"(x)>0 (см. рис.1 а или 1г) приближение к корню А происходит со стороны точки a, а последовательность точек определяется по рекуррентной формуле
(10)
Для оценки погрешности метода касательных будем предполагать, что обе производные функции f(x) непрерывны и сохраняют свой знак на отрезке [a,b]. Тогда, если А - точный корень, а xn - приближение к корню, то оценка погрешности может быть вычислена по одной из следующих формул
(11)
(12)
где
Если , то формула (12) упрощается
(12*)
Пример. Вычислить приближенно по методу касательных с точностью до пяти десятичных знаков после запятой больший отрицательный корень уравнения x3 -12x-8 = 0.
Решение. Графически отделяя корни данного уравнения, заключаем, что уравнение имеет три действительных корня; больший отрицательный корень принадлежит отрезку [-1,0]. Укажем отрезок меньшей длины, на котором находится корень; это отрезок [-0,7;-0,65], поскольку f(-0,7) = 0,057, f(-0,65)= -0,474625.
Т.к. f"(x) = 6x, то очевидно для отрицательных x имеем f"(x)<0, то выполнено условие f(a)f"(x)<0, поэтому в качестве нулевого приближения берем b0=-0,65, а последующие приближения вычисляем по формуле (9).
Описанный выше процесс решения удобно оформить в виде вычислительного бланка.

Процесс прервался, после третьего шага, т.к. мы достигли нужную точность, которая в данном случае может быть рассчитана по формуле (12*). Приближенное значение корня А=-0,694593±0,000001.
Задания. Выполнить задание 2.2 ИДЗ№1.
Занятие № 5. Комбинированный метод.Цель - ознакомить студентов с комбинированным методом решения алгебраических уравнений.
Данный метод основан на построении схематического графика функции, определении интервалов его пересечения с осью абсцисс и последующим «сжатием» этого интервала при помощи строимых хорд и касательных к графику этой функции.
Преимущество комбинированного метода заключается в «двустороннем сжатии» рассматриваемого отрезка, что наглядно показано на рис.3. Из рис.3 также видно, что истинный корень А лежит между приближениями со стороны точки a и со стороны точки b, т.е.а n < A < bn .

Также как, и отдельно в методе хорд и в методе касательных, для расчета приближений мы рассмотрим два случая.
1. Случай
(13)
2. Случай
(14)
Если задана допустимая погрешность ε, то процесс уточнения корня прекратится, как только За значение корня А следует взять
Задания. Выполнить задание 2.2 ИДЗ№1.
Занятие № 6. Метод итераций.Цель - ознакомить студентов с методом итераций решения алгебраических уравнений.
Если каким-либо способом получено приближенное значение x0 корня уравнения f(x)=0, то уточнение корня можно осуществить методом последовательных приближений или методом итераций. Для этого уравнение
f(x)=0 представляют в виде
x = φ (x) (15)
что всегда можно сделать, и притом многими способами, например,
x = x + c • f(x) (16)
где с - произвольная постоянная.
Пусть число x1 есть результат подстановки x0 в правую часть уравнения (15), т.е. x1=φ(x0). Далее, x2=φ(x1), x3=φ(x2),...,
xn+1 =φ (xn) (17)
Процесс последовательного вычисления чисел xn (n=1,2,3,...). по формуле (17) называется методом последовательных приближений или методом итераций.
Утверждение. Итерационный процесс сходится (здесь А- точный корень исходного уравнения f(x)=0), если на отрезке [a,b] содержащем корень А и его последовательные приближения xn, выполнено условие
(18)
Замечание. В качестве первого приближения x0 можно взять произвольное значение из интервала, содержащего корень A, причем такой интервал можно сделать достаточно малым.
Для оценки погрешности метода итераций используется следующая формула
(19)
При оценка погрешности упрощается
(20)
Пример 1 . Методом итераций найти меньший положительный корень уравнения x3 -5x+1=0.
Решение. Графически отделяя корни данного уравнения, заключаем, что уравнение имеет три действительных корня, лежащих в отрезках [-3;-2], [0;1], [2;3]. Найдем корень принадлежащий отрезку [0; I]. Укажем отрезок меньшей длины, на котором находится корень. Поскольку, f(x)=x3 -5x+1, f(0)=1>0,, то корень принадлежит отрезку [0; 0,5]. Данное уравнение приведем к виду (15), разрешая его относительно x:
т.е.
Так как, то условие (18) выполняется, следовательно процесс итераций будет сходиться, а для оценки погрешности можно использовать (20).
Взяв в качестве начального приближения середину отрезка, т.е. приняв x0=0,25, вычисление последующих приближений проведем по формуле

Результаты этих вычислений представлены в таблице_1, из которой видно, что искомый корень A= 0,20164, вычислен с заданной точностью ε=0,0001, т.к. согласно формуле (20) имеем


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

т.е. условие (18) не выполняется. В этом случае данное уравнение следует представить в другом виде, например, для функции условие (18) на отрезках [-3;-2], [2;3] будет выполняться.
Пример 2. Методом итераций найти отрицательный корень уравнения
x4 + x -3 = 0
Решение. Данное уравнение имеет два действительных корня; отрицательный корень находится на отрезке [-1,5;-1,4], так как для его концов выполняется условие

Уравнение запишем в виде (16):
x = x + c • ( x 4 + x — 3 )
где с - произвольная постоянная. Выберем значение с так, чтобы для функции на отрезке [-1,5;-1,4] выполнялось условие (18). В качестве такого значения можно взять с=0,1, тогда
φ (x) = x + 0,1 • (x4 + x — 3 ) ,

т.е. выполнено условие (18).
Возьмем x0=-1,45. Вычисление последующих приближений осуществим по формуле

и представим результаты в таблице_2, из которой видно, что A=-1,45262 является приближенным корнем данного уравнения.
Таблица_2

Задания: Выполнить задание 2.2 и 3 ИДЗ№1.
Занятие № 7. Прямые методы решения систем линейных алгебраических уравнений.Существует несколько методов, позволяющих получить решение системы n линейных уравнений с n неизвестными приближенно, с заданной точностью.
Способы решения систем линейных уравнений делятся на две группы:
точные методы, представляющие собой конечные алгоритмы для вычисления корней системы (решение систем с помощью обратной матрицы, правило Крамера, метод Гаусса и др.),
итерационные методы, позволяющие получить решение системы с заданной точностью путем сходящихся итерационных процессов (метод итерации, метод Зейделя и др.).
Вследствие неизбежных округлений результаты даже точных методов являются приближенными. При использовании итерационных методов, сверх того, добавляется погрешность метода.
Эффективное применение итерационных методов существенно зависит от удачного выбора начального приближения и быстроты сходимости процесса.
Метод Гаусса
1. Схема единственного деления. Рассмотрим сначала простейший вариант метода Гаусса, называемый схемой единственного деления.
Прямой ход состоит из n - 1 шагов исключения.
1-й шаг. Целью этого шага является исключение неизвестного x1 из уравнений с номерами i = 2, 3, n. Предположим, что коэффициент a11 ≠ 0. Будем называть его главным элементом 1-го шага.
Найдем величины

называемые множителями 1-го шага. Вычтем последовательно из второго,третьего, …, n-го уравнений системы первое уравнение, умноженное соответственно на q21, q31, …, qn1. Это позволит обратить в нуль коэффициенты при x1 во всех уравнениях, кроме первого. В результате получим эквивалентную систему


в которой вычисляются по формулам

2-й шаг. Целью этого шага является исключение неизвестного х2 из уравнений с номерами i = 3, 4, n. Пусть a22(1) ≠ 0, где a22(1) - коэффициент, называемый главным (или ведущим) элементом 2-го шага. Вычислим множители 2-го шага

и вычтем последовательно из третьего, четвертого, ..., n-го уравнения системы второе уравнение, умноженное соответственно на q32, q42, …, qm2. В результате получим систему

Здесь коэффициенты вычисляются по формулам

Аналогично проводятся остальные шаги. Опишем очередной k-й шаг.
к-й шаг. В предположении, что главный (ведущий) элемент k-го шага akk(k-1) отличен от нуля, вычислим множители k-го шага

и вычтем последовательно из (k + 1)-го, …, n-го уравнений полученной на предыдущем шаге системы k-e уравнение, умноженное соответственно на
После (n - 1)-го шага исключения получим систему уравнений

матрица которой является верхней треугольной. На этом вычисления прямого хода заканчиваются.
Обратный ход. Из последнего уравнения системы находим xn. Подставляя найденное значение xn в предпоследнее уравнение, получим xn-1. Осуществляя обратную подстановку, далее последовательно находим xn-1, xn-2, ., x1. Вычисления неизвестных здесь проводятся по формулам

Необходимость выбора главных элементов. Заметим, что вычисление множителей, а также обратная подстановка требуют деления на главные элементы akk(k-1). Поэтому, если один из главных элементов оказывается равным нулю, то схема единственного деления не может быть реализована. Здравый смысл подсказывает, что и в ситуации, когда все главные элементы отличны от нуля, но среди них есть близкие к нулю, возможен неконтролируемый рост погрешности.
Метод Гаусса с выбором главного элемента по всей матрице (схема полного выбора).
В этой схеме допускается нарушение естественного порядка исключения неизвестных.
На 1-м шаге метода среди элементов aiJ определяют максимальный по модулю элемент Первое уравнение системы и уравнение с номером i1 меняют местами. Далее стандартным образом производят исключение неизвестного из всех уравнений, кроме первого.
На k-м шаге метода среди коэффициентов при неизвестных в уравнениях системы с номерами i = k, ., n выбирают максимальный по модулю коэффициент Затем k-е уравнение и уравнение, содержащее найденный коэффициент, меняют местами и исключают неизвестное из уравнений с номерами i = k + 1, ., n.
На этапе обратного хода неизвестные вычисляют в следующем порядке:

Пример. Решить систему методом Гаусса с выбором главного элемента:

Решение системы удобно оформить в виде таблицы

Используя выделенные строки, получим треугольную систему для определения неизвестных

Откуда последовательно находим
Задания: выполнить задание 4 ИДЗ№1.
Занятие № 8. Метод итераций для систем уравнений.Цель - ознакомить студентов с методом итераций решения систем линейных алгебраических уравнений.
Рассмотрим систему линейных алгебраических уравнений (СЛАУ):
a11 x1 + a12 x2 + … + a1n xn = b1
a21 x1 + a22 x2 + … + a21 xn = b2
……………………………………. .
an1 x1 + an1 x2 +… + ann xn = bn,
или в векторно-матричном виде:
Ax = B, (1)
где
А =a11 a12 … a1na21 a22 … a2n… an1 an2 … ann , B=b1b2…bn, x=x1x2…xn
Итерационные методы решения СЛАУ используются для решения СЛАУ большой размерности с разреженными матрицами, а также для уточнения решения СЛАУ, полученного с помощью прямого метода. Формулировка и применение итерационных методов требует определенных знаний и определенного опыта. Выбор эффективного итерационного метода решения конкретной задачи существенно зависит от ее характерных свойств и от архитектуры вычислительной машины, на которой будет решаться задача. Поэтому никаких общих правил выбора наилучшего итерационного метода решения не существует. Метод простой итерации приведен здесь как иллюстрация действия механизма вычисления решения на основе итерационной процедуры.
Суть метода состоит в следующем. От системы уравнений вида (1) переходят к системе уравнений
x=Dх + С (2)
Например, от системы уравнений
а11х1+а12х2+а13х3=в1 а21х1+ а22х2+а23х3=в2 (3)
а31х1+а32х2+а33х3=в3
можно перейти к виду (2), выразив из первого уравнения х1, из второго - х2, из третьего - х3:
х1= - а12/а11х2 - а13/а11х3+в1/а11
х2= - а21/а22х1 - а23/а22х3+в2/а22 (4)
х3= - а31/а33х1 - а32/а33х2+ в3/а33
Приведение исходной системы уравнений в виду (2) можно осуществить различными способами. Например, в СЛАУ (3) из первого уравнения можно выразить х2, из второго - х1, из третьего - х3 и, переставив уравнения для сохранения порядка следования переменных в векторе решения х, снова прийти к виду (2). Естественно, что матрица D и вектор с будут уже иными. Возможны и другие способы преобразования исходных уравнений.
После преобразования (1) к виду (2) назначается нулевое приближение решения х (0):
х1 (0)
х (0) = х2 (0)
х3 (0).
Если приблизительно известны значения хi вектора решения х, то они выбираются в качестве нулевого приближения, если нет, то в качестве вектора х (0) выбирается любой вектор, например х (0) =О.
Первый шаг итерационного процесса состоит в вычислении приближения х(1):
х (1) = Dx (0) +С.
Например, назначив х (0) и подставив его в систему уравнений (4), получим:
х1 (1) = - а12/а11х (0) - а13/а11х3 (0) +в1/а11
х2 (1) = - а21/а22х1 (0) - а23/а22х3 (0) +в2/а22
х3 (1) = - а31/а33х1 (0) - а32/а33х2 (0) +в3/а33.
Далее вычисляем:
х (2) = Dx (1) +C
х (к) =Dх (к-1) +С
и т.д.
Достаточное условие сходимости метода итерации заключается в следующем, если норма матрицы D (обозначается ║D║) меньше 1, то система уравнений (2) имеет единственное решение х* и итерации сходятся к этому решению со скоростью геометрической прогрессии. Иными словами, если
║D║<1, (5)
то
limк→∞║х (к) - х*║ = 0
и выполняется тождество
х*=Dх*+С.
В качестве нормы матрицы D используются нормы ║D║1 или ║D║∞:
║D║1 = max1≤i≤nj=1n| dij |,
║D║∞= max1≤j≤ni=1n| dij ||.
Аналогично вводятся нормы вектора х:
║х║1 = i=1n| xi |║х║∞= max1≤i≤n|xi|.
Из условия сходимости (5) ясно, что не всякое преобразование исходной системы к виду (2) позволит получить решение уравнения на основе итерационного процесса, а только такое, которое обеспечит выполнение условия ║D║<1. Важно иметь в виду, что при выполнении этого условия итерационный процесс сходится для любого начального приближения х (0) и выбор х (0) =O диктуется просто соображениями удобства назначения х (0).
Условие (5) выполняется, если модули диагональных коэффициентов для каждого уравнения системы больше суммы модулей недиагональных коэффициентов, т. е.

Если задана допустимая погрешность вычислений Δ, то для оценки погрешности к - го приближения широко используется следующее неравенство:
║х (к) - х*║≤║D║ ∕ (1-║D║) •║х (к) - х (к-1) ║<Δ (6)
Из этого неравенства следует критерий окончания итерационного процесса
║х (к) - х (к-1) ║ < (1-║ D║) •∆ ∕ ║D║ (7)
Каждый раз при вычислении очередного приближения х (k) проверяется выполнение неравенства (7).
Выполнение неравенства (7) означает выполнение неравенства
║х (к) - х*║ < ∆
и, следовательно, прекращение итерационного процесса.
Метод Зейделя.
Метод Зейделя представляет собой некоторую модификацию метода итераций. Основная его идея заключается в том, что при вычислении (k + 1)-го приближения неизвестной xi учитываются уже вычисленные ранее (k + 1)-е приближения неизвестных x1, x2, …, xi - 1.
Пусть для системы (1) получена эквивалентная система

Выберем произвольно начальные приближения корней . Далее, предполагая, что k-ые приближения корней известны, согласно Зейделю будем строить (k + 1)-е приближения корней по формулам:

Метод прогонки.
Метод прогонки является модификацией метода Гаусса для частного случая разреженных систем - системы уравнений с трёхдиагональной матрицей. Такие системы получаются при моделировании некоторых инженерных задач, а также при численном решении краевых задач для дифференциальных уравнений.
Запишем систему уравнений в виде:

На главной диагонали матрицы этой системы стоят элементы над ней - элементыпод ней - элементы При этом обычно все коэффициенты bi не равны нулю.
Метод прогонки состоит из двух этапов - прямой прогонки (аналога прямого хода метода Гаусса) и обратной прогонки (аналога обратного хода метода Гаусса). Прямая прогонка состоит в том, что каждое неизвестное хi выражается через xi+1 с помощью прогоночных коэффициентов AAi, ВВi:
(1)
Из первого уравнения системы найдём

С другой стороны, по формуле (1) Приравнивая коэффициенты в обоих выражениях для х1, получаем
(2)
Из второго уравнения системы выразим х2 через х3, заменяя xi по формуле (1):

Отсюда найдём

или

Аналогично можно вычислить прогоночные коэффициенты для любогономера i:
(3)
Обратная прогонка состоит в последовательном вычислении неизвестных xi. Сначала нужно найти хn. Для этого воспользуемся выражением (1) при i=n-1 и последним уравнением системы. Запишем их:

Отсюда

Далее, используя формулы (1) и выражения для прогоночных коэффициентов (2), (3), последовательно вычисляем все неизвестные xi для i=n-1, n-2,.. .,1.
При анализе алгоритма метода прогонки надо учитывать возможность деления на нуль в формулах (3). Можно показать, что при выполнении условия преобладания диагональных элементов, т.е. если

причём хотя бы для одного значения i имеет место строгое неравенство, деления на нуль не возникает, и система имеет единственное решение.
Приведённое условие преобладания диагональных элементов обеспечивает также устойчивость метода прогонки относительно погрешностей округлений. Последнее обстоятельство позволяет использовать метод прогонки для решения больших систем уравнений. Заметим, что данное условие устойчивости прогонки является достаточным, но не необходимым. В ряде случаев для хорошо обусловленных систем метод прогонки оказывается устойчивым даже при нарушении условия преобладания диагональных элементов.
Задания: выполнить задание 5 ИДЗ№1.
Занятие № 9. Равномерное приближение функций многочленами. Интерполяционный многочлен Лагранжа.В основе большинства численных методов математического анализа лежит подмена одной функции f(x) (известной, неизвестной или частично известной) другой функцией φ(х) близкой к f(x) и обладающей «хорошими» свойствами, позволяющими легко производить над нею те или иные аналитические или вычислительные операции. Будем называть такую подмену аппроксимацией или просто приближением функции f(x) функцией φ(х). Для того, чтобы построить какую-то разумную теорию таких приближений и предложить конкретные способы получения аппроксимирующих функций φ (х) по заданным тем или иным образом аппроксимируемым функциям f(x), предварительно следует ответить на ряд вопросов.
Что известно о функции f(x)? Задана ли она своим аналитическим выражением или таблицей своих значений, какова степень ее гладкости и доступны ли значения ее производных, как расположены точки в интересующей части области определения f(x), где известны ее значения, и можно ли их задавать по своему усмотрению, и т.п.
Какому классу {семейству) функций должна принадлежать функция φ(х)? Какие дополнительные требования предъявляются к φ(х), выделяющие ее из заданного класса?
Что понимать под близостью между f(x) и φ(х); иначе, какой принять критерий согласия между ними? Говоря языком функционального анализа, по метрике какого пространства должно быть малым расстояние между f(x) и φ(х)?
Как видим, задача аппроксимации функции f(x) функцией φ(х) состоит в построении для заданной функции f(x) такой функции φ(х), что
(1)
причем левая часть приближенного равенства (1) должна быть обусловлена ответами на вопросы первой группы, правая часть — второй группы, а ответ на вопрос 3) должен уточнить значение связывающего f(x) и φ(х) символа «≈» .Прежде всего, определимся с ответом на второй вопрос. Договоримся использовать в качестве аппроксимирующих функций φ(х) только многочлены или функции, составленные из многочленов в таком случае будем говорить о полиномиальной аппроксимации или кусочно-полиномиальной аппроксимации соответственно.
По сравнению с другими семействами функций, пригодных для построения теории приближений, например, таких, как тригонометрические или показательные функции, рациональные функции, для вычислительной математики многочлены привлекательны тем, что они являются линейными функциями своих параметров (коэффициентов), и их вычисление сводится к выполнению конечного числа простейших арифметических операций — сложения и умножения.
Будем считать, что аппроксимация функции f(x) производится с помощью многочленов степени п∈N0. Тогда в зависимости от выбора критерия согласия и, в частности, от количества точек согласования f(x) с φ(х) (будем называть их узлами), т.е. точек, в которых известна информация об f(x) и, возможно, ее производных, можно рассмотреть разные конкретные способы аппроксимации.
Таблица 1

Интерполяционный многочлен Лагранжа
Пусть в точках таких, что известны значения функции у=f(x), т.е. на отрезке [а, b] задана табличная (сеточная) функция
(2)
Функция φ(х) называется интерполирующей или интерполяционной для f(x) на [а, b], если ее значения в заданных точках х0, х1, …, хп, называемых узлами интерполяции, совпадают с заданными значениями функции f(x), т.е. с у0, у1, …, уn соответственно. Геометрически факт интерполирования означает, что график функции φ(х) проходит так, что, по меньшей мере, в n+1 заданных точках он пересекает или касается графика функции f(х) (рис.1).

Рис.1. Геометрическая интерпретация задачи интерполирования
Легко представить, что таких графиков φ(х), проходящих через заданные точки, можно изобразить сколько угодно, и они могут отличаться от графика f(x) сколь угодно сильно, если не накладывать на φ(х) и f(x) определенных ограничений.
Будем считать, что интерполяционная функция φ(х) есть многочлен степени п. Тогда задача интерполяции, точнее, полиномиальной, алгебраической или параболической интерполяции (поскольку график любого многочлена называют параболой) формулируется так:
для функции f(x), заданной таблицей (2), найти многочлен Рп(х) такой, что выполняется совокупность условий интерполяции
(3)
Найти многочлен Рп(х) — это значит, учитывая его каноническую форму
(4)
найти его n+ 1 коэффициент. Для этого имеется как раз n+ 1 условие (3). Таким образом, чтобы многочлен (4) был интерполяционным для функции (2), нужно, чтобы его коэффициенты удовлетворяли системе уравнений

Из курса алгебры известно, что определитель этой линейной системы (так называемый определитель Вандермонда) отличен от нуля, т.е. решение этой системы существует и единственно. Однако практическое построение интерполяционного многочлена таким путем малоэффективно. Поэтому изберем другой путь.
Будем строить многочлен п-й степени Ln(x) в виде линейной комбинации многочленов п-й же степени li(x) (i = 0,1,..., п). Для того, чтобы такой многочлен был интерполяционным для функции f(x), достаточно зафиксировать в качестве коэффициентов сi, этой линейной комбинации заданные в табл.2 значения yi =f(xi), а от базисных многочленов li(x) потребовать выполнения условия
(5)
В таком случае для многочлена в каждом узле xj (j=0, 1, …, п), в силу (5), справедливо.

т.е. выполняются условия интерполяции (3).
Чтобы конкретизировать базисные многочлены li(x), учтем, что они должны удовлетворять условиям (5). Равенство нулю i-го многочлена во всех узлах, кроме i-го, означает, что li(x) можно записать в виде

а коэффициент Аi этого представления легко получается из содержащегося в (5) требования li(x)=1. Подставляя в выражение li(x) значение х=хi и приравнивая результат единице, получаем

Таким образом, базисные многочлены Лагранжа имеют вид

а искомый интерполяционный многочлен Лагранжа есть
(6)
Пример. Рассмотрим квадратичную интерполяцию функции y=sinx на отрезке [0, π/2] по ее трем значениям: sin0 = 0, sinπ/4=22, sinπ/2=1.
По формуле (8) строим многочлен Лагранжа второй степени

или после преобразований

Подставим в полученный интерполяционный многочлен контрольную точку Получим приближенное значение отличающееся от значения sinπ/6= 0.5 на величину ≈0.017.
Интерполяционная схема Эйткена
Пусть функция f(x) и расположение узлов x0, x1, ..., xn на отрезке интерполяции [a, b] таковы, что имеет место сходимость процесса интерполяции. Пусть требуется найти не общее выражение Ln(x), а лишь его значения при конкретных x, т.е. решается частная задача вычисления отдельных приближенных значений функции f(x) с помощью вычисления соответствующих им значений интерполяционного многочлена Лагранжа Ln(x). Для построения эффективного способа решения такой частной задачи интерполяции учтем следующее:
использование многочлена Лагранжа в виде (6) неудобно из-за его громоздкости, что ведет к большим вычислительным затратам;
заранее не известно, многочлены какой степени нужно использовать для интерполирования данной функции с требуемой точностью. А постепенное улучшение точности за счет вычислений Ln(x) с большим показателем степени n невыгодно, т.к. Ln-1(x) плохо перестраивается в Ln(x);
функция f(x) задается таблицей своих приближенных значений. Процесс сходимости Ln(x) к f(x) при больших значениях n будет нарушен влиянием на результат исходных ошибок.
Построим вычислительную схему для получения приближенного значения сеточной функции f(x) в заданной точке , в основу которой будет положена интерполяция Лагранжа на сетке узлов Организация вычислений по этой схеме будет иметь итерационный характер. Каждый шаг заключается в вычислении некоторого определителя второго порядка.
Пусть даны две точки на кривой Построим функцию

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

Она тоже является многочленом Лагранжа первой степени, построенным по двум точкам
Если на кривой y=f(x) заданы три точки (xo,yo), (x1, y1), (x2, y2), то, используя введенные линейные функции образуем новую функцию:

Эта функция есть многочлен второй степени, решающий задачу параболической интерполяции по трем точкам (x0, y0), (x1,y1), (x2,y2). Но такой многочлен единственный, следовательно, где - многочлен Лагранжа.
Продолжая процесс рассуждения, получим рекуррентное задание последовательности интерполяционных многочленов Лагранжа, которое составляет суть интерполяционной схемы Эйткена:
(7)
где i=1, 2,…, n и по определению P0(x)=y0, P1(x)=y1.
Схема Эйткена легко реализуется на ЭВМ. Организация вычислений по формуле (7) должна быть такова, что если заранее неизвестна степень интерполяционного многочлена, который нужно использовать для вычисления y(x), то должно происходить постепенное повышение степени k интерполирующих ее многочленов за счет подключения новых узлов. Счет ведется до тех пор, пока идет уточнение приближенного значения y(x).
Об этом можно судить по уменьшению величины при увеличении k и подходящем фиксировании i.
Пример. Пусть некоторая функция y=y(x) задана таблицей своих значений, округленных до двух знаков после запятой:

Рассмотрим процесс вычисления двух значений этой функции по схеме Эйткена в точках: а) б) Результаты промежуточных вычислений (в которых один знак запасной) сведем в табл. 3, 4. Числа в столбцах, помеченных посредством представляют собой значения многочленов Лагранжа k-ой степени, построенных по узлам от i-го до (i+k)-го рекуррентно по формуле:
(8)
где k=1, 2, …, в соответствии с интерполяционной схемой Эйткена.
Порядок получения этих значений показан проставленными в каждой клетке номерами.
Таблица 3.
Вычисление по схеме Эйткена значения y(0.1)

Таблица 4.
Вычисление по схеме Эйткена значения y(0.25)

Процесс вычисления значений многочленов Лагранжа ведется до тех пор, пока идет уточнение приближенного значения т.е. уменьшается величина ) при увеличении k и подходящем фиксировании i.
Например, для подсчета приближенного значения данной функции в точке расположенной между узлами xo=0.0 и xi=0.2, целесообразно в качестве основной последовательности значений многочленов Лагранжа брать строку таблицы 3, соответствующую значению i=0, т.е. числа
Вычислив разности между последующими и предыдущими числами этой строки, а именно:
0.0050.0040.010,
видим, что дальнейший счет бессмыслен; разность перестала уменьшаться. Т.е. по данной информации о функции y(x) более точное значение y(0.1), чем y(0.1)=1.001, получаемое с помощью L3(0.1), найти не удастся.
В случае б) вычисление по схеме Эйткена дает следующий результат: y(0.25)≈1.038, полученный с помощью L3(0.25).
Задания: выполнить задание 1 ИДЗ№2.
Занятие № 10. Интерполяционный многочлен Ньютона. Многочлены Чебышева.
Конечные разности
Зададимся целью придать интерполяционной формуле более простой вид, подобный виду широко используемой в математическом анализе формулы Тейлора. Если в интерполяционном многочлене Лагранжа (6) все слагаемые однотипны и играют одинаковую роль в образовании результата, хотелось бы иметь такое представление интерполяционного многочлена, в котором, как и в многочлене Тейлора, слагаемые располагались бы в порядке убывания их значимости. Такая структура интерполяционного многочлена позволила бы более просто перестраивать его степень, добавляя или отбрасывая удаленные от начала его записи члены.
Поставленной цели будем добиваться сначала для несколько суженной постановки задачи интерполяции. А именно, будем считать, что интерполируемая функция у = f(x) задана своими значениями на системе равноотстоящих узлов т.е. таких, что любой узел хi, этой сетки можно представить в виде
хi = х0 + ih,
где i=0, 1, …, п, a h>0 — некоторая постоянная величина, называемая шагом сетки (таблицы).
Прежде чем строить желаемые интерполяционные формулы, рассмотрим элементы теории конечных разностей.
Вычитая из каждого последующего члена конечной последовательности из n+1 чисел предыдущий, образуем п конечных разностей первого порядка

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

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

где к = 1, 2, …, п и
В некоторых случаях требуется знать выражения конечных разностей непосредственно через значения функции лежащей в их основе. Для нескольких первых порядков разностей их можно получить прямой подстановкой:

Подметив закономерность в коэффициентах рассмотренных представлений конечных разностей, записываем общую формулу
(1)
которая может быть строго обоснована методом математической индукции и которая напоминает биномиальное разложение для (у - 1)к.
Привлекая определение производной, можно обнаружить прямую связь между конечными разностями и производными. А именно, если учесть, что

то можно сказать, что при малых h имеет место приближенное равенство

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

т.е. и, вообще,
(2)
Таким образом, на конечные разности можно смотреть как на некоторый аналог производных. Отсюда справедливость многих их свойств, одинаковых со свойствами производных.
Отметим лишь простейшие свойства конечных разностей:
1)конечные разности постоянной равны нулю;2)постоянный множитель у функции можно выносить зазнак конечной разности.
3)конечная разность от суммы двух функций равна суммеих конечных разностей в одной и той же точке.
Свойства 2 и 3 характеризуют операцию взятия конечной разности как линейную операцию.
Учитывая роль, которую играют многочлены в теории интерполирования, посмотрим, что представляют собой конечные разности многочлена.
Поскольку многочлен в своей канонической форме есть линейная комбинация степенных функций, положим сначала у=хп. Используя биномиальное разложение п-й степени двучлена, получим:

т.е. первая конечная разность степенной функции у=хп есть многочлен степени п-1 со старшим членом Если взять теперь конечную разность от функции
(3)
то, в силу линейных свойств Δу, можно записать

Первое слагаемое в этой сумме, как выяснено, есть многочлен (n-1)-й степени, второе, аналогично, — многочлен степени п -2, и т.д. Следовательно, первая конечная разность многочлена (3) в точке х с шагом h есть тоже многочлен со старшим членом a0nhxn-1, вторая конечная разность — многочлен со старшим членом a0n(n-1)h2xn-2,…, к-я разность — многочлен со старшим членом а0п(п-1)...(n-к + 1)hkxn-k.
При к = п получаем постоянную разность п-го порядка

для многочлена (3); конечные разности более высоких порядков, естественно, равны нулю.
Итак, главный вывод из предыдущих рассуждений: п-е конечные разности многочлена п-й степени постоянны, а (п+1)-е и все последующие равны нулю.
Более важным для понимания сути полиномиального интерполирования является утверждение, обратное сделанному выше выводу. А именно, что если конечные разности п -го порядка некоторой функции у=у(х) постоянны в любой точке х при различных фиксированных шагах h, то эта функция у(х) есть многочлен степени п.
Для функции у = f(x)y заданной таблицей своих значений в узлах где конечные разности разных порядков удобно помещать в одну общую таблицу с узлами и значениями функции (последние можно интерпретировать как конечные разности нулевого порядка). Эту общую таблицу называют таблицей конечных разностей.

Пример. Составим таблицу конечных разностей для функции у =xln2x по ее значениям, вычисленным с тремя знаками после запятой в точках xj=0.4 + 0.2j, где j=0, 1, ..., 10.

Абсолютные величины конечных разностей сначала убывают с увеличением их порядка, а затем начинают увеличиваться. Это типичное поведение конечных разностей при ограниченной точности задания значений сеточной функции.
Природу наблюдаемого в примере поведения модулей конечных разностей нетрудно понять. Если шаг достаточно мал, а данная табличная функция — достаточно гладкая, то сначала с увеличением k в силу упомянутой связи (2). Когда эти величины становятся достаточно малыми, большую роль начинают играть продукты взаимодействия исходных ошибок округления (так называемый шум округлений).
Что происходит с одной отдельно взятой ошибкой величины ε у значения yi, можно проследить по таблице. Как видим, с ростом порядка разностей она «расползается» по таблице и увеличивается по абсолютной величине.

Погрешности, имеющиеся у каждого из данных значений функции, с ростом порядка разностей все больше взаимодействуют.
Из сделанных наблюдений напрашивается следующий вывод. Если какой-то столбец в таблице конечных разностей (в ее эксплуатируемой части) состоит из чисел, абсолютные величины которых составляют всего несколько единиц десятичного знака, являющегося последним в записи исходных значений функции, скажем, не превосходят величины 10ε, где ε — абсолютная погрешность исходных данных, то эти конечные разности и разности всех последующих порядков не несут практически никакой информации о функции, и их не следует использовать. Разности же предшествующего столбца называются практически постоянными, и их порядок определяет степень многочлена, которую можно и должно использовать для идеальной в данных условиях полиномиальной интерполяции.
Вспоминая о том, что многочлен k-й степени имеет k-е разности постоянными, а все последующие — нулевыми, приходим к заключению, что если k-е разности таблицы конечных разностей некоторой функции практически постоянны, то эта функция ведет себя в рассматриваемой области, как многочлен k-й степени; эту степень и следует применять для интерполирования с наибольшей для данных реалий точностью.Обратимся к таблице нашего примера. Видим, что если исключить из рассмотрения верхнюю диагональную строку, то для всей остальной части таблицы третьи разности удовлетворят условию (где 0.0005 — предельная абсолютная погрешность значений yi). В такой ситуации разности более высоких порядков не следовало вообще вычислять, а разности второго порядка можно считать практически постоянными, т.е. для подсчета любых промежуточных значений данной функции, за исключением, быть может, тех, которые находятся вблизи узла x0= 0.4, нужно применять квадратичную интерполяцию.
Конечноразностные интерполяционные формулы
Пусть функция у=f(x) задана на сетке равноотстоящих узлов хi=х0+ ih, где i= 0, 1, …, п, и для нее построена таблица конечных разностей.
Будем строить интерполяционный многочлен Pn(х) в форме
(4)
Его п+1 коэффициент а0, а1, ….,ап будем находить последовательно из п+1 интерполяционных равенств

А именно, полагая i=0, т.е. х=х0 в (4) имеем а по условию интерполяции Рп(х0)=у0, следовательно, ао =у0.
Далее, при i = 1 аналогично получаем равенство

в которое подставляем уже найденное значение Разрешая это равенство относительно а1 и используя обозначение конечной разности, получаем

Следующий шаг, при i = 2, дает:

Полной индукцией можно показать справедливость выражения
(5)
Подставляя найденные коэффициенты (4), получаем многочлен
(6)
который называют первым интерполяционным многочленом Ньютона.
Учитывая, что каждое слагаемое многочлена (6), начиная со второго, содержит множитель x-x0, естественно предположить, что этот многочлен наиболее приспособлен для интерполирования в окрестности узла х0. Будем называть узел x0 базовым для многочлена (6), и упростим (6) введением новой переменной q равенством или (что то же) равенством х=x0+qh. Так как h при любых i=0, 1, ..., n

то в результате подстановки этих разностей в (6) приходим к первой интерполяционной формуле Ньютона в виде
(7)
где обозначениеуказывает не только на п-ю степень многочлена, но и на базовый узел x0 и связь переменных х и q.
Первая формула Ньютона (7) обычно применяется при значениях |q|<1, а именно, для интерполирования вперед, (при т.е. при ) и экстраполирования назад (при х < х0, т.е. при q < 0).
Так как, реально степени интерполяционных многочленов бывают не так велики, в то время как таблицы значений функций достаточно обширны, и так как в реальной числовой таблице никаких индексов — номеров узлов нет, то за базовый для формулы (7) узел x0 можно принимать узел, ближайший к заданной фиксированной точке х, если за ним имеется достаточное число узлов для построения необходимых для (7) разностей. Поскольку в первой формуле Ньютона используются нисходящие диагонали таблицы конечных разностей, то такое смещение узла, принимаемого за базовый, в конце таблицы будет неприемлемо.Учет этого обстоятельства приводит к потребности в симметричной, в определенном смысле, для (7) формулы, которая была бы пригодной для интерполирования в конце таблицы. Для этого, в отличие от (4), форма интерполяционного многочлена Рп(х) берется такой, которая предусматривает поочередное подключение узлов в обратном порядке: сначала последний, потом предпоследний и т.д., т.е.

Коэффициенты этого многочлена находятся аналогично тому, как они находились для многочлена (4), только здесь подстановка узловых точек вместо х и рассмотрение интерполяционных равенств производится тоже в обратном порядке. Полагая имеем:

и т.д. В общем случае

Таким образом, получаем второй интерполяционный многочлен Ньютона

(8)
в котором базовым является узел хп и коэффициенты которого определяются конечными разностями, расположенными на восходящей от уп диагонали.
Положим в (8) х=хп+qh, иначе, введем новую переменную и преобразуем к ней входящие в (8) разности:

В результате приходим ко второй интерполяционной формуле Ньютона вида
(9)
Ее также целесообразно использовать при значениях |q|<1 т.е. в окрестности узла хn для интерполирования назад (при q∈(-1, 0)) н экстраполирования вперед (при q> 0).
Наряду с выведенными специально для начала и конца таблицы первой и второй интерполяционными формулами Ньютона, имеется еще несколько формул, рассчитанных на их применение в центральной части таблицы и потому называемых центральными интерполяционными формулами. Прежде, чем определять эти формулы, введем понятие центральных разностей.
Будем считать, что узел x0 расположен в середине таблицы, и нумерация остальных узлов производится, начинаясь с х0, с использованием как положительных, так и отрицательных индексов, т.е. считаем хi=х0+ih, где i= 0, ±1, ±2,.... Тогда центральная часть таблицы конечных разностей будет проиндексирована так, как это показано в таблице.
Все подчеркнутые в ней конечные разности (находящиеся с в одной строке и на полстроки выше и ниже) называются центральными разностями.

Интерполяционный многочлен ищем в форме
(10)
предполагающей постепенное подключение узлов хi: сначала при i=0, затем при i=1, потом при i= -1 и т.д., т.е. с двух сторон от х0. При этом здесь и далее не будем фиксировать степени многочленов и не будем стремиться выписывать общие и, тем более, последние члены таких многочленов. Как и в предыдущих случаях, коэффициенты находим один за другим последовательной подстановкой в Р(х) и в интерполяционные равенства значений

и т.д. Введя новую переменную и выразив через нее разности для всех i = 0, ± 1, ± 2, ..., в результате подстановки этих разностей и выражений коэффициентов в шаблон (10), приходим к первой интерполяционной формуле Гаусса
(11)
Записанные слагаемые легко дополнить следующими, если знать, что в этой формуле используются нижние центральные разности все возрастающих порядков, т.е. те, которые подчеркнуты в табл.1 сплошной чертой.
Таблица 1

Совершенно аналогично, подключая узлы в другом порядке после x0 сначала предшествующий, затем последующий и т.д., т.е. можно построить вторую интерполяционную формулу Гаусса
(12)
использующую верхние центральные разности (подчеркнутые в табл.1 пунктирной линией).
Интерполяционные формулы Гаусса служат полуфабрикатами для получения более симметричных, использующих все центральные разности интерполяционных формул.
Так, полусумма первого и второго интерполяционных многочленов Гаусса после преобразований приводит к формуле
(13)
называемой интерполяционной формулой Стирлинга К.
Если же взять полусумму второго интерполяционного многочлена Гаусса и такого же многочлена, но с нижними индексами, увеличенными на единицу (т.е. с базовой точкой х1 вместо x0), то придем к интерполяционной формуле Бесселя
(14)
В последней формуле обращает на себя внимание тот факт, что она сильно упростится, если в нее подставить значение q=1/2, соответствующее значению аргумента
Этот частный случай формулы Бесселя называют формулой интерполирования на середину:
(15)
Итак, если точка х, в которой нужно найти приближенное значение таблично заданной функции f(х), находится в начале или в конце таблицы, применяется соответственно первая (7) или вторая (9) формулы Ньютона с таким выбором базовой точки, чтобы значение |q| было как можно меньше. Если точка х находится в середине таблицы, то всегда можно зафиксировать точку x0 в таблице центральных разностей так, чтобы либо было по модулю меньше 0.25 и тогда применять интерполяционную формулу Стирлинга (13), либо чтобы q ∈[0.25, 0.75] и использовать формулу Бесселя (14).
Пример. Пусть требуется для функции у = f(x), заданной в предыдущем примере таблицей нескольких своих значений с тремя знаками после запятой, найти приближенные значения: а) f(0.5); б) f(1.22); в) f(0-5); г) f(1.94); д) f(2.5), записав предварительно соответствующие каждому случаю интерполяционные формулы.
Для решения поставленной задачи учитываем, что значения f(x) заданы на сетке равноотстоящих узлов, поэтому здесь можно применить конечноразностную интерполяцию. При этом будем пользоваться уже составленной таблицей конечных разностей и проведенным ранее ее анализом на выявление оптимальной степени многочлена. Для случаев б-д фиксируем вторую степень, для а — третью. В каждом случае, т.е. для конкретного значения аргумента, выбираем базовый узел, подсчитываем значение вспомогательной переменной q и, в зависимости от положения базового узла и значения q, пользуясь представленными в таблице числами, записываем требуемую интерполяционную формулу. Подстановка в нее значения q приводит к искомому значению
а) При x=0.5 (начало таблицы) полагаем х0=0.4; тогда Соответствующую интерполяционную формулу для аппроксимации f(x) при х=0.4 + 0.2q с q∈(-1,1) записываем, глядя на первую интерполяционную формулу Ньютона (7) и таблицу значений:

Отсюда получаем искомое значение

б)Точка х=1.22 находится в средней части таблицы. Поэтому здесьцелесообразно применить формулу Стирлинга или Бесселя. Полагая и найдя останавливаемся на формуле Стирлинга (13), которая в данном случае имеет вид

и при q = 0.1 приводит к искомому значению

в)Здесь, очевидно, напрашивается применение формулы (15)интерполирования на середину. Полагая имеем:

г)Глядя на положение точки х=1.94 в заданной системе узлов, видим, что для вычисления f(1.94) также возможно применениецентральных интерполяционных формул. Положив х0=1.8 и вычислив на основе (14) записываем интерполяционную формулу Бесселя

Из нее получаем

д)Точка х = 2.5 расположена за последним узлом, поэтому для экстраполяции f(x) здесь однозначно следует применить вторую интерполяционную формулу Ньютона (9). Считая хп =2.4 (индекс п здесь ис-пользуется условно, без придания ему конкретного значения), записываемформулу экстраполяции

откуда при находим

Определение и свойства многочленов Чебышева
Многочленом Чебышева называется функция
(15)
где
Функция Тп(х), представленная с помощью тригонометрических функций, на самом деле является многочленом при любом n=0, 1, 2,.... и справедливо равенство
(16)
Формула (16) определяет при n=1, 2, 3,... последовательность функций начинающуюся с Т0(х)=1, T1(x) = x, рекуррентно; при этом нужно иметь в виду, что здесь x∈[-l, 1], как и в (15).
Подставляя в (16) заданные начальные члены последовательности (Тn(х)), найдем несколько ее последующих членов:
и т. д.
Графики нескольких многочленов Чебышева (с первого по четвертый) изображены на рис. 1.

Рис. 1. Графики многочленов
Анализ рекуррентной формулы (16) позволяет считать очевидными следующие факты:
1)все функции Тn(х), определенные в (15), являются многочленами при любом натуральном n;
2)степени этих многочленов возрастают с увеличением n, причем старший член многочлена Тn(х) равен 2n-1 хn;
3) многочлены Тn(х) при четных n выражаются через степенные функции только четных степеней, при нечетных — только нечетных.
Наряду с многочленами Чебышева Тn(х), часто используют многочлены, получаемые из Тn(х) делением на старший коэффициент, т.е.

— многочлены со старшим коэффициентом 1. Будем называть их нормированными многочленами Чебышева.
Многочлены Чебышева обладают рядом замечательных свойств. Рассмотрим некоторые их свойства, имеющие отношение к поставленной выше проблеме аппроксимации функций.
1. Многочлен Чебышева Тn(х) имеет на отрезке [-1, 1] ровно n различных действительных корней; все они задаются формулой
(17)
2. Корни многочленов Чебышева перемежаются с точками их наибольших и наименьших значений, равных соответственно +1 и -1 и для . а именно, при j=0, 1, …, n имеют место экстремумы
(18)
3. (теорема Чебышева). Из всех многочленов степени n со старшим коэффициентом 1 нормированный многочлен Чебышева наименее уклоняется от нуля на отрезке [-1, 1].
Это свойство означает, что среди всех многочленов степени n вида
(19)
именно нормированный многочлен Чебышева минимизирует максимальное расстояние от графика многочлена при х∈[-1,1] до оси абсцисс.
Интерполяция по Чебышевским узлам
Желая, чтобы интерполяционный многочлен Лагранжа Ln(х)в целом хорошо приближал функцию у= f(x) на отрезке [а, b], поставим вопрос: как расположить на нем n+1узлов интерполяции хi, (i= 0, 1, …, n), чтобы при этом минимизировать максимальную на [а,b] погрешность?
Максимальная погрешность интерполирования достаточно гладкой функции на отрезке [-1, 1] многочленом n-й степени будет минимальной, когда в качестве узлов интерполяции берутся корни многочлена Чебышева Tn+1(t). Будем называть их чебышевскими узлами интерполяции.
Знание экстремальных значений многочлена Чебышева позволяет уточнить величину максимального отклонения Ln(x) от f(x) при таком выборе узлов, т.е. когда точки ti есть корни a А именно,

Эта оценка называется наилучшей равномерной оценкой погрешности интерполяции.
Если функция f(x) бесконечно дифференцируема на [а, b] и в качестве узлов интерполяции берутся корни многочленов Чебышева (приведенные к отрезку [а, b], то
Обобщением этого факта для непрерывных (необязательно дифференцируемых) функций и произвольных (не обязательно интерполяционных) многочленов является широко известная в математическом анализе теорема.
Теорема (Вейерштрасса). Для любой непрерывной на [а, b] функции f(x) найдется многочлен Qn(x) такой, что

Доказано, что для любой функции f(x) существует единственный многочлен такой, который из всех многочленов Qn(x) степени n наилучшим образом аппроксимирует на [а, b] функцию f(x), минимизируя максимальное расстояние между f(x) и Qn(x). Этот многочлен, т.е. многочлен такой, что

называется многочленом наилучшего равномерного приближения для f(x) на [а, b] или ее чебышевским приближением.
Одним из характеристических свойств многочленов наилучших равномерных приближений является критерий Чебышева. Его отражает следующая теорема.
Теорема (Чебышева). Многочлен является многочленом наилучшего равномерного приближения для функции f(x) тогда и только тогда, когда на [а, b] существует не менее n + 2 точек хi таких, что в них поочередно принимаются наибольшие положительные и отрицательные отклонения, т. е. поочередно разность f (xi) -равна Е или - Е, где

Эта теорема говорит о том, что максимальная ошибка аппроксимации функции многочленом наилучшего равномерного приближения реализуется в числе точек, большем, по меньшей мере, на 2, чем степень многочлена, причем знаки ошибки чередуются. Точки хi, в которых реализуется максимальное отклонение многочлена от f(x) на [а, b], называются точками чебышевского альтернанса.
К сожалению, неизвестны ни общий вид многочленов наилучших равномерных приближений, ни способы их построения, имеются лишь некоторые методики построения многочленов, близких к наилучшим равномерным, а также способы построения чебышевских приближений невысокого порядка для нескольких весьма узких классов функций. Последние существенно опираются на приведенную теорему о чебышевском альтернансе, что демонстрируется в следующих двух простейших случаях.
Случай А. Пусть функция f(х) непрерывна на [a,b], и пусть для нее требуется построить многочлен наилучшего равномерного приближения нулевой степени. Обозначим этот приближающий многочлен через Он определяется всего одним параметром: . Чтобы найти значение этого параметра для заданной функции f(х), воспользуемся тем свойством непрерывной на замкнутом промежутке функции, согласно которому на нем всегда найдутся, по крайней мере, две точки, в которых она принимает свои наименьшее и наибольшее значения.
Пусть Тогда совершенно очевидно, что полагая т.е. подменяя функцию f(x) функцией будем иметь максимальное отклонение

причем точки отрезка [а, b], в которых оно реализуется — это точки, где принимаются значения m и М. В силу непрерывности f(х), локальные минимумы и максимумы должны чередоваться; по меньшей мере, два из них определяют точки чебышевского альтернанса (с чередованием знаков разностей f(x)-m+M2 ). Поэтому не существует другой постоянной, которая приближала бы f(x) на [а, и] лучше, чем постоянная m+M2.

Случай Б. Пусть аппроксимируемая функция f(x) дифференцируема и выпукла (в широком смысле) на отрезке [а, b], a аппроксимирующая ее функция — многочлен наилучшего равномерного приближения первой степени. Чтобы найти его коэффициенты следует изучить разность между f(x) и φ(х), т. е. функцию
Так как функция f(x) по предположению выпукла, а сдвиг на линейную функцию A0 + А1х не изменяет выпуклости, то и функция u(х) выпукла на [а, b]. Следовательно, существует единственная точка с∈[а, b], в которой u(х) имеет минимум; если u(х) выпукла вверх, то в точке с должен быть максимум u(х). В любом случае, точка с∈[а, b] есть точка экстремума u(х), и за счет возможности варьирования коэффициентов функции φ(х) (точнее, коэффициента А1) можно считать, что точка с является внутренней точкой отрезка [а, b].
Потребуем, чтобы точки а, с и b в указанной последовательности составляли чебышевский альтернанс, т. е. чтобы в них последовательно принимались значения Е, - Е, Е или - Е, E, -E, где — максимальная погрешность аппроксимации функции f(x) функцией φ(х). Добавляя к этим требованиям еще необходимое условие экстремума дифференцируемой функции u(х) в точке с и возвращаясь к исходным функциям, приходим к системе четырех уравнений относительно четырех неизвестных A0, А1, Е и с (из которых, в основном, лишь первые три неизвестные представляют интерес):

Получить решение такой системы в общем виде не представляется возможным, поскольку неизвестная величина с входит в нее нелинейным образом (в каждом конкретном случае подобная система без проблем решается численно).
Пример. Построим многочлены наилучшего равномерного приближения нулевой и первой степеней для функции y = sinx на отрезке [0, π/6].
Сразу заметим, что данная функция всюду дифференцируема (а значит, непрерывна) и выпукла вверх на заданном отрезке. При этом

Следовательно, согласно рассмотренному выше случаю А, найдя

при х∈[0, π/6] можно считать, что sinх≈ 0.25 с предельной погрешностью 0.25.
Далее, в соответствии со случаем Б, продифференцировав данную функцию, составляем систему:

Из нее последовательно находим:

Таким образом, функцию y = sinx на отрезке [0, π/6] можно подменить линейной функцией у = 0.0045 + 0.9549x, и наибольшая ошибка при этом не будет превышать величины ≈0.0045.
Задания: выполнить задание 2 и 3 ИДЗ№2.
Занятие № 11. Интерполяция с кратными узлами и формула Эрмита.Рассмотрим задачу полиномиальной интерполяции функции у= f(x) в более общей постановке.
Пусть на промежутке расположены m+1 несовпадающих узлов и пусть в этих точках известны значения у0 = f(x0), у1=f(x1,),..., ут = f(xm) данной функции, а также некоторые ее производные (максимальный порядок производных в разных узлах различен; в каких-то узлах производные могут быть вовсе неизвестны). Такие узлы будем называть кратными узлами. Конкретнее, будем считать, что заданы:
(1)
тогда кратность узла x0 считается равной k0, узла x1 - k1, …, узла xm - km.
Предполагая, что суммарная кратность узлов есть
(2)
ставим задачу построения многочлена Нп(х) степени n (не выше п) такого, что
(3)
где — заданные посредством (1) значения функции f(x) и ее производных и по определению считается Многочлен будем называть интерполяционным многочленом Эрмита, а совокупность требований (3) — условиями эрмитовой интерполяции.
Формально можно считать, что нахождение такого многочлена состоит в том, чтобы однозначно определить п +1 коэффициентов a0, а1,..., ап его канонического представления
(4)
из условий (3). В силу предположения (2) о суммарной кратности узлов эрмитовой интерполяции, совокупность требований (3) можно рассматривать как систему из п+1 уравнений относительно п+1 неизвестных — коэффициентов ak многочлена (4):

Выявление общего вида интерполяционных многочленов Эрмита Нп(х) представляет непростую задачу и требует привлечения определенных сведений из теории функций комплексной переменной. Рассмотрим одну из возможных процедур фактического построения таких многочленов, не требующую знания их общего вида.
Пусть Lm(x) — интерполяционный многочлен Лагранжа, построенный по данным т+1 значениям yi = f(x i = 0,1,..., т. Будем пользоваться обозначением Так как по условию т заведомо не превосходит п, то по теореме о делении многочлена с остатком искомый многочлен Эрмита Нп(х) можно представить в виде
(5)
где Hn-(m+1)(x) — некоторый неизвестный пока многочлен степени п-т-1.
Для построения многочлена Hn-(m+1)(x) будем привлекать информацию о производных данной функции, т.е. равенства в тех узлах хi где первые производные, в соответствии с (1), заданы (информация о самих значениях функции уже полностью исчерпана: в силу для всех хi от х0 до хm, согласно (5), будет и при любых i∈{ 0,1,..., т}).
Продифференцировав равенство (5), имеем
(6)
Поскольку в тех узлах хi, где по условию эрмитовой интерполяции справедливо можно записать

Отсюда выражаем значения многочлена Hn-(m+1)(x) в этих узлах:

Правая часть этого равенства может быть вычислена; обозначим ее через . Таким образом, в ряде узлов хi известны значения многочлена Hn-(m+1)(x)= , по которым этот многочлен однозначно восстанавливается обычной лагранжевой интерполяцией, если в условиях (1) не содержится производных порядка, выше первого (т.е. нет ни одного узла кратности больше 1); подстановка найденного многочлена Hn-(m+1)(x)в (5) приводит к искомому интерполяционному многочлену Эрмита. Если же в исходной информации (1) об f(x) имеются значения производных более высокого порядка, чем первый, то для восстановления многочлена Hn-(m+1)(x) ставится задача эрмитовой же интерполяции, для чего наряду с полученными его значениями , находят значения его производных путем дифференцирования равенства (6) (возможно неоднократного, в зависимости от максимального порядка заданных производных функции f(x)). Эта процедура построения интерполяционных многочленов Эрмита все более низких степеней продолжается до исчерпывания всей информации (1) о функции и ее производных.
Рассмотрим реализацию описанного процесса эрмитовой интерполяции на простом примере, демонстрирующем возможность восстановления многочлена n-й степени по его значениям и значениям некоторых его производных при суммарной кратности узлов п+1.
Пример. Пусть сведения о некоторой функции у= f(x) представлены следующей дискретной информацией:

В соответствии с обозначениями (1) здесь: т=2; к0-1 =1, k1 -1 = 2, к2-1=1 => п+1=к0+к1+к2=7 => n=6. Таким образом, по данным сведениям о функции у = f(x), сосредоточенным в трех узлах кратности, соответственно, 2, 3, 2, следует строить интерполяционный многочлен Эрмита Н6( х).
Согласно предложенной выше схеме, сначала, пользуясь столбцами таблицы данных, записываем интерполяционный многочлен Лагранжа второй степени

Далее по формуле (5) представляем Н6(х) через L2(x),H3(x) и П3(х):
(7)
и дифференцируем этот многочлен дважды:

Подстановкой в значений х=-1, х=0 и х=1 и приравниванием заданным значениям получаем значения

Учитывая их, из условия Я находим
Итак, для выявления многочлена H3(х) в (7) снова имеем задачуэрмитовой интерполяции с данными, содержащимися в следующейтаблице:

Здесь: В соответствии с (5), для этого случая записываем:

(где — многочлен Лагранжа, интерполирующий функцию Н3(х))-Остается найти постоянную H0, для чего воспользуемся условием Имеем:

Следовательно,
Подставив это в (7), получаем окончательное выражение искомого многочлена:
Занятие № 12. Понятие сплайна. Виды сплайнов.Цель - ознакомить студентов с понятием сплайн, видами сплайнов.
Помимо «классических» интерполяционных полиномов, в XX веке получила широкое распространение в практике еще одна разновидность интерполянтов, называемых сплайнами. Исторически они возникли в авиастроении при конструировании обтекаемых профилей.
Сплайном S(x) является интерполирующая функция, которая обладает наименьшей кривизной и имеет квадратично интегрируемую вторую производную. Так как сплайн является интерполянтом, то для него должно выполняться условие (2.1):
(1)
Вид функции сплайна S(x) находится из условия минимальности функционала
(2)
Эта функция физически моделируется гибким и упругим стержнем, который проходит через все заданные точки интерполяции. Минимум функционала (2) соответствует минимуму упругой потенциальной энергии этого стержня.
Особенно часто используется естественный кубический сплайн - полином 3-й степени, удовлетворяющий требованиям (1) и (2). Такой естественный кубический сплайн является непрерывной функцией и имеет непрерывные первую и вторую производную на интервале интерполяции (x0, xn).
Как и в предыдущих параграфах данной главы, интерполянт строится на базе таблицы данных типа

Таблица с n+1 узлами определяет n подинтервалов между этими узлами. Кубический полином (естественный сплайн) необходимо построить для каждого подинтервала. Таким образом, интерполяция кубическими сплайнами относится к классу кусочно-многочленной интерполяции. Так как каждый кубический полином характеризуется четырьмя константами, поэтому для построения сплайна S(x) на всем интервале [x0; xn] необходимо определить 4n параметров.
Требование непрерывности функции, ее первой и второй производных во внутренних (n-1) узлах дает 3(n-1) уравнения для определения неизвестных параметров. Условие (1) для всех (n+1) узлов дает еще (n+1) уравнений. Недостающие два условия для однозначного определения естественного сплайна задаются приравниванием нулю вторых производных в точках x0 и xn
(3)
что означает нулевую кривизну сплайна на концах интервала.
Иногда при построении сплайнов задаются другие граничные условия на концах интервала интерполяции для самих функций или первых производных.
После формулировки общих положений займемся конструированием явного вида коэффициентов кубического сплайна. Для этого вначале рассмотрим отдельный подинтервал (xi, xi+1). Введем обозначения
(4)
(5)
(6)
Сплайн на подинтервале (xi, xi+1) удобно представить формулой:
(7)
где σi+1 и σi - константы, которые необходимо выразить через известные величины xi, yi. Очевидно, что функция (7) - полином 3-й степени относительно аргумента x.
Легко убедиться непосредственной подстановкой, что
(8)
Получим выражения для трех производных функции (7):
(9)
(10)
(11)
Формула (9) позволяет получить выражение для первой производной на концах любого подинтервала. Декларированная непрерывность первой производной сплайна S(x) на всем интервале (x0 , xn) позволяет приравнять правые и левые производные в точках узлов xi. В результате получается соотношение
(12)
В уравнении (12) и последующих этого параграфа использовано новое обозначение:
(13)
Заметим, что в этом параграфе величины определенные формулой (13), называются разделенными разностями, в отличие от конечных разностей, введенных в раннее.
Соотношения (12) получены для всех n-1 внутренних узлов. Индекс i пробегает значения 1, ... , n-1.
Таким образом, получена система (n-1) линейных уравнений относительно неизвестных величин σi.
В качестве недостающих условий можно взять уравнения, выражающие требования (3). Можно поступить иначе - использовать единственные кубические полиномы, которые проходят точно через 4 первые и 4 последние точки из заданных xi, yi (i = о,. . ., n). Тогда два дополнительных условия на концах интервала (x0 , xn) выразятся следующим образом:
(14)
где
(15)
Полученная система линейных уравнений (12) и (14) для нахождения n+1 параметров σi ( i = 0, ..., n) имеет единственное решение. Матрица коэффициентов при неизвестных σi может быть записана в следующем виде:
(16)
Правые части уравнений (свободные члены) соответственно равны:

Решение полученной системы линейных уравнений удобно реализовать, пользуясь трехдиагональностью матрицы коэффициентов при неизвестных, применяя метод прогонки. В результате искомые величины σi представятся в виде рекуррентных выражений:
(17)
где коэффициенты выражаются так:
(18)
(19)
При необходимости многократно вычислять сплайн его удобнее представить в несколько иной форме:
(20)
где xi < x< xi+1, i = 0, …, n - 1.
Преобразование функции (7) к форме (20) дает следующие выражения для искомых коэффициентов:
(21)
где i = 0,…, n - 1.
Подстановка коэффициентов (21) в выражение (20) дает функцию сплайна в явном виде для аргументов xi<x<xi+1, лежащих в i-м поддиапазоне интерполяции.
Пример. Пусть необходимо построить сплайн-интерполянт по набору шести точек.

Для вычисления коэффициентов сплайна целесообразно по формулам (4) - (6), (11) и (13) вычислить промежуточные данные и свести их в таблицу:

Пользуясь полученными величинами можно составить систему линейных уравнений (12 и (14). Запишем эту систему в матричном виде:
(22)
Решением системы (22) является следующий набор значений

которые мы представили в виде обыкновенных дробей.
Вычисленные значения σi и hi (i = 0, 1,. . 5) подставляются в формулы (4) - (7) или (20), (21), по которым можно получить явный вид сплайна S(x) на каждом i-м поддиапазоне xi < x < xi+1.
Например, на интервале х∈[3, 6] (т.е. между узлами i=1 и i=2) сплайн-интерполянт может быть выражен аналитически следующим кубическим полиномом:
S12(x) = 1,6102 + 0,0226 х + 0,3701 х2 - 0,0374 х3 ,
где коэффициенты вычислены с точностью до 0,0001. График полученного сплайна представлен на рис. 1.
Аналогичным образом можно построить кубические интерполирующие полиномы на всех интервалах, ограниченных узлами интерполяции. В результате получается сплайн-интерполянт для диапазона 1 < x < 9, график которого изображен на рис. 2.

Рис. 1. График сплайн-интерполянта на интервале х∈[3, 6].
Кружки - исходные данные, линия - график полинома (23)

Рис. 2. График сплайн-интерполянта на интервале х∈ [1, 9].
Задания: выполнить задание 4 ИДЗ№2.
Занятие № 13. Численное дифференцирование.Численное дифференцирование применяется в тех случаях, когда: функция f(x) задана таблично и, следовательно, методы дифференциального исчисления неприменимы; аналитическое выражение f(x) столь сложно, что вычисления производной представляют значительны трудности.
В основе численного дифференцирования лежит следующий прием: исходная функция f(x) заменяется на рассматриваемом отрезке [a,b] интерполяционным полиномом Pn(x) и считается, что f'(x) и P'n(x) примерно равны, т.е. f’(x)=P'n(x), .(a < x < b ).
Рассматривается следующая задача: дана таблица функции, требуется построить таблицу ее производной с тем же шагом.
Всегда, когда это возможно, для численного дифференцирования используется интерполяционный многочлен с равноотстоящими узлами, так как это значительно упрощает формулы численного дифференцирования. При равноотстоящих узлах в начале строится интерполяционный полином Ньютона, а затем он дифференцируется.
Если f(x) задана таблицей с неравноотстоящими узлами, то вначале строится интерполяционный полином Лагранжа, а затем он дифференцируется.
Будем изучать численное дифференцирование, основанное на первой интерполяционной формуле Ньютона. Из-за больших погрешностей численное дифференцирование практически используется для вычисления производных не выше второго порядка. Обычно при численном дифференцировании интерполяционный полином Ньютона строится не по всем узлам таблицы, а по трем-пяти узлам, близлежащим к точке, в которой требуется вычислить производную. Если требуется вычислить производную во всех узлах, то вначале полином Ньютона строится по первым 3-5 узлам и в них вычисляется производная, потом полином Ньютона строится по следующим 3-5 узлам и в них вычисляется производная и т. д. до тех пор, пока не будет просчитана вся таблица.
Если по узлам x0, xb x2, x3, x4 построить первый интерполяционный полином Ньютона и его продифференцировать с учетом

то получим приближенное выражение для в виде:
(1)
Вторая производная f"(x) в результате дифференцирования формулы (1) с учетом равенства

имеет следующее приближенное выражение:
(2)
Если требуется вычислить f'(x) в узлах x0, x1, …, xn, то в формулу (1) необходимо поочередно поставить q=0, q=1, q=2, q=3, q=4 соответственно.
В случае выбора узлов xj, xj+1, …, xj+4 вместо узлов x0, x1, ...,x4 в формулах (1), (2) следует заменить конечные разности y0 на конечные yj, которые берутся из строчки с номером j таблицы конечных разностей. Если требуется найти производные функции f(x) в основных табличных точка xj то каждое табличное значение считается за начальное и тогда xj=x0, q=0, а формулы (1), (2) будут соответственно для всех f'(x) и f'' (x) иметь вид: (1')
(2')
Если интерполяционный полином строится не по пяти, а по четырем или трем узлам, то в формулах (1), (1'), (2), (2') отбрасываются соответственно одно или два последних слагаемых.
Чтобы вычислить значение производной в точке, соответствующей концу таблицы, следует воспользоваться формулой, полученной дифференцированием второй интерполяционной формулы Ньютона, т.е. формулой

Метод неопределенных коэффициентов
Аналогичные формулы можно получить и для случая произвольного расположения узлов. Использование многочлена Лагранжа в этом случае приводит к вычислению громоздких выражении, поэтому удобнее применять метод неопределенных коэффициентов. Он заключается в следующем. Искомое выражение для производной k-го порядка в некоторой точке х=хi представляется в виде линейной комбинации заданных значений функции в узлах х0, x1, ..хn:
(3)
Предполагается, что эта формула имеет место для многочленов y=1, y=x-xi, у=(х-xi)n. Подставляя последовательно эти выражения в (3), получаем систему n+1 линейных алгебраических уравнений для определения неизвестных коэффициентов с0, с1, …, сn.
Пример. Найти выражение для производной ух в случае четырех равноотстоящих узлов (n = 3).
Равенство (3) запишется в виде
(4)
Используем следующие многочлены:
(5)
Вычислим их производные:
(6)
Подставляем последовательно соотношения (5) и (6) соответственно в правую и левую части (4) при х=х1:

Получаем окончательно систему уравнений в виде

Решая эту систему получаем
Подставляя эти значения в равенство (4), находим выражение для производной:

Задания: выполнить задание 1 и 2 ИДЗ№3.
Занятие № 14. Численное интегрирование. Квадратурные формулы Ньютона-Котеса.Задача численного интегрирования. Квадратурные формулы прямоугольников
Пусть требуется найти значение I интеграла Римана для некоторой заданной на отрезке [а, b] функции f(x). Хорошо известно, что для функций, допускающих на промежутке [а, b] конечное число точек разрыва первого рода, такое значение существует, единственно и может быть формально получено по определению:
(1)
где — произвольная упорядоченная система точек отрезка [а, b] такая, что

а ξi — произвольная точка элементарного промежутка [хi-1,_ хi].
В математическом анализе обосновывается аналитический способ нахождения значения I с помощью знаменитой формулы Ньютона-Лейбница
I = F(b)-F(a), (2)
где F(x) — некоторая первообразная для данной функции f(x).
К сожалению, применение этого весьма привлекательного подхода к вычислению I наталкивается на несколько серьезных препятствий. Самое главное из них — это несуществование первообразной среди элементарных функций для большинства элементарных функций f(x); например, таким способом не удается вычислить
Если первообразная F(x) для заданной функции f(x) все же найдена, то вычисление двух ее значений F(a) и F(b) может оказаться более трудоемким, чем вычисление существенно большего количества значений f(x).
Поскольку в общем случае значения функций находятся лишь приближенно, использование точной формулы (2) приводит к приближенному результату, который может быть более эффективно получен с помощью какой-либо специальной приближенной формулы на основе значений подынтегральной функции f(x). Такие специальные приближенные формулы для вычисления определенных интегралов называют квадратурными формулами (механическими квадратурами) или формулами численного интегрирования. Первый из этих терминов можно связать с геометрическим смыслом определенного интеграла: вычисление при f(х)>0 равносильно построению квадрата, равновеликого криволинейной трапеции с основанием [а,b] и «крышей» f(x).
Простые квадратурные формулы можно вывести непосредственно из определения интеграла, т.е. из представления (1). Зафиксировав там некоторое п≥ 1, будем иметь
(3)
Это приближенное равенство назовем общей формулой прямоугольников (площадь криволинейной трапеции приближенно заменяется площадью ступенчатой фигуры, составленной из прямоугольников, основаниями которых служат отрезки [хi-1, xi], а высотами — ординаты f(ξi), рис. 1).

Рис. 1. Геометрическая интерпретация общей формулы прямоугольников (.3)
Чтобы из общей формулы (3) получить конструктивное правило приближенного вычисления интеграла, воспользуемся свободой расположения точек х,, разбивающих промежуток интегрирования [а, b] на элементарные отрезки [xi-1, xi], и свободой выбора точек ξi на этих отрезках.
Условимся в дальнейшем пользоваться равномерным разбиением отрезка [а, b] на п частей точками хi с шагом полагая
(4)
При таком разбиении формула (3) приобретает вид
(5)
Теперь дело за фиксированием точек ξi на элементарных отрезах [xi-1, xi]. Рассмотрим три случая.
1)Положим ξi = xi-1. Тогда из (5) получаем
(6)
2)Пусть в (5) ξi = xi. Тогда имеем
(7)
Формулы (6) и (7) называются квадратурными формулами левых и правых прямоугольников соответственно. Совершенно очевидно (рис. 2), что дают двусторонние приближения к значению I интеграла от монотонной функции.

Рис. .2. Геометрическое оценивание интеграла от монотонной функции с помощью:
Можно рассчитывать на большую точность получения значения интеграла, если взять точку ξi посередине между точками xi-1 и xi. Отсюда приходим к следующему случаю.
3) Фиксируем В результате имеем квадратурную формулу средних прямоугольников или, иначе, формулу средней точки
(8)
Учитывая априорно большую точность формулы (8) по сравнению с формулами (6) и (7) при том же объёме и характере вычислений, эту симметричную формулу будем впредь называть просто формулой прямоугольников.
Полученные из определения интеграла квадратурные правила (6)—(8) не содержат в себе сведений, позволяющих сказать, каким нужно взять число п элементарных промежутков [xi-1, xi] или, иначе, каким должен быть шаг h разбиения (4) отрезка интегрирования [а, b], чтобы гарантированно найти значение I интеграла с наперед заданной точностью ε. Ответ на этот вопрос дает формула остаточного члена (глобальной погрешности) квадратурной формулы прямоугольников:
(9)
Как видно из формулы (9), при увеличении числа п элементарных отрезков, на которые разбивается промежуток интегрирования [а, b), ошибка численного интегрирования по формуле средней точки (8) убывает пропорционально квадрату шага h. Погрешность численного интегрирования непрерывно дифференцируемой функции по формулам левых и правых прямоугольников (6), (.7) убывает лишь по линейному закону.
Семейство квадратурных формул Ньютона-Котеса
Подстановка в интеграл вместо функции f(x) её интерполяционного многочлена Лагранжа той или иной степени п приводит к семейству квадратурных формул, называемых формулами Ньютона-Котеса.
Функция f(x) может быть единственным образом представлена в виде

где — интерполяционный многочлен Лагранжа, — остаточный член, Если система узлов интерполирования совпадает с точками разбиения (4) отрезка [а, b] с шагом h, то замена переменной x=x0+qh трансформирует многочлен Лагранжа к виду
(10)
Для того, чтобы использовать такое выражение Ln(x) вместо f(х) в нужно изменить границы интегрирования (значению х=а соответствует значение q= 0,ах=b — значение q = п) и учесть, что dx=hdq. Таким образом, получаем

Это равенство, переписанное в виде
(11)
и есть квадратурная формула Ньютона-Котеса, где
(12)
— коэффициенты Котеса.
На самом деле, формулы (11)-(12) определяют семейство квадратурных формул. Параметром этого семейства является число п — степень интерполяционного многочлена, которым заменяется подынтегральная функция.
Для практического применения целесообразно рассмотреть некоторые частные случаи.
Пусть n = 1. Это значит, что диапазон интегрирования [a, b] содержит только два узла: x0 = a и x1 = b.
Вычислим коэффициенты Котеса для этого случая. Подставляя в (12) n = 1 и i = 0, 1 получим:
(13)
Подстановка коэффициентов (13) в (12) дает квадратурную формулу Ньютона-Котеса первого порядка:
(14)
В данном случае величина шага h равна длине всего диапазона интегрирования [a, b].
При таком приближении подынтегральная функция f (x) заменяется линейной (полиномом Лагранжа 1-й степени). Графически значение определенного интеграла I(f, a, b) изображается площадью криволинейной трапеции, ограниченной сверху (или снизу) графиком функции f (x) (см. пример на рис. 3).

Рис.3. К выводу формулы трапеций.
Сплошная линия - график интегрируемой функции f (x),
штрихпунктирная - график линейного приближения интегрируемой функции на интервале [x0, х1].
Крестиками заполнена область R, площадь которой равна погрешности квадратурной формулы (14).
Приближенное значение интеграла, выраженное квадратурной суммой (14), равно площади обычной трапеции с той же высотой h. Площадь криволинейного сектора дает погрешность R вычисления определенного интеграла с помощью квадратурной формулы (14).
Когда данный способ вычисления определенного интеграла используют на практике, то обычно диапазон интегрирования [a, b] разбивают на N одинаковых поддиапазонов шириной
h = (b - a) /N (15)
и к каждому из них применяют формулу (14). Тогда искомый интеграл представится суммой, каждый член которой будет содержать общий множитель (15). Каждое значение интегрируемой функции yk (k = 0, 1, …, N) в эту сумму будет входить дважды, кроме двух крайних y0 и yN. В результате интеграл I(f, a, b) выразится следующей суммой:
(16)
Последняя формула называется формулой трапеций.
Теперь пусть n=2. Это значит, что диапазон интегрирования [a, b] содержит три узла: х0=a, x1=(a+b)/2, x2=b. В данном случае шаг интегрирования равен половине диапазона интегрирования h=(b-a)/ 2
Коэффициенты Котеса в данном случае получаются вычислением интегралов (11) с подстановкой n = 2 для i = 0, 1, 2:
(7)
Искомый интеграл представится в виде:
(18)где (b - a) = 2h.
Формула (18) была получена заменой подынтегральной функции f(x) на квадратичную, проходящую через заданные точки (xi, yi) , i = 0, 1, 2. Полученное значение искомого интеграла графически изображается площадью под параболой (см. рис. 4).

Рис. 4. К выводу формулы Симпсона.
Жирная сплошная линия - график интегрируемой функции f (x),
штриховая - график квадратичного приближения L2(x) интегрируемой функции на интервале [x0, x2]. Заштрихованная область R выражает погрешность квадратурной формулы (18).
Видно, что на интервале [x0, x1] квадратичное приближение дает заниженный результат, на интервале [x1, x2] - завышенный.
На практике весь диапазон интегрирования [a, b] обычно предварительно разбивается на четное число N = 2M равных интервалов, которые попарно объединяются в M поддиапазонов. Каждый поддиапазон содержит три узла: два крайних и один внутренний. К каждому из M поддиапазонов применяется формула (18), т.е. каждой тройке узлов соответствует своя интерполирующая парабола. Параболы стыкуются в совпадающих крайних узлах отдельных поддиапазонов. Номера всех таких узлов четные от i = 2 до i = 2 M - 2, поэтому при суммировании выражений (18) по поддиапазонам следует учесть, что слагаемые с этими номерами встретятся дважды.
В результате суммирования с учетом значений коэффициентов Котеса (17) получается рабочая формула
(19)
где для краткости введены обозначения следующих сумм:
(20)
Заметим, что сумма So содержит M слагаемых, а сумма Se состоит из M-1 слагаемых. Шаг интегрирования равен h=(b - a)/(2M)=(b - a)/N.
Формула (19) называется квадратурной формулой Симпсона.
При n = 3 аналогичным путем получается квадратурная формула с тремя узлами на интервале интегрирования:
(21)
где h = (b - a) / 3.
При решении прикладных задач диапазон интегрирования [a, b] предварительно разбивается на равные поддиапазоны, число которых кратно трем (N = 3M), и к каждому из поддиапазонов применяется формула (20). Тогда получится рабочая формула вида:
(22)
Шаг интегрирования равен h=(b-a)/(3M) = (b-a)/N, а отдельные суммы выражены следующими формулами:
(23)
Сумма (22) называется квадратурной формулой Ньютона или формулой «трех восьмых». В практике квадратурная формула Ньютона употребляется гораздо реже, чем формулы Симпсона и трапеций.
Задания: выполнить задание 3 ИДЗ№3.
Занятие № 15. Квадратурные формулы Гаусса.Все формулы численного интегрирования дают результат с некоторой погрешностью. Здесь заметим, что уменьшение погрешности расчета является одной из важнейших задач численных методов.
В предыдущем параграфе рассматривались различные квадратурные формулы для равномерного расположения узлов в интервале интегрирования [a, b]. Оказывается, при фиксированном количестве узлов можно добиться значительного уменьшения погрешности вычисления определенного интеграла, если отказаться от их равномерного расположения.
Для нахождения оптимального расположения узлов (которое обеспечивает уменьшение погрешности численного интегрирования) прежде всего сведем интервал интегрирования общего вида [a, b] к стандартному интервалу [-1, 1] с помощью линейной замены переменной интегрирования:
x = (а + b) / 2 + t (b - а) / 2. (25)
Исходный интеграл приобретет вид:
(26)
где t - новая переменная интегрирования.
В соответствии с общей квадратурной формулой запишем:
(27)
Нашей целью является выбор таких п узлов ti внутри интервала интегрирования [-1, 1] и таких п коэффициентов Ai (i = 1, п), чтобы формула (27) была точной для всех полиномов максимально большой степени. Мы имеем возможность варьировать 2n величин ti и Ai (i=1, …, п), следовательно, эта максимальная степень полинома равна N = 2n - 1.
Далее требуется доказать следующую лемму.
Лемма. Для достижения поставленной цели необходимо и достаточно, чтобы формула (27) выполнялась точно для следующих 2n подынтегральных функций
Доказательство.
По условиям леммы выполняются следующие равенства для степенных функций f(t) = tk (k = 0, 1, …, 2п- 1):
(28)
Как было принято выше, подынтегральную функцию f(t) представим в виде полинома степени 2п - 1:
(29)
где Ck - коэффициенты полинома.
Проведем интегрирование полинома (29) по интервалу [-1, 1], используя условие (28) и изменяя порядок суммирования:
(30)
Последнее равенство обусловлено определением (29).
Мы получили, что равенство (27) является точным для полинома вида (29), при этом в ходе преобразований мы воспользовались условиями (28).
Лемма доказана.
Эта лемма обеспечивает справедливость следующего важного утверждения. Пусть нам удастся найти такие числа ti и Ai (i = 1, …, n), обращающие формулу (27) в точную, используя в качестве подынтегральных степенные функции вида tk (k = 0, 1, ..., 2n - 1). Тогда эти найденные числа ti и Ai обеспечат точность формулы (27) и для любой подынтегральной функции, которая представляется в виде полинома степени 2n - 1.
Теперь в уравнениях (28) проинтегрируем аналитически левые части:
(31)
Сравнение (31) с правыми частями уравнений (28) дает нам систему 2n уравнений для 2n искомых неизвестных ti , Ai (i = 1, ..., n):
(32)
Система является нелинейной относительно неизвестных ti, Ai (i = 1,…, n). Для ее решения целесообразно воспользоваться свойствами полиномов Лежандра. Краткие сведения о полиномах Лежандра приведены в приложении 1.
Выберем подынтегральные функции f(x) в виде:
(33)
где Pn (t) - полином Лежандра степени n.
Это значит, что все n функций (33) являются полиномами степени ≤ 2n-1. Для таких подынтегральных функций, согласно доказанному выше, формула (27) является точной и коэффициенты Ai (i = 1,…, n) удовлетворяют системе (32).
Применим квадратурную формулу (27) для функций (33):
(34)
В силу свойства ортогональности полиномов Лежандра левые части всех уравнений (34) равны нулю. Следовательно, и правые части (34) будут равны нулю
(35)
при произвольных значениях Ai, если положить Pn (ti) = 0.
Таким образом, оказывается, что в качестве узлов ti следует взять точки нулей полиномов Лежандра степени n. Точки нулей полиномов Лежандра различных степеней рассчитаны с большой точностью и сведены в специальные таблицы. Характерным свойством узлов ti и коэффициентов Ai в данном методе является симметрия их значений относительно центра таблицы t = 0.
Если числа ti известны, то система (32) становится линейной относительно величин Ai. Для ее решения можно применить любой из описанных методов. Численные значения узлов ti и коэффициентов Ai с точностью до восьми знаков для различных n приведены в таблице приложения 2.Вычислив значения узлов ti и коэффициентов Ai для выбранного n, можно вернуться к исходной переменной интегрирования x и записать квадратурную формулу Гаусса в следующем виде:
(36)
где
(37)
Точность интегрирования методом Гаусса резко возрастает с увеличением числа узлов в интервале интегрирования [a, b]. Теоретические оценки и расчеты показывают, что для интервалов небольшой величины обычно достаточная точность обеспечивается при n = 8. Соответствующая таблица узлов ti и коэффициентов Ai приведена в приложении 2.
Если диапазон интегрирования [a, b] имеет значительную ширину, его можно разбить на несколько равных поддиапазонов и провести процедуру вычисления интеграла методом Гаусса для каждого поддиапазона отдельно. Интеграл по j-му поддиапазону с границами [aj , bj ] с помощью линейного преобразования (37) представляется в виде суммы, аналогичной (36):
(38)
где
Искомое значение определенного интеграла по полному диапазону [a, b] вычисляется как сумме интегралов (38) по всем поддиапазонам.
Задания: выполнить задание 3 ИДЗ№3.
Занятие № 16. Численные методы решения задачи Коши для обыкновенных дифференциальных уравнений. Классификация методов. Метод Пикара.Будем рассматривать обыкновенное дифференциальное уравнение (ОДУ) первого порядка
(1)
с начальным условием
y(х0) = у0, (2)
где f(x) — некоторая заданная, в общем случае, нелинейная функция двух переменных. Будем считать, что для данной задачи (1)-(2), называемой начальной задачей или задачей Коши, выполняются требования, обеспечивающие существование и единственность на отрезке [х0,b] ее решения у=у(х).
Несмотря на внешнюю простоту уравнения (1), решить его аналитически, т.е. найти общее решение у=у(х, С) с тем, чтобы затем выделить из него интегральную кривую у=у(х), проходящую через заданную точку (х0;у0), удается лишь для некоторых специальных типов таких уравнений. Поэтому, как и в родственной для (1)-(2) задаче вычисления интегралов, приходится делать ставку на приближенные способы решения начальных задач для ОДУ, которые можно разделить на три группы:
приближенно-аналитические методы;
графические или машинно-графические методы;
численные методы.
К методам первой группы относят такие, которые позволяют находить приближение решения у(х) сразу в виде некоторой «хорошей» функции φ(х). Например, широко известен метод степенных рядов, в одну из реализаций которого заложено представление искомой функции у(х) отрезком ряда Тейлора, где тейлоровские коэффициенты, содержащие производные высших порядков, находят последовательным дифференцированием самого уравнения (1). Другим представителем этой группы методов является метод последовательных приближений, суть которого приведена чуть ниже.
Название графические методы говорит о приближенном представлении искомого решения у(х) на промежутке [x0,b] в виде графика, который можно строить по тем или иным правилам, связанным с графическим толкованием данной задачи. Физическая или, возможно, точнее будет сказать, электротехническая интерпретация начальных задач для определенных видов уравнений лежит в основе машинно-графических методов приближенного решения. Реализуя на физико-техническом уровне заданные электрические процессы, на экране осциллографа наблюдают поведение решений дифференциальных уравнений, описывающих эти процессы. Изменение параметров уравнения приводит к адекватному изменению поведения решений, что положено в основу специализированных аналоговых вычислительных машин (АВМ).
Наконец, наиболее значимыми в настоящее время, характеризуемое бурным развитием и проникновением во все сферы человеческой деятельности цифровой вычислительной техники, являются численные методы решения дифференциальных уравнений, предполагающие получение числовой таблицы приближенных значений yi искомого решения у(х) на некоторой сетке значений аргумента х. Этим способам и будет посвящено дальнейшее изложение. Что делать с получаемыми численными значениями решения, зависит от прикладной постановки задачи. Если речь идет о нахождении только значения у(b), тогда точка b включается как конечная в систему расчетных точек хi, и все приближенные значения yi ≈y(xi), кроме последнего, участвуют лишь как промежуточные, т.е. не требуют ни запоминания, ни обработки. Если же нужно иметь приближенное решение у(х) в любой точке х, то для этого к получаемой числовой таблице значений yi можно применить какой-либо из способов аппроксимации табличных функций, рассмотренных ранее, например, интерполяцию или сплайн-интерполяцию. Возможны и другие использования численных данных о решении.
Коснемся одного приближенно-аналитического способа решения начальной задачи (1)-(2), в котором искомое решение у=у(х) в некоторой правой окрестности точки х0 является пределом последовательности получаемых определенным образом функций уп(х).
Проинтегрируем левую и правую части уравнения (1) в границах от х0 до х:

Отсюда, с учетом того, что одной из первообразных для у'(х) служит у(х), получаем

или, с использованием начального условия (2),
(3)
Таким образом, данное дифференциальное уравнение (1) с начальным условием (2) преобразовалось в интегральное уравнение (неизвестная функция здесь входит под знак интеграла).
Полученное интегральное уравнение (3) имеет вид задачи о неподвижной точке для оператора Формально к этой задаче можно применить метод простых итераций
(4)
достаточно обстоятельно рассматривавшийся применительно к системам линейных и нелинейных алгебраических и трансцендентных уравнений. Беря в качестве начальной функции y0(х) заданную в (2) постоянную y0, по формуле (4) при п=0 находим первое приближение

Его подстановка в (4) при п=1 дает второе приближение

и т.д. Таким образом, этот приближенно-аналитический метод, называемый методом последовательных приближений или методом Пикара определяется формулой
(5)
где n=0,1, 2,... и у0(х)=y0.
Отметим две характеристики метода последовательных приближений Пикара, которые можно отнести к негативным. Во-первых, в силу известных проблем с эффективным нахождением первообразных, в чистом виде метод (5) редко реализуем. Во-вторых, как видно из вышеприведенного утверждения, этот метод следует считать локальным, пригодным для приближения решения в малой правой окрестности начальной точки. Большее значение метод Пикара имеет для доказательства существования и единственности решения задачи Коши, нежели для его практического нахождения.
Занятие № 17. Методы Эйлера.Цель - ознакомить студентов с методами Эйлера решения задачи Коши для обыкновенных дифференциальных уравнений.
Метод Эйлера — разные подходы к построению
Учитывая ключевую позицию, которую занимает метод Эйлера в теории численных методов ОДУ, рассмотрим несколько способов его вывода. При этом будем считать, что вычисления проводятся с расчетным шагом расчетными точками (узлами) служат точки хi = х0 + ih (i=0,1,..., п) промежутка [x0, b] и целью является построение таблицы

приближенных значений yi решения у = у(х) задачи (1)-(2) в расчетных точках хi.
Геометрический способ. Пользуясь тем, что в точке х0 известно и значение решения у(х0) = у0 (согласно (2)), и значение его производной у'(xо)=f(x0, y0) (согласно (1)), можно записать уравнение касательной к графику искомой функции y= у(х) в точке (х0;y0):
(6)
При достаточно малом шаге h ордината
(7)
этой касательной, полученная подстановкой в правую часть (6) значения по непрерывности должна мало отучаться от ординаты у(х1) решения у(х) задачи (1)-(2). Следовательно, точка (х1, у1) пересечения касательной (6) с прямой x=x1 может быть приближенно принята за новую начальную точку. Через эту точку снова проведем прямую

которая уже приближенно отражает поведение касательной к у=у(х) в точке (х1,у(х1)). Подставляя сюда х=х2(=х1+h), иначе, пересекая эту «касательную» прямой х = х2, получим приближение значения у(х2) значением

и т.д. В итоге этого процесса, определяемого формулой
i = 0,1,..., п (8)
и называемого методом Эйлера, график решения у = у(х) данной задачи Коши (1)-(2) приближенно представляется ломаной, составленной из отрезков приближенных касательных (рис. 1), откуда происходит другое название метода (8) — метод ломаных.

Рис. 1. Геометрическая интерпретация метода Эйлера
Применение формулы Тейлора. Описываемый здесь способ вывода метода Эйлера тесно связан с предыдущим. Линеаризуя решение в окрестности начальной точки по формуле Тейлора, имеем

Отсюда при х = х1 получаем
(9)
Точное равенство (9), переписанное в виде

говорит о том, что здесь мы имеем одновременно как саму формулу Эйлера для вычисления значения (сравните с формулой (7)), так и ее остаточный член
(10)
где ξi — некоторая точка интервала (x0, x1).
Остаточный член (10) характеризует локальную (или, иначе, шаговую) ошибку метода Эйлера, т.е. ошибку, совершаемую на одном шаге. Очевидно, что от шага к шагу, т.е. при многократном применении формулы (8), возможно наложение ошибок. За п шагов, т.е. в точке b, образуется глобальная ошибка. Порядок глобальной ошибки (относительно шага h) на единицу ниже, чем порядок локальной ошибки, а порядком глобальной ошибки и определяется порядок соответствующего численного процесса решения задачи Коши. Таким образом, локальная ошибка метода Эйлера, согласно (10), есть глобальная — 0(h), т.е. метод Эйлера относится к методам первого порядка.
Разностный способ. Рассматривая уравнение (1) в точке x0 с учетом (2) имеем равенство

Применяя к его левой части аппроксимацию производной правым разностным отношением первого порядка точности

получаем

что идентично равенству (9), поставляющему формулу для вычисления у1 вида (7) и локальный остаточный член (10). Ясно, что для получения общей расчетной формулы (8) можно было сразу применить аппроксимацию по формуле (6.16) в равенстве
(11)
заменив неизвестное точное значение y(xi) известным приближенным значением yi.
Заметим, что порядок получающегося таким способом метода численного интегрирования дифференциальной задачи (1)-(2) совпадает с порядком аппроксимации производной в левой части уравнения (1).
Квадратурный способ. Как было показано начальную задачу для ОДУ (1)-(2) можно заменить эквивалентным интегральным уравнением (3). При х=х1 из него получится равенство
(12)
Применение к интегралу в правой части равенства (12) простейшей (одноточечной) формулы левых прямоугольников дает приближенную формулу

правая часть которой, очевидно, совпадает с выражением (7) для подсчета значения у1. В общем случае расчетная формула (8) метода Эйлера получается численным интегрированием посредством простейшей формулы левых прямоугольников в равенстве
(13)
в предположении, что на каждом i-м шаге в роли начальной точки (x0,y0) выступает точка (xi, yi). Зная точность используемой в (13) квадратурной формулы, легко прийти к тому же выражению локальной ошибки метода Эйлера, что и при других способах его построения.
Существуют и другие подходы к выводу метода Эйлера. В частности, он будет возникать далее как частный случай некоторых семейств численных методов решения задачи (1)-(2).
Несколько простых модификаций метода Эйлера
Разовьем последний из подходов к построению метода Эйлера. Очевидно, применение к интегральному равенству (13) других простейших квадратурных формул будет порождать новые методы численного интегрирования задачи Коши (1)-(2).
Так, если в (13) использовать простейшую квадратурную формулу правых прямоугольников, придем к методу
(14)
Этот метод называют неявным (или обратным) методом Эйлера, поскольку для вычисления неизвестного значения по известному значению требуется решать уравнение, в общем случае нелинейное.
Применение к интегралу в (13) простейшей квадратурной формулы трапеций приводит тоже к неявному методу
(15)
который будем называть методом трапеций. Квадратурная формула трапеций на порядок точнее формул левых и правых прямоугольников. Отсюда вытекает более высокий (на единицу) порядок точности метода трапеций (15) по сравнению с явным и с неявным методами Эйлера (8) и (14), т.е. метод трапеций (15) — это метод второго порядка.
Некоторый интерес представляет совместное применение явного метода Эйлера и неявного метода трапеций.
По форме равенство (15) представляет собой скалярную задачу о неподвижной точке относительно неизвестного значения уi+1. Поэтому, если в правую часть (15) подставить хорошее начальное приближение , подсчитываемое по формуле (14), то тогда само это равенство можно считать шагом метода простых итераций для уточнения этого значения. Таким образом, получаем гибридный метод
i=0,1,..., n, (16)
который называют методом Хойна.
Ясно, что можно достичь большей точности, если, исходя из того же начального приближения

сделать не одну, а несколько итераций по методу трапеций:
(17)
Такой вариант совместного применения метода Эйлера и метода трапеций называют усовершенствованным методом Эйлера-Коши с итерационной обработкой. Делать много итераций по формуле (17) не рекомендуется (обычно их выполняют не более трех-четырех). Совпадение определенного числа разрядов в полученных числах говорит о точности, с которой решено методом простых итераций уравнение (15) относительно уi+1, а вовсе не о том, что с такой точностью найдено значение у(хi+1).
Чтобы получить следующую модификацию метода Эйлера,проинтегрируем уравнение (1) по отрезку Имеем

откуда следует равенство
(18)
Применяя к последнему интегралу одноточечную квадратурную формулу средних прямоугольников и заменяя значения y(xi-1) и y(xi) известными приближенными значениями уi-1 и уi соответственно, из (18) выводим формулу для подсчета приближенного значения y(xi+1)
(19)
которую будем называть уточненный методом Эйлера.
Как известно, квадратурная формула прямоугольников (средней точки) имеет тот же порядок точности, что и квадратурная формула трапеций, так что уточненный метод Эйлера (19) тоже является методом второго порядка. Подтверждением этого факта может служить вывод метода (19) на разностной основе. Применив к равенству (11) формулу симметричной аппроксимации y'(xi) второго порядка точности, получим

откуда после приближенной замены следует(19).
Обратим внимание на одно принципиальное отличие метода (19) от всех других рассмотренных до этого момента методов: метод (19) является двухшаговым. Здесь для вычисления значения уi+1 привлекаются два предыдущих значения уi и уi-1. Двухшаговость накладывает определенные ограничения, по крайней мере, на начало численного процесса: значение не может быть найдено непосредственно этим методом с тем же шагом h. Поэтому недостающую вторую начальную для процесса (19) точку приходится получать другим путем, например, явным методом Эйлера, а чтобы не сделать сразу большой ошибки, применяя на старте метод более низкого порядка точности, рекомендуется осуществлять постепенное вхождение в процесс (19). Так, «разгон» можно выполнить по формулам
(20)
а далее уже переключаться на счет по формуле (19).
Пример 1. Рассмотрим простое линейное уравнение с начальным условием y(0)=1.
На этой задаче легко проследить за вычислениями, реализующими различные выведенные выше методы. Знание ее точного решения позволяет провести сравнение результатов приближенных вычислений по разным формулам с истинным решением и проверить, насколько соответствуют представления о точности тех или иных методов тому, что наблюдается в данном, наверное, далеко не самом типичном частном случае.
Сначала сделаем несколько последовательных приближений по методу Пикара. Его итерационная формула (5) для данной начальной задачи имеет вид

Подставляя сюда у0 =1 при п = 0,1, 2 последовательно получаем:

Эти результаты удобно сравнить с точным решением, если в последнем разложить в ряд по степеням х фигурирующую там функцию e-3x. Тогда получим представление решения в виде ряда

с которым, как видим, хорошо согласуются приближения y1, y2, y3, определяемые методом Пикара.
Теперь проведем подсчет приближенных значений решения у(х) данной задачи в точке х=0.2 численным методом Эйлера и его модификациями, принимая h = 0.1 (т.е. за два шага). Результаты этих вычислении и фактические ошибки, найденные сравнением с точным значением y(0.2) = 0.581881..., отражены в следующей таблице.
Метод
Эйлера (8) 0.7 0.51 ≈0.07
Неявный Эйлера (14) ≈0.7846 ≈0.6343 ≈-0.05
Трапеций (15) ≈0.7478 ≈0.5788 ≈0.003
Хойна(16) 0.755 ≈0.5895 ≈-0.008
Уточненный Эйлера (19)-(20) 0.755 0.587 ≈-0.005
Исправленный метод Эйлера
Пусть найдено приближенное значение решения у = у(х) задачи (1)-(2) и требуется вычислить где Запишем разложение решения по формуле Тейлора р-го порядка, принимая за базовую точку xi (т.е. по степеням х-хi) и положим в этом разложении x= xi+1. Имеем
(21)
Если ограничится двумя слагаемыми в правой части разложения (21), то, получим обычный метод Эйлера (8). Посмотрим, что дает учитывание третьего слагаемого.
При р = 2 из (21) следует равенство
(22)
Значение первой производной в точке xi в силу связи (1), приближенно известно:
(23)
Дифференцируя (1), по формуле полной производной

находим приближенное значение второй производной:
(24)
Подставляя приближенные выражения ) в равенство (22), получаем следующую формулу для вычисления при i=0, 1, …, n: (25)
Определяемый ею метод будем называть исправленным методом Эйлера.
Так как при i=0 формулы (23) и (24) точны, а y0=y(x0),согласно начальному условию (2), то на первом шаге вычислений по формуле (25) будет совершаться ошибка, связанная только с усечением ряда Тейлора. Следовательно, локальная ошибка или, иначе, шаговая погрешность метода (25) составляет величину o(h3), а это означает, что исправленный метод Эйлера относится к методам второго порядка.
Занятие № 18. Методы Рунге-Кутта.О семействе методов Рунге-Кутта. Методы второго порядка
Недостатком исправленного метода Эйлера (25) и других методов более высоких порядков, основанных на пошаговом представлении решения у(х) задачи (1)-(2) по формуле Тейлора и последовательном дифференцировании уравнения (1) для получения тейлоровых коэффициентов, является необходимость вычисления на каждом шаге частных производных функции f(x,у).
Идея построения явных методов Рунге-Кутты p-го порядка заключается в получении приближений к значениям f(хi+1) по формуле вида
(26)
где φ(х, у, h) — некоторая функция, приближающая отрезок ряда Тейлора до р-го порядка и не содержащая частных производных функции f(х, у).
Так, полагая в (26) φ(х, у, р) = f(x, у) приходим к методу Эйлера (8), т.е. метод Эйлера можно считать простейшим примером методов Рунге-Кутты, соответствующим случаю р = 1.
Для построения методов Рунге-Кутты порядка, выше первого, функцию φ(х, у, h) берут многопараметрической, и подбирают ее параметры сравнением выражения (26) с многочленом Тейлора для у(х) соответствующей желаемому порядку степени.
Рассмотрим случай р=2. Возьмем функцию φ в (26) следующей структуры:

Ее параметры с1, с2, а и b будем подбирать так, чтобы записанная, согласно (26), формула
(27)
определяла метод второго порядка, т.е. чтобы максимальная локальная ошибка составляла величину o(h3).
Разложим функцию двух переменных f(x + ah, у + bhf(x, у)) по формуле Тейлора, ограничиваясь линейными членами:

Ее подстановка в (27) дает
(28)
Сравнение последнего выражения с тейлоровским квадратичным представлением решения у(х) (22) с точностью до о(h3) равносильно сравнению его с выражением уi+1 по формуле (25), т е. с исправленным методом Эйлера. Очевидно, чтобы (28) и (25) совпадали с точностью O(h3), от параметров нужно потребовать выполнение следующей совокупности условий:
(29)
Полученная система условий содержит три уравнения относительно четырех параметров метода. Это говорит о наличии одного свободного параметра. Положим с2 = α(≠ 0). Тогда из (29) имеем:

В результате подстановки этих значений параметров в формулу (27) приходим к однопараметрическому семейству методов Рунге-Кутты второго порядка.
(30)
Выделим из семейства методов (30) два наиболее простых и естественных частных случая:
при α=1/2 получаем формулу

в котором узнаём метод Хойна (16), полученный ранее из других соображений;
при α = 1 из (30) выводим новый простой метод
(31)
который назовем методом средней точки.
Методы Рунге-Кутты произвольного и четвертого порядков
Любой метод из семейства методов Рунге-Кутты второго порядка (30) реализуют по следующей схеме. На каждом шаге, т.е. при каждом i=0, 1, 2,..., вычисляют значения функции

а затем находят шаговую поправку

прибавление которой к результату предыдущего шага дает приближенное значение решения у(х) в точке xi+1 = xi + h:

Метод такой структуры называют двухэтапным по количеству вычислений значений функции — правой части уравнения (1) — на одном шаге.
Анализ устройства методов Рунге-Кутты второго порядка позволяет представить, в какой форме следует конструировать явный метод Рунге-Кутты произвольного порядка. По аналогии с предыдущим для семейства методов Рунге-Кутты p-го порядка используется запись, состоящая из следующей совокупности формул:
(32)
где к = 2, 3, …, р (для р-этапного метода). Многочисленные параметры сk, ак, bkj, фигурирующие в формулах (32), подбираются так, чтобы получаемое методом (32) значение уi+1 совпадало со значением разложения y(xi+1) по формуле Тейлора с погрешностью O(hp+1) (без учета погрешностей, совершаемых на предыдущих шагах).
Наиболее употребительным частным случаем семейства методов (32) является следующий метод Рунге-Кутты четвертого порядка относящийся к четырехэтапным и имеющий вид:
(33)
Не пытаясь воспроизвести выкладки, приводящие от общей записи семейства (32) при р=4 к конкретному методу (33), дадим геометрическое толкование последнего.
Обратив внимание на то, что шаговая поправка Δуi, есть средневзвешенная величина поправок каждого этапа (с весовыми коэффициентами1/6, 2/6, 2/6, 1/6 соответственно), проанализируем, как получаются эти поправки этапов. На первом этапе создается приращение соответствующее шаговой поправке Эйлера, — этоочевидно. На рис 2 ему отвечает отрезок ВС вертикалих=xi+1 (точка В получена ортогональным проектированиемточки А на эту вертикаль).

Рис. 2. Геометрическая иллюстрация одного шага методов Рунге-Кутта четвертого порядка
Так как точка М, благодаря свойству средней линии треугольника (см. ΔАВС), имеет ординату определяет значение f(M), служащее (согласно связи у=f(x, у) и геометрическому смыслу производной) тангенсом угла А в новом треугольнике с противолежащим этому углу катетом Далее, аналогично, подсчитав на вертикали x=xi+1 откладываем следующую промежуточную (этапную) поправку Вычислив величину f(E)= являющуюся значением тангенса угла А во вновь получаемом ΔABG, имеем поправку последнего этапа. Итоговая шаговая поправка есть продукт усреднения с указанными коэффициентами четырех этапных поправок — длин отрезков ВС, BD, BE и BG. Точка Н будет стартовой для следующего, i+1-го, шага метода (33).
Заметим, что если первый этап, как уже упоминалось, соответствует применению явного метода Эйлера, то четвертый — неявного, а второй и третий — уточненного методов Эйлера. Последний имеет более высокий порядок точности, отсюда и больший вес отвечающих ему значений этапных поправок.
Занятие № 19. Пошаговый контроль точности. Метод Кутты-Мерсона.Цель - ознакомить студентов с пошаговым контролем точности и методом Кутты-Мерсона решения задачи Коши для обыкновенных дифференциальных уравнений.
Нетрудно понять, что выведение надежных и, в то же время, простых и эффективных оценок погрешности, гарантирующих получение таблицы значений решения у=у(х) заданной точности, является делом малоперспективным, особенно для методов более-менее высоких порядков. Поэтому главным способом отслеживания точности при реализации численных процесcов решения задачи Коши остается применение различных полуэмпирических правил, основанных на принципе Рунге.
Будем считать, что при использовании метода р-го порядка абсолютная шаговая погрешность должна находиться в пределах ε>0. Тогда, согласно принципу Рунге, осуществляется счет по системе узлов с шагом h и по системе узлов — с шагом h/2. При четных j вторая система будет совпадать с первой, т.е. Переход от расчетной точки xi с приближенным значением решения в ней уi к расчетной точке один раз совершается за один шаг длины h и приводит к значению другой раз — за два шага длины h/2 («транзитом» через точку со значением и дает значение
Поправка Ричардсона в таком случае будет составлять величину
(34)
Если величина меньше заданного ε, то можно считать, что ошибка приближенного равенства не превосходит ε. Если же то следует уменьшить расчетный шаг h. При условии стоит попытаться двигаться дальше с более крупным шагом (например, удвоить h).
Пример 2. (продолжение примера 1). Посмотрим, что дает применение принципа Рунге к нескольку простым методам численного решения того же уравнения y'=2х-3y начальным условием у(0) = 1. Из точки х=0 перейдем в точку х=0.2 за один шаг h =0.2 четырьмя одношаговыми методами: явным и неявным методами Эйлера, методом трапеций и методом Хойна — частным случаев метода Рунге-Кутты второго порядка. С помощью полученных значений и найденных ранее теми же методами в примере 1 значений подсчитаем поправки Ричардсона

при р=1 для методов Эйлера и р=2 для методов трапеций и Хойна. Эти результаты, а также уточненные прибавлением к значениям поправок Ричардсона приближенные значения решения y(0.2) и их истинные погрешности сведем в следующую таблицу.
Метод
Эйлера (8) 0.4 0.11 0.62 ≈-0.04 ≈0.07
Неявный Эйлера (14) 0.675 ≈ -0.0407 ≈ 0.5936 ≈-0.01 ≈ -0.05
Трапеций (15) ≈0.5692 ≈0.0032 ≈0.5820 ≈-0.0001 ≈0.003
Хойна (16) 0.62 ≈ -0.0102 ≈0.5793 ≈ 0.0026 ≈ -0.008
В эту таблицу последним столбцом помещен последний столбец из таблицы результатов примера 1, содержащий погрешности значений Сравнение с ним столбца со значениями поправок Ричардсона показывает, что эти поправки хорошо отражают поведение погрешностей методов (хотя и не дают основания считать их модули оценками погрешностей), а предпоследнего — эффективность уточнения по правилу Рунге-Ричардсона.
Грубо обозначенная здесь технология пошагового контроля точности численного интегрирования дифференциальных уравнений и автоматического выбора расчетного шага при этом на основе двойного счета в такой непосредственной форме говорит о ее значительной «дороговизне». Действительно, предположим, что для решения задачи (1)-(2) применяется четырехэтапный метод Рунге-Кутты четвертого порядка (33). Тогда выполнение одного его шага с контролем точности по правилу Рунге потребует 11 вычислений правой части уравнения (1) (по четыре получения каждого из значений и минус одно общее для весьма затратно.Более «дешевый», но, возможно, менее строгий способ судить о том, достаточно ли малым выбран шаг h расчетов по методу Рунге-Кутты четвертого порядка (33), — это вычисление при каждом i= 0,1, 2,... величин

Считается, что если величина Θi, не превосходит нескольких сотых, то можно продолжить вычисления с данным шагом или пытаться при переходе от i к i+1 его увеличить; в противном случае шаг следует уменьшить, например, вдвое.
Стремление повысить вычислительную эффективность привело к появлению различных вычислительных версий методов Рунге-Кутты, благо для этого в семействе методов (32) имеется значительное число свободных параметров. Основные соображения, положенные в основу этих версий, таковы: нужно получить формулы из семейства методов Рунге-Кутты (32), которые использовали бы одни и те же значения функции — правой части уравнения (1) — и определяли бы разные конкретные методы одного порядка (или смежных порядков, например, четвертого и пятого); при этом, чтобы по разности результатов подсчета приближенных значений решения по выведенным близким формулам (с одним и тем же шагом h(!)) можно было судить о точности одного из них.
Приведем один из таких методов, который называется методом Кутты-Мерсона или, иначе, пятиэтапным методом Рунге-Кутты четвертого порядка, а также методов вложенных форм.
На i-ом шаге решения задачи (1)-(2) последовательно вычисляют:

После этого подсчитывают величину

и проводят сравнения. Если значение R окажется больше заданного допустимого уровня абсолютных погрешностей ε, то шаг уменьшают вдвое и возвращаются к началу второго этапа, т.е. заново вычисляют и т.д. Если R≤ε, то считают с точностью ε. При переходе к следующему шагу делается проверка на возможность увеличить расчетный шаг: если то далее расчет ведется с шагом h:= 2h.
Занятие № 20. Многошаговые методы Адамса.Цель - ознакомить студентов с многошаговыми методами Адамса решения задачи Коши для обыкновенных дифференциальных уравнений.
Будем строить численные методы решения начальной задачи
(1)
(2)
Будем считать, что уже найдено несколько приближенных значений (j = 0,1,..., i) решения у=у(х) задачи (1)-(2) на равномерной сетке xj=x0+jh, и нужно получить правило для вычисления очередного значения Для вывода таких правил используем интегро-интерполяционный подход. А именно, проинтегрировав левую и правую части уравнения (1) по промежутку[xi , xi+1], в полученном равенстве
(3)под интеграл вместо функции f(x,у(х)) подставим интерполирующий ее многочлен Хотя выражение функции f(x,y(x))как функции одной переменной х, вообще говоря, неизвестно, ее дискретные приближенные значения обозначаемые в дальнейшем для краткости fj, при j=1, 2, ..., j можно считать известными. В таком случае, дополняя эти известные значения пока что неизвестным значением можно построить таблицу конечных разностей (табл. 1), служащую основой для образованияинтерполяционных многочленов к-й степени для интерполирования назад из точек
Таблица 1
Таблица конечных разностей для построения конечноразностных формул Адамса

При интерполировании назад из узла хi, по второй интерполяционной формуле Ньютона имеем
(4)
(см.конечные разности, подчеркнутые в табл. 1 сплошной линией), а из узла xi+1 по той же формуле получаем многочлен
(5)
(использующий разности, подчеркнутые пунктиром).
Подстановка многочленов в равенство (3) приводит к формулам для вычисления очередного значения вида
(6)
и
(7)
В результате применения к интегралам в (6) и (7) формулы Ньютона-Лейбница получается два семейства методов (с параметром k∈N0), называемых многошаговыми методами Адамса. Рассмотрим по отдельности каждое из этих семейств.
Экстраполяционные методы Адамса-Башфорта. Чтобы подставить в (6) многочлен (4), зависящий от переменной сделаем в интеграле замену переменной в соответствии с которой
Тогда формула (6) может быть переписана в виде
(8)
где

(9)
Таким образом, на основе (8) получается следующая конечно-разностная формула, определяющая экстраполяционный метод Адамса-Башфорта:
(10)
Посмотрим, что представляют собой наиболее простые частные случаи метода Адамса-Башфорта, соответствующие нескольким первым значениям параметра к в формуле (8). Сразу заметим, что при фиксировании к = 0,1, 2,... в (8) тем самым задается степень интерполяционного многочлена (нулевая, первая, вторая и т.д.) и, соответственно, число слагаемых, равное 1, 2, 3, в правой части (9) (или, что то же, в скобках формулы (10)). Конечные разности в получающихся при этом конкретных формулах будем раскрывать через значения функции, приводя формулы к виду, называемому иногда ординатным. Имеем:
при к = 0 (11)
при к = 1
(12)
при к =2
(13)
при k= 3

(14)
Формулы (11), (12), (13) и (14) определяют методы Адамса-Башфорта соответственно первого, второго, третьего и четвертого порядков. Относительно порядка метода (11) сомнений нет: мы узнаём метод Эйлера (8).
Интерполяционные методы Адамса-Моултона. В интеграле, фигурирующем в формуле (7), делаем заменуx=xi+1+qh и подставляем в него выражениеопределяемое формулой (5). Приходим к аналогичному (8) равенству

где

(15)
Отсюда следует конечноразностная формула интерполяционного метода Адамса-Моултона
(16)
Аналогично тому, как это делалось для методов Адамса-Башфорта, при к= 0,1, 2,3, т.е. фиксированием одного, двух, трех, четырех членов в представлении (15) интеграла получаем следующие частные формулы:
при к=0 (17)
при к= 1 (18)
при к=2 (19)
при к=3
(20)
Формулы (17) и (18) определяют уже известные нам методы, а именно, неявный метод Эйлера (14) и метод трапеций (15), имеющие первый и второй порядки точности соответственно. Заметим, что оба эти метода являются одношаговыми, а следующие за ними методы Адамса-Моултона (19) и (20) третьего и четвертого порядков относятся, как легко видеть, соответственно к двухшаговым и трехшаговым методам. Таким образом, для интерполяционных методов Адамса-Моултона порядок шаговости на единицу ниже порядка точности метода (за тривиальным исключением, отвечающим случаю к = 0).
Важное различие в экстраполяционных и интерполяционных методах Адамса заключается в том, что первые из них являются явными, а вторые — неявными. Эти термины однозначно определяют, о каком из двух семейств методов Адамса идет речь, а их сущность диктует особенности использования методов Адамса при практических расчетах, что найдет отражение в следующем параграфе.
Занятие № 21. Методы прогноза и коррекции.Под названием методы прогноза и коррекции (иначе методы предсказания и уточнения, предиктор-корректорные методы) понимается совместное применение явных и неявных методов одинакового или смежных порядков. По явной формуле значение решения у(х) задачи (1)-(2) в текущей (расчетной) точке хi+1 прогнозируется, т.е. находится его, быть может, достаточно грубое приближение, а с помощью неявной формулы, в правую часть которой подставляется спрогнозированное значение, оно уточняется (корректируется). Пример приближенного вычисления у(хi+1) по такой явно-неявной схеме у нас уже есть: парное использование явного метода Эйлера для предсказания и метода трапеций для уточнения (17).
Остановимся подробнее на методах прогноза и коррекции, базирующихся на парах явных и неявных методов Адамса одинакового порядка. Обозначим через приближенное значение решения y(xi+1), подсчитываемое по явной экстраполяционной формуле Адамса-Башфорта, и составим несколько пар из рассмотренных в предыдущем параграфе частных формул Адамса-Башфорта (11), (12), (13), (14) и Адамса-Моултона (17), (18), (19), (20).
Имеем следующие предиктор-корректорные методы Адамса:
первого порядка (он же явно-неявный метод Эйлера)

второго порядка

третьего порядка

четвертого порядка
(21)
Одним из главных достоинств методов прогноза и коррекции является возможность контролировать шаговую погрешность сравнением двух полученных по явной и неявной формулам приближений к y(xi+1 ). Покажем, как реализуется эта возможность для наиболее употребительного предиктор-корректорного метода Адамса четвертого порядка (21).
Вспомним, что первая из формул (21) была получена из общей формулы Адамса-Башфорта (10), а вторая — из общей формулы Адамса-Моултона (16), в которых последними брались разности третьего порядка (подынтегральная функция в равенстве (3) аппроксимировалась интерполяционным многочленом третьей степени). Считая, что расчетный шаг h достаточно мал и конечные разности с ростом их порядка убывают, главные части шаговых погрешностей формул Башфорта и Моултона четвертого порядка, в соответствии с (10) и (16), характеризуются величинами для явной и для неявной формул. Таким образом, если наряду с введенным обозначением обозначить через приближенное значение у(хi+1), получаемое по формуле Адамса-Моултона четвертого порядка, то можно записать два приближенных представления y(xi+1):
(22)
и
(23)
Отсюда видно, что если четвертые разности функции f(x, у(х)) в используемой части табл. 1 конечных разностей практически постоянны ( а это можно связать с удачным выбором величины шага h при достаточном запасе знаков в значениях f(xj,yj)), то во-первых, значения и дают двусторонние приближения к точному решению y(xi+1), а во-вторых, через разность между значениями и можно оценить точность каждого из них.
Действительно, приравнивая правые части приближенных равенств (22) и (23) и отождествляя имеем:

откуда

Подставляя последнее в (23), получаем приближенное равенство
(24)
Использование приближенной формулы (24) может быть двояким. Переписав ее в виде

применяем это для пошагового контроля точности:
если то полагаем с точностью ε и переходим к следующему шагу (i:=i+1), иначе уменьшаем шаг h и снова подсчитываем
Другое назначение формулы (24) — это прямое применение ее правой части для получения уточненного значения: полагаем
(25)
Наверное, есть смысл контроль точности делать на каждом шаге, а к уточнению по формуле (25) прибегать при выводе окончательных результатов.
Замечание 1. При выводе формулы (24) под мы понимаем значение, соответствующее «чистому» методу Адамса-Моултона четвертого порядка, т.е. — это точная реализация неявной формулы (20). Вторая же формула предиктор-корректорного метода (21) соответствует лишь одному приближению к по методу простых итераций, где в качестве начального приближения берется . Поэтому применения формулы (24) к методу прогноза и коррекции (21) будут убедительны в том случае, если его вторая формула итерируется хотя бы один-два раза. Однако, чем больше таких итераций, тем ниже вычислительная эффективность этого метода, в целом весьма высокая по сравнению с многоэтапными методами Рунге-Кутты.
Занятие № 22. Метод Милна четвертого порядкаЦель - ознакомить студентов с методом Милна четвертого порядка решения задачи Коши для обыкновенных дифференциальных уравнений.
Рассмотрим еще один широко известный метод прогноза и коррекции — метод Милна.
Для вывода первой формулы Милна (т.е. формулы предсказания) проинтегрируем данное уравнение (1) на промежутке [xi-3, xi+1] и в полученном интегральном равенстве
(26)
подынтегральную функцию f(x,у(х)) заменим первым интерполяционным многочленом Ньютона Р3(х), построенным по четырем узлам с предполагающимися уже известными приближенными значениями

Тогда, после замены переменной на основании (26) имеем:

Отсюда, выразив конечные разности через значения функции, получаем первую формулу Милна (предсказания)
(27)
которую, очевидно, следует отнести к экстраполяционным.
Главный член локальной погрешности формулы (27) находим интегрированием следующего (первого из неучтенных) слагаемого интерполяционного многочлена Ньютона. Именно:

Считая четвертые разности примерно одинаковыми, опустим индекс у функции f в записи в результате получаем следующее приближенное представление решения в точке xi+1 : (28)
Вывод второй формулы Милна более прост. Проинтегрируем уравнение (1) теперь на промежутке [xi-1, xi+1] и в полученном равенстве

применим к интегралу простейшую формулу Симпсона. Имеем
(29)
Отбрасывая здесь остаточный член и заменяя значения решения y(xi-1) и y(xi;) известными приближенными значениями yi-1 и yi, а стоящее в правой части под знаком функции f неизвестное значение у(хi+1) тем значением которое получается в результате вычислений по явной первой формуле Милна (27), приходим ко второй формуле Милна (уточнения)
(30)
являющейся интерполяционной.Для вывода приближенной оценки шаговой погрешности воспользуемся приближенным равенством где так же, как и в (28), — условная запись практически постоянных четвертых разностей. Исходя из точного равенства (29), локальную погрешность получаемого с помощью формулы (30) (возможно, с итерационной обработкой, см. замечание) приближенного значения yi+1 можно приближенно охарактеризовать величиной , т.е.
(31)
Сравнение (28) и (31) дает:

следовательно,
(32)
Таким образом, при численном интегрировании начальной задачи (1)-(2) методом Милна четвертого порядка, определенным формулами (27) и (30), на каждом i-м шаге следует вычислять величину

и сравнивать ее модуль с величиной ε > 0 допустимой шаговой погрешности. Если то за у(хi+1) принимается полученное по второй формуле Милна значение уi+1 (или его уточненное значение ); иначе шаг должен быть уменьшен.
Фигурирующая в приближенном равенстве (32) постоянная 1/29 примерно вдвое меньше постоянной 19/270≈1/14 в аналогичном равенстве (24) для предиктор-корректорного метода Адамса четвертого порядка (22), что характеризует метод Милна как несколько более точный при одинаковых вычислительных затратах.
Приложение 1ПОЛИНОМЫ ЛЕЖАНДРА
Полиномы Лежандра являются специальными функциями, которые применяются при решении многих теоретических и прикладных задач. Полином Лежандра n-й степени можно определить с помощью производной n-го порядка следующим образом:
(1)
где z - комплексная переменная.
В данном учебном пособии рассматриваются и используются полиномы Лежандра для действительного аргумента x, лежащего в интервале x∈[-1, 1].
С помощью определения (1) легко получить явные выражения полиномов Лежандра действительного аргумента низших степеней:
(2)
Графики перечисленных полиномов приведены на рис.1.
Все полиномы Лежандра Pn(x) имеют следующие граничные значения:
(3)
Нетрудно убедиться, что полиномы Лежандра четной степени являются четными функциями и наоборот.
Важным для практических применений является свойство ортогональности полиномов Лежандра:
(4)
где Qk(x) - любой полином степени k, меньшей n (k < n).
Полиномы Лежандра подчиняются рекуррентному соотношению
(5)
которое, в частности, удобно для последовательного вычисления полиномы высоких степеней.


Рис.1. Графики полиномов Лежандра а) n = 0, 1, 2, б) n = 3, 4.
Приложение 2.Параметры квадратурных формул Гаусса

Индивидуальное домашнее задание № 11. В приведенных задачах числа m, n, k вычислены с некоторой погрешностью. Необходимо вычислить и определить погрешность результата для Х.
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
2. Ознакомьтесь с методами приближенного вычисления корней уравнений. Найдите один действительный корень уравнения с точностью 10-5. В ходе решения осуществить следующие шаги:
2.1. Отделить корень уравнения.
2.2. Вычислить с помощью программы значение отдельного корня методами: деление отрезка пополам, хорд, касательных, комбинированным методом, методом итераций. При использовании метода простых итераций найти решение при разных начальных приближениях. Результаты вычислений занести в таблицу.
Вариант задания выбрать из табл. 1.1.

3. Найдите действительный корень уравнения с точностью 10-4, на интервале [a,b]. На первом этапе решения методом деления пополам, уменьшать интервал, содержащий корень, до тех пор, пока его длина не станет меньше 0,2. Потом, применить один из «более» быстрых методов.
1 11
2 12
3 13
4 14
5 15
6 16
7 17
8 18
9 19
10 20
ОТЧЕТ О РАБОТЕ
Отчет должен содержать:
1. График исследуемой функции с интервалами отделения корней.
2. Таблицы пошаговых расчетов корня уравнения.
3. Обоснованное заключение о преимуществах и недостатках использования исследованных методов решения применительно к заданному уравнению (для задания 1).
4. Используя схему Гаусса (схема единственного деления и схема полного выбора) решить систему уравнений


5. Решить систему уравнений двумя способами — методом итераций и методом Зейделя. Продолжать итерации до тех пор, пока точность приближенного решения не станет меньше 0,01.


Индивидуальное задание №2.1. Функция y=f(x) задана таблицей значений:

Указания. Для вариантов 10 - 12 значения аргумента x предварительно перевести из градусов в радианы.
Даны контрольные значения аргумента x1=12; x2=26; x3=42.
Написать подходящие для приближенного вычисления значений y1=f(x1), y2=f(x2), y3=f(x3) интерполяционные многочлены Лагранжа первой и второй степени. Получить эти значения.
Составить алгоритм и написать программу на языке высокого уровня, реализующую схему Эйткена вычисления с максимально возможной точностью значения y=f(x) в произвольной точке x промежутка Пользуясь этим алгоритмом, вычислить приближенные значения
Сделать анализ результатов заданий 1, 2.
2. Для заданной таблично функции построить все возможные интерполяционные многочлены Ньютона максимальной степени, пригодные для определения значения функции в указанных промежуточных точках Для всех вариантов

3. Вычислить значения данной функции и ее прозводной с помощью интерполяционного полинома Лагранжа Ln(x). В качестве узлов интерполяции взять:
1) равномерно распределенные точки на отрезке [a; b];
2) чебышевский набор узлов на отрезке [a; b].
При табулировании функции вычислять ряд с точностью 10-6.



Замечание. При вычислении ряда учесть, что каждый последующий член ряда an+1 получается из предыдущего члена an умножением на величину qn, т.е. Это позволит избежать переполнения при вычислении факториалов.
4. Найти приближенные значения функции при данных промежуточных значениях аргумента с помощью кубического сплайна и визуализируйте результаты сплайн-интерполяции.


Отчет должен содержать:
постановку задачи и исходные данные;
описание методов решения;
графики, полученных интерполяционных многочленов;
листинг программы.
Индивидуальное задание №31. Используя данные таблицы 1, вычислить производную указанной функции в точке х (точка х не является узлом таблицы)
2. Используя данные таблицы 1, вычислить производную указанной функции в точке х (точка х – узел таблицы)
Таблица 1.
Вариант Задание 1. Задание 2.
Таблица № х Таблица № х Используемая формула
1 2 3.0 2 (взять 5 последних значений) 5,3 Лагранжа
2 3 3.5 3 (взять 4 последних значения) 6,7 Лагранжа
3 4 2.5 5 2,6 Ньютона
4 5 5.8 4 -3,2 Ньютона
5 2 3.1 3 2,3 Ньютона
6 3 3.9 2 2,1 Ньютона
7 4 3.3 4 (взять 5 первых значений) -0,8 Лагранжа
8 5 6.0 5 (взять 4 первых значения) 3,8 Лагранжа
9 2 3.2 2 (взять 5 первых значений) 2,9 Лагранжа
10 3 5.3 4 1,6 Ньютона
11 4 3.9 3 3,4 Ньютона
12 5 7.2 5 (взять 4 первых значения) 5 Лагранжа
13 2 4.4 5 6,2 Ньютона
14 3 3.6 3 (взять 5 последних значений) 4,5 Лагранжа
15 4 2.2 4 (взять 5 последних значений) 4 Лагранжа
16 5 6.8 2 3,7 Ньютона
17 2 3.4 3 5,6 Ньютона
18 3 3.7 4 (взять 4 последних значения) 6,4 Лагранжа
19 4 1.8 5 (взять 5 первых значений) 7,4 Лагранжа
20 5 7.6 2 4,5 Ньютона
Таблица 2.
x f(x)=1/x·lgx+x2
1,3 1,7776
2,1 4,5634
2,9 8,5694
3,7 13,8436
4,5 20,3952
5,3 28,2267
6,1 37,3387
Таблица 3.
x f(x)=ln2,3x-0,8/x
1,2 0,3486
2,3 1,3180
3,4 1,8214
4,5 2,1592
5,6 2,4128
6,7 2,6156
7,8 2,7845
Таблица 4.
X f(x)=2,1·sin0,37x
-3,2 -1,9449
-0,8 -0,6126
1,6 1,1718
4 2,0913
6,4 1,4673
8,8 -0,2397
11,2 -1,7698
Таблица 5.
x f(x)=1,7-cos(0,4-0,7x)
2,6 2,1874
3,8 3,2888
5 3,9061
6,2 3,8209
7,4 3,2452
8,6 2,6949
9,8 2,6535
3. Вычислить значения интеграла, используя квадратурные формулы:
левых прямоугольников,
правых прямоугольников,
центральных прямоугольников,
трапеции,
Симпсона,
Ньютона,
Гаусса с двумя узлами.
Интеграл вычислить с точностью ε=10-6. Точность вычисления интеграла определяется сравнением результатов при различном числе разбиений отрезка интегрирования. Именно, точность ε считается достигнутой, если

здесь - значение составной квадратурной формулы при разбиении отрезка на N частей.
Отчет должен содержать:
постановку задачи и исходные данные,
описание методов решения и расчетные формулы,
таблицы значений интегралов с указанием числа разбиений, потребовавшихся для достижения заданной точности,
листинг программы.
Варианты заданий.
1. 2. 3. 4.
5. 6. 7. 8.
9. 10. 11. 12.
13. 14. 15. 16.
17. 18. 19. 20.
Индивидуальное домашнее задание № 41. а) Найти первые пять, отличных от нуля, членов разложения решения задачи Коши в ряд Тейлора;
б) методом последовательных приближений найти три последовательные приближения решения задачи Коши, оценить точность.
Варианты:
y'=x2+y2, y(0)=0
y'=y+x, y(0)=1
y'=x2-y2, y(0)=0
y'=x2-y2, y(0)=1
y'=x+y2, y(0)=0
y'=x+y, y(0)=1
y'=2y-2x2-3, y(0)=2
y'=x+y2, y(0)=1
y'=x+y, y(0)=0
y'=x2+y2, y(0)=1
y'=x-y, y(0)=1
y'=4x(1+x), y(0)=1
y'=2x-y, y(0)=1
y'=2x2+y2, y(0)=1
y'=x+2y2, y(0)=1
y'=x2-y, y(0)=2
y'=x+3y, y(0)=-1
y'=x-2y, y(0)=0
y'=2x-y, y(0)=1
y'=y-x, y(0)=1.5
y'=x-3y, y(0)=0
y'=x+2y, y(0)=0
y'=3x+y, y(0)=0
y'=xy, y(0)=0
y'=-x+y, y(0)=1
2. Методом Эйлера (используя все его виды, кроме усовершенствованного метода Эйлера-Коши с итерационной обработкой) найти значение решения задачи Коши на отрезке [1,2] для вариантов 3, 4, 5, 9, 10, 11, 16, 17, для остальных – на отрезке [0,1] с шагом h =0,1 в виде таблицы и построить ломаную Эйлера.
Варианты:
y'=x+y, y(0)=0
y'=2+x, y(0)=0
y'= - y x, y(1)=1
y'=yx+y, y(1)=1
y'=y-3xx+3y, y(1)=1
y'=x+1, y(0)=0
y'=-x-1, y(0)=0
y'=x+2y, y(0)=0
y'=- xy, y(1)=1
y'= y x, y(1)=1
y'=-3x+1, y(0)=1
y'=1+2x, y(0)=0
y'=2x+y, y(0)=0
y'=-2x+1, y(0)=0
y'=1-2x, y(0)=0
y'=- 2xy, y(1)=1
y'= xy, y(1)=1
y'=2x-y, y(0)=0
y'=x-y, y(0)=0
y'=y, y(0)=0
y'=-x+2y, y(0)=0
y'=x+2y, y(0)=0
y'=3x+y, y(0)=0
y'=2+y, y(0)=0
y'=2x+y, y(0)=1
3. Найти с точностью до 0,001 решение дифференциального уравнения первого порядка с указанными начальными условиями на заданном отрезке:
а) методом Эйлера-Коши с итерационной обработкой;
б) метод Рунге-Кутта;
в) методом Кутты-Мерсона;
г) методом Адамса;
д) методом Милна.
Сравнить результаты.
Варианты:
y'=- 1y2- x, y(1)=0, [1;2]
y'= xyx2+y2, y(0)=1, [0;1]
y'=y3ex-2y, y(0)=1, [0;1]
y'-2y=3ex, y(0,3)=1,415, [0,3;0,6]
y'=xy3-y, y(0)=1, [0;1]
y'=x+y2, y(0)=0, [0;0,3]
y'=x3+y2, y(0)=0,5, [0;0,5]
y'=-x2+y2, y(1)=1, [1;2]
y'=2x+cosy, y(0)=0, [0;0,1]
y'=x2+y2, y(0)=0,27, [0;1]
y'=ex-y2, y(0)=0, [0;0,4]
y'+xy(1-y2)=0, y(0)=0,5, [0;1]
y'=x+sin(y/3), y(0)=1, [0;2]
y'=x2-xy3+y2, y(0)=0,1, [0;1]
y'=2xy+x, y(0)=0, [0;0,5]
y'= 2y-xy, y(1)=2, [1;2]
y'=xy+ey, y(0)=0, [0;0,1]
y'=x2+xy+y2+1, y(0)=0, [0;1]
y'=x3-y, y(1)=-1, [1;2]
4y'=y2+4x2, y(0)=1, [0;1]
y'=xy+0,5y, y(0)=1, [0;1]
y'=y- 2xy, y(0)=1, [0;1]
x2y'-xy=1, y(1)=0, [1;2]
y'=x2+y2, y(0)=0, [0;1]
y'=x+y, y(0)=1, [0;1]
ЛитератураБахвалов Н. С. Численные методы: учеб. пособие для вузов/ Н. С. Бахвалов, Н. П. Жидков, Г. М. Кобельков. - 4-е изд.; Гриф МО. - Москва: БИНОМ. Лаборатория знаний, 2006. - 636 с. - (Классический университетский учебник). - Библиогр.: с. 624-628. - Предм. указ.: с. 629-632.
Вержбицкий В. М. Численные методы: Линейная алгебра и нелинейные уравнения: учеб. пособие для вузов / В. М. Вержбицкий. - 2-е изд., испр.; Гриф МО. - Москва: ОНИКС 21 век, 2005. - 431 с.: ил. - Библиогр.: с. 419-424. - Предм. указ. с. 425-429.
Демидович Б. П. Численные методы анализа: приближение функций, дифференциальные и интегральные уравнения: учеб. пособие / Б. П. Демидович, В. А. Марон, Э. З. Шувалова; под ред. Б. П. Демидовича. - Изд. 4-е, стер. - Санкт-Петербург: Лань, 2008. - 400 с. - (Классическая учебная литература по математике). - Библиогр. в конце гл.
Киреев В. И. Численные методы в примерах и задачах: учеб. пособие для вузов/ В. И. Киреев, А. В. Пантелеев. - Изд. 2-е, стер.; Гриф УМО. - Москва: Высш. шк., 2006. - 480 с.: ил. - (Прикладная математика для ВТУЗов). - Библиогр.: с. 477-480 .
Пантина И. В. Вычислительная математика [Электронный ресурс] учебник для вузов/ И. В. Пантина, А. В. Синчуков. - 2-е изд., перераб. и доп. - Москва: Синергия, 2012. - 175 с.
Пирумов У. Г. Численные методы: учеб. пособие для вузов/ У. Г. Пирумов. - 2-е изд., испр. и доп.; гриф МО. - Москва: Дрофа, 2003. - 221 с.: ил. - (Высшее образование). - Библиогр.: с. 216. - Имен. указ.: с. 217.
Срочко В. А. Численные методы: курс лекций / В. А. Срочко. - Гриф УМО. - Санкт-Петербург [и др.]: Лань, 2010. - 202 с. - Библиогр.: с. 200.
Турчак Л. И. Основы численных методов: учеб. пособие / Л. И. Турчак, П. В. Плотников. - Изд. 2-е, перераб. и доп.; Гриф МО. - Москва: ФИЗМАТЛИТ, 2005. - 300 с.: ил. - Библиогр.: с. 290-292. - Прил.: с. 286-289. - Предм. указ.: с. 293-300.
Шевцов Г. С. Численные методы линейной алгебры: учеб. пособие для мат. напр. и спец. / Г. С. Шевцов, О. Г. Крюкова, Б. И. Мызникова. - Гриф УМО. - Москва: Финансы и статистика: ИНФРА-М, 2008. - 479 с. - (Финансы и статистика). - Библиогр.: с. 473-474. - Предм. указ.: с. 475-479.

Приложенные файлы

  • docx 23676701
    Размер файла: 4 MB Загрузок: 0

Добавить комментарий