Математические основы в энергетике

Математические основы в энергетике Курс лекций Оглавление 1. Элементы теории погрешностей 4 2. Приближенное решение нелинейных и трансцендентных уравнений 8 2.1. Постановка задачи 8 2.2. Графическое решение уравнений 8 2.3. Метод половинного деления (дихотомии) 9 2.4. Метод хорд 10 2.5. Метод Ньютона (метод касательных) 11 2.6. Комбинированный метод 13 3. Численное решение систем линейных алгебраических уравнений (СЛАУ) 15 3.1. Решение задач линейной алгебры 15 3.2. Метод Гаусса 17 3.3. Схема Гаусса с выбором главного элемента 21 3.4. Вычисление обратной матрицы методом Гаусса 22 3.5. Вычисление определителей методом Гаусса 24 3.6. Метод простой итерации (метод Якоби) 26 3.7. Метод Зейделя 27 3.8. Метод скорейшего спуска (градиента) для случая системы линейных алгебраических уравнений 29 4. Приближенное решение систем нелинейных уравнений 34 4.1. Метод Ньютона 34 4.2. Метод градиента (метод скорейшего спуска) 37 5. Интерполирование и приближение функций 40 5.1. Задача интерполирования и аппроксимации функций 40 5.2. Интерполирование алгебраическими многочленами 40 5.3. Интерполяционная формула Ньютона 41 5.4. Сходимость интерполяционного процесса 43 5.5. Задача обратного интерполирования 44 5.6. Отыскание параметров эмпирических формул методом наименьших квадратов 45 5.7. Суть метода наименьших квадратов 45 6. Приближенное дифференцирование и интегрирование функции одной переменной … 50 6.1. Численное интегрирование 50 6.1.1. Постановка задачи 50 6.1.2. Формула прямоугольников 51 6.1.3. Формула трапеций 52 6.1.4. Формула Симпсона 53 6.1.5. Вычисление определенных интегралов методами Монте-Карло 55 6.2. Решение обыкновенных дифференциальных уравнений 58 6.2.1. Методы решения задачи Коши 58 6.2.2. Метод рядов, не требующий вычисления производных правой части уравнения 60 6.2.3. Метод Рунге-Кутта 61 6.2.4. Многошаговые методы 63 6.2.5. Экстраполяционные методы Адамса 64 6.2.6. Интерполяционные методы Адамса 65 7. Статистические и динамические модели. Метод монте карла……………………………67 7.1. Статистические и динамические модели. …………………..……………………………67 7.2. Метод монте карла…………………………………………………………………………68 Литература 69 1. Элементы теории погрешностей Численное решение любой задачи, как правило, осуществляется приближенно, с различной точностью. Это может быть обусловлено неточностью исходных данных, конечной разрядностью вычислений (вручную или на ЭВМ) и т. п. Главная задача численных методов – фактическое нахождение решения с требуемой или, по крайней мере, оцениваемой точностью. Отклонение истинного решения от приближенного называется погрешностью. Полная погрешность вычислений состоит из двух составляющих: 1) неустранимая погрешность; 2) устранимая погрешность. Неустранимая погрешность обусловлена неточностью исходных данных и никаким образом не может быть уменьшена в процессе вычислений. Устранимая погрешность состоит из двух составляющих: а) погрешность аппроксимации (метода); б) погрешность вычислений. Эти составляющие могут быть уменьшены выбором более точных методов и увеличением разрядности вычислений. Как правило, в дальнейшем нас будут интересовать корректно поставленные задачи вычисления. Задача вычисления y = A(x) называется корректно поставленной, если для любых входных данных из некоторого класса решение задачи существует, единственно и устойчиво по входным данным (т. е. непрерывно зависит от входных данных задачи). В сформулированном понятии корректности учтены достаточно естественные требования, т. к. чтобы численно решать задачу, нужно быть уверенным, что ее решение существует. Столь же естественны требования единственности и устойчивости решения. Рассмотрим подробнее понятие устойчивости. Обычно нас интересует решение y, соответствующее входным данным x. Реально мы имеем возмущенные входные данные с погрешностью δx, т.е. x + δx, и находим возмущенное решение: y + δy = A(x+δx). Эта погрешность входных данных порождает неустранимую погрешность решения: δy = A(x+δx) - A(x). Если решение непрерывно зависит от входных данных, то всегда при , и задача устойчива по входным данным. Здесь символ || || - норма. Если X – точное значение величины, а X* – приближенное значение, то абсолютная погрешность приближения определяется выражением . Относительной погрешностью величины называется отношение абсолютной погрешности к модулю ее точного значения: δ = Δ / |X|. Достаточно часто точное значение величины неизвестно, поэтому указывают границы погрешности: (1.1) (1.2) Рассмотрим подробнее погрешность округления чисел, участвующих в вычислениях. В позиционной системе счисления с основанием r запись (1.3) означает, что . Здесь r – целое число, большее единицы. Каждое из чисел может принимать одно из значений {0, 1, …, r-1}. Числа называются разрядами, например: – третий разряд перед запятой, – второй разряд после запятой. Запись вещественного числа в виде (1.3) называется также его представлением в форме числа с фиксированной запятой. В ЭВМ чаще всего используется представление чисел в форме с плавающей запятой. Так как наиболее часто в компьютерах применяется двоичная система с плавающей запятой, то вещественное число можно представить виде (1.4) где Здесь p - целое число называется порядком числа a, а – мантиссой. Если исходная величина или промежуточный результат требуют большего числа разрядов, то производится округление до t – го разряда. Значащие цифры называются верными до t – го разряда, если абсолютная величина разности между a* и a меньше или равна половине единицы младшего разряда: . (1.5) Ограничения на порядки чисел, представляемых в ЭВМ , порой приводят к прекращению вычислений (так называемое исчезновение порядка); в других случаях относительно небольшая разрядность представления чисел в ЭВМ приводит к недопустимому искажению результата вычислительной погрешностью. Приведем несколько примеров иллюстрирующих это и способы (приемы) уменьшения вычислительной погрешности за счет несложных алгебраических преобразований. Рассмотрим типичный пример, в котором порядок выполнения операций существенно влияет на погрешность результата. Пример 1.1. Необходимо отыскать минимальный корень уравнения. Вычисления производим в десятичной системе счисления, причем в числе после округления оставляем четыре действующие цифры (разряда): Рассмотрим другой алгоритм вычисления корня, для чего избавимся от иррациональности в числителе: Как видно из сравнения полученных результатов, применение "неудач­ного" алгоритма завышает результат на 30 %. Это явление в прикладной математике (в практике вычислений) называется потерей значащих цифр, и часто наблюдается при вычитании близких величин. Потеря значащих цифр, например, довольно часто приводит к существенному искажению результатов при решении даже сравнительно небольших систем линейных алгебраических уравнений. Пример 1.2. На машине с плавающей запятой необходимо вычислить значение суммы . Эту сумму можно вычислить двумя способами: 1. 2. Оказывается, для второго алгоритма вычислительная погрешность будет существенно меньше. Тестовые расчеты на конкретной ЭВМ по первому и второму алгоритмам показали, что величина погрешности для обоих алгоритмов составляет 2∙10-4 и 6∙10-8соответственно. Причина этого ясна, если вспомним, как числа представлены в ЭВМ (см. рис. 1.1). Рис. 1.1. Представление чисел с плавающей точкой в типичном 32-бит (4 - байт) формате: а) Число ½; b) число 3; c) число ¼; d) число 10-7, представленное в машине с максимальной точностью (нормализованное); e) то же самое число 10-7, но представленное с тем же порядком, что и число 3; f) сумма чисел 3 + 10-7, которая эквивалентна 3. Этот пример ярко иллюстрирует тот факт, что даже если оба слагаемых представлены в компьютере точно, то их сумма может быть представлена с погрешностью, особенно если слагаемые различаются на порядки. При оценке погрешностей арифметических действий следует учитывать следующее: а) абсолютная погрешность алгебраической суммы не превышает суммы абсолютных погрешностей ее членов: (1.6) б) относительная погрешность произведения не превышает суммы относительных погрешностей сомножителей: ; (1.7) в) относительная погрешность частного не превышает суммы относительных погрешностей делимого и делителя: (1.8) Вопросы для самопроверки • Дайте определения и приведите примеры устранимой и неустранимой погрешностей. • Что такое погрешность округления? Какова ее связь с разрядностью ЭВМ? • Как вычислить относительную погрешность, зная абсолютную? • Как по абсолютной погрешности вычислить относительную погрешность? 2. Приближенное решение нелинейных и трансцендентных уравнений 2.1. Постановка задачи Пусть дано уравнение f(x) = 0, (2.1) где функция f(x) определена и непрерывна в конечном или бесконечном интервале a < x < b. Всякое значение ξ, обращающее функцию f(x) в нуль, то есть такое, что f(ξ) = 0, называется корнем уравнения (2.1) или нулем функции f(x). Предположим, что уравнение (2.1) имеет лишь изолированные корни, то есть для каждого корня существует окрестность, не содержащая других корней этого уравнения. Приближенное нахождение изолированных действительных корней уравнения (2.1) складывается обычно из двух этапов: 1. Отделение корней, то есть установление возможно тесных промежутков [α, β], в которых содержится один и только один корень исходного уравнения (2.1). 2. Уточнение приближенных корней, то есть доведение их до заданной степени точности. 2.2. Графическое решение уравнений Действительные корни уравнения f(x) = 0 приближенно можно определить как абсциссы точек пересечения графика функции y = f(x) с осью ОХ (см. рис. 2.1, а). На практике часто бывает удобнее уравнение (2.1) заменить равносильным ему уравнением , (2.2) где функции φ(x) и ψ(x) более простые, чем функция f(x). Тогда, построив графики этих функций, искомые корни получим как абсциссы точек пересечения этих графиков (смотри рис. 2.1, б). а) б) Рис. 2.1. Графический метод нахождения корней уравнения. 2.3. Метод половинного деления (дихотомии) Сформулируем без доказательства очень важную для рассмотрения дальнейших вопросов теорему. Теорема: Если непрерывная функция f(x) принимает значения разных знаков на концах отрезка [α, β], то есть f(α)·f(β) < 0, то внутри этого отрезка содержится по меньшей мере один корень уравнения f(x) = 0, а именно: найдётся хотя бы одно число такое, что f(ξ) = 0. Пусть дано уравнение f(x) = 0, (2.3) где функция f(x) определена и непрерывна на интервале [a, b] и f(a)·f(b) < 0. Для нахождения корня уравнения делим отрезок [a, b] пополам: • если f((a + b)/2) = 0, то ξ = (a + b)/2 является корнем уравнения (2.3); • если , то выбираем ту половину отрезка [a, (a + b)/2] или [(a + b)/2, b], на концах которого функция f(x) имеет противоположные знаки. Новый суженный отрезок [a1, b1] снова делим пополам и проводим тот же анализ и т.д. Очевидно, что закончить уточнение значения корня можно при достижении условия |аj – bj| < ε , где ε > 0 - сколь угодно малое число. Второй способ закончить вычисления - задать максимальное значение невязки: f((aj + bj)/2) < ε. Замечания • Метод половинного деления очень прост, здесь нет вычислительной формулы и можно обеспечить практически любую точность. • Как недостаток метода можно отметить его медленную сходимость (за один шаг интервал, где находится корень, сужается всего в два раза). 2.4. Метод хорд Пусть дано уравнение f(x) = 0, (2.4) где функция f(x) определена и непрерывна на интервале [a, b] и выполняется соотношение f(a)·f(b) < 0. Пусть для определенности f(a) < 0, f(b) > 0. Тогда вместо того, чтобы делить отрезок [a, b] пополам, более естественно разделить его в отношении - f(a):f(b). При этом новое значение корня определяется из соотношения x1 = a + h1, (2.5) где . (2.6) Далее этот прием применяем к одному из отрезков [a, x1] или [x1, b], на концах которого функция f(x) имеет противоположные знаки. Аналогично находим второе приближение x2 и т.д. (см. рис. 2.2.). Геометрически этот способ эквивалентен замене кривой y = f(x) хордой, проходящей через точки А(a, f(a)) и B(b, f(b)). Рис. 2.2. Уточнение корня уравнения методом хорд Действительно, уравнение хорды АВ имеет вид (2.7) Учитывая, что при х = х1 => y = 0, получим (2.8) Полагая, что на отрезке [a, b] вторая производная f''(x) сохраняет постоянный знак, метод хорд сводится к двум различным вариантам: 1. Из рис. 2.2,a видно, что неподвижна точка а, а точка b приближается к ξ, то есть (2.9) Преобразовав выражение (2.9), окончательно получим (2.10) 2. Из рис. 2.2,b видно, что точка b остается неподвижной, а точка а приближается к ξ, тогда вычислительная формула примет вид (2.11) Таким образом, для вычисления корня уравнения имеем две различные вычислительные формулы (2.10) и (2.11). Какую точку брать за неподвижную? Рекомендуется в качестве неподвижной выбирать ту точку, в которой выполняется соотношение f(x)·f”(x) > 0. (2.12) 2.5. Метод Ньютона (метод касательных) Пусть корень ξ уравнения f(x) = 0, (2.13) отделен на отрезке [a, b], причем первая и вторая производные f(x) и f(x) непрерывны и сохраняют определенные знаки при . Найдя какое-нибудь n-ое приближение корня , мы можем уточнить его по методу Ньютона следующим образом. Пусть ξ = xn + hn, (2.14) где hn - величина малая. Отсюда по формуле Тейлора получим (ограничиваясь первым порядком малости относительно hn) f(xn + hn) = f(xn) + hn f(xn) = 0. (2.15) Следовательно, hn = - f(xn) / f (xn). (2.16) Подставив полученное выражение в формулу (2.14), найдем следующее (по порядку) значение корня: (2.17) Проиллюстрируем графически нахождение корня методом Ньютона (рис. 2.3.). Рис. 2.3. Уточнение корня методом касательных Если в качестве начального приближения выбрать точку х0 = В0 , то процесс быстро сходится. Если же выбрать точку х0 = А0, то х1 [a, b], и процесс нахождения корня расходится. Рекомендуется: в качестве х0 выбрать точку, где f(x)·f(x) > 0. 2.6. Комбинированный метод Пусть f(a)·f(b) < 0, а f(x) и f(x) сохраняют постоянные знаки на отрезке [a¸ b]. Соединяя метод хорд и метод касательных, получаем метод, на каждом шаге которого находим значения по недостатку и значения по избытку точного корня ξ уравнения f(x) = 0. Теоретически здесь возможны четыре случая: • f(x) > 0; f(x) > 0; • f(x) > 0; f(x) < 0; • f(x) < 0; f(x) > 0; • f(x) < 0; f(x) < 0. Рассмотрим только первый случай, так как остальные три ведут себя аналогично и могут быть сведены к первому. Итак, пусть f(x) > 0 и f(x) > 0 при . Полагаем, что (для метода хорд), (для метода касательных). Тогда новые значения корня вычисляем по формулам (2.18) Рис. 2.4 наглядно иллюстрирует суть комбинированного метода. Рис. 2.4. Уточнение корня комбинированным методом Доказано, что . Следует обратить внимание на то, что на каждом шаге метод хорд применяется к новому отрезку . Если задать максимальное значение погрешности ε > 0, процесс уточнения значения корня продолжаем до тех пор, пока не выполнится условие . (2.19) Пример 2.1. Вычислить с точностью до 0.0005 положительный корень уравнения f(x) = x5 – x – 0.2 = 0. На первом этапе отделения корней выбрали интервал [1.0, 1.1], на концах которого функция имеет противоположные знаки. Действительно, f(1) = – 0.2 < 0, f(1.1) = 0.31051 > 0. В выбранном нами интервале f(x) > 0, f(x) > 0, то есть знаки производных сохраняются. Применим комбинированный метод, приняв . По формулам (2.18) вычислим . Так как точность недостаточная (погрешность велика), вычислим следующие значения: Таким образом, за два шага мы обеспечили требуемую точность. Замечания • Комбинированный метод наиболее трудоемок. • Метод, как и метод Ньютона не всегда сходится (почему?). • Комбинированный метод сходится быстрее всех ранее рассмотренных, (если он сходится). Вопросы для самопроверки • Какие точные методы решения нелинейных уравнений вы знаете? • Для чего нужен первый этап - отделение корней? • Сформулируйте условия существования решения уравнения. Являются ли эти требования необходимыми и достаточными? • Что можно сказать о точности методов половинного деления, хорд, касательных и комбинированного? По каким параметрам их еще можно сравнить? • В соответствии с известной теоремой на отрезке [a, b] существует решение. Всегда ли его можно найти методом половинного деления, методом хорд, и т.п.? 3. Численное решение систем линейных алгебраических уравнений (СЛАУ) 3.1. Решение задач линейной алгебры Линейные системы имеют в вычислениях очень большое значение, так как к ним может быть приведено приближенное решение широкого круга задач. Так, основными источниками возникновения СЛАУ являются теория электрических цепей, уравнения балансов и сохранения в механике, гидравлике и т.п. Пусть дана система n линейных алгебраических уравнений с n неизвестными: (3.1) Или в матричной форме: ; (3.2) где (3.3) - матрица коэффициентов системы (3.1); - вектор неизвестных; - вектор свободных членов. Если матрица A неособенная, т.е. (3.4) то система (3.1) или эквивалентное ей матричное уравнение (3.2) имеют единственное решение. Действительно, при условии, что detA  0, существует обратная матрица A-1. Умножая обе части уравнения (3.2) слева на A-1, получим: (3.5) Формула (3.5) даёт решение уравнения (3.2), причём единственное. Пример 3.1. Для матрицы A порядка n > 4 непосредственное нахождение обратной матрицы A-1 требует много времени (операций). Поэтому формула (3.5) на практике употребляется достаточно редко. Обычно значения неизвестных xi (i = 1,2, ... n) могут быть получены по известным формулам Крамера: (3.6) Здесь матрица Ai получается из матрицы A заменой её i-го столбца столбцом свободных членов. Пример 3.2. Решим вышеприведенную систему по формулам Крамера: Применяемые в настоящее время методы решения СЛАУ можно разбить на две группы: точные и приближённые. Точными методами называются такие методы, которые в предположении, что вычисления ведутся точно (без округлений), за конечное число действий позволяют получить точные значения неизвестных xi. Приближенными методами называются такие методы, которые даже в предположении, что вычисления ведутся без округлений, позволяют получить решение системы (x1, x2, ..., xn) лишь с заданной точностью. Точное решение СЛАУ в этих случаях может быть получено теоретически как результат бесконечного процесса. К приближенным методам относятся метод простой итерации, метод Зейделя и т.п. 3.2. Метод Гаусса Наиболее распространенным методом решения СЛАУ является метод Гаусса, в основе которого лежит идея последовательного исключения неизвестных. Существуют различные схемы, реализующие данный метод. Рассмотрим одну из них – схему единственного деления. Для простоты ограничимся рассмотрением СЛАУ с четырьмя неизвестными: (3.7) Пусть a11  0 (ведущий элемент). Разделив первое уравнение на a11, получим первую главную строку: (3.8) где (j = 2,3,4,5). Используя уравнение (3.8), можно исключить неизвестные x1 из 2-го, 3-го и 4-го уравнений системы (3.7). Для этого последовательно умножаем уравнение (3.8) на a21; a31; a41 и вычитаем результат из 2-го, 3-го и 4-го уравнений системы (3.7) соответственно. В результате получим систему из трех уравнений: (3.9) где коэффициенты вычисляются по формуле (i = 2, 3, 4; j = 2, 3, 4, 5). (3.10) Далее первое уравнение системы (3.9) делим на ведущий элемент и получаем (3.11) где , (j = 3, 4, 5). Аналогично предыдущему шагу, исключая x2, как и x1, получим систему (3.12) Здесь (i = 3, 4; j = 3, 4, 5). Разделив первое уравнение системы (3.12) на , получим: (3.13) где (j = 4, 5). Теперь с помощью уравнения (3.13) исключим x3 из второго уравнения системы (3.12), окончательно получим: , (3.14) где (j=4, 5). Таким образом, исходную систему (3.7) привели к составленной из главных строк (3.8), (3.11), (3.13) и (3.14) эквивалентной системе с треугольной матрицей(3.15): (3.15) Из (3.15) последовательно находим (3.16) Итак, решение СЛАУ (3.7) распадается на два этапа: • прямой ход (приведение системы (3.7) к треугольному виду (3.15)); • обратный ход (определение неизвестных по формуле (3.16)). Пример 3.3. Прямой ход: Из выражений (3.10) вычислим коэффициенты : Аналогично вычислим коэффициенты при (i = 3, 4) и составим систему Разделив первое уравнение системы на , получим Значит, Из (3.12) вычислим для i = 3 и j = 3, 4, 5: Аналогично, вычислив коэффициенты для i = 4, получим: Разделив первое уравнение на a(2)33 = 16.425, получим: где По формуле (3.14) находим коэффициенты : и записываем одно уравнение с одним неизвестным: 1.1199786x4 = -1.1199768. x1 + 0.5x2 - 0.05x3 + 0.5x4 = 1.35; x2 + 13.4x3 - 29x4 = 71.2; x3 - 1.72298x4 = 4.72298; 1.11998x4 = -1.11998. На этом закончен прямой ход. Обратный ход: x4 = -1.000; x3 = 4.72298 - 1.72298 = 3; x2 = 71.2 - 13.4 * 3-29 = 2; x1 = 1.35 - 0.5 * 2 + 0.05 * 3 + 0.5 = 1. 3.3. Схема Гаусса с выбором главного элемента Рассмотрим СЛАУ (3.17) Запишем расширенную прямоугольную матрицу коэффициентов системы (3.17): . (3.18) Среди элементов матрицы aij (i,j = 1, ...n) выберем наибольший по модулю, называемый главным элементом. Пусть им будет, например, элемент apq. Строка, содержащая главный элемент, называется главной строкой. Далее вычисляем множители mi = aiq / apq для всех i  p.Затем преобразуем матрицу (3.18) следующим образом: из каждой i-ой неглавной строки вычитаем почленно главную строку, умноженную на mi. В результате получим матрицу, у которой все элементы q-го столбца за исключением apq, равны 0. Отбрасывая этот столбец и главную строку, получим новую матрицу M1 с числом строк и столбцов на 1 меньше. Над матрицей М1 повторяем те же операции, после чего получим матрицу M2 и т.д. Таким образом продолжаем до тех пор, пока не получим матрицу, содержащую одну строку из двух элементов, которую тоже считаем главной. Затем объединим все главные строки, начиная с последней. После некоторой перестановки они образуют треугольную матрицу, эквивалентную исходной. На этом заканчивается этап вычислений, называемый прямым ходом. Решив систему с полученной треугольной матрицей коэффициентов, найдём последовательно значения неизвестных xi (i = 1, 2, ..., n). На этом заканчивается обратный ход. Смысл выбора главного элемента состоит в том, чтобы сделать возможно меньшими числа mi и тем самым уменьшить погрешность вычислений. Пример 3.4. Рассмотрим СЛАУ, состоящую из трех уравнений. Запишем расширенную матрицу m2 = -1/6; m3 = -2/3. m2 = -5/16. M2 = [ 87/96 174/32]. x3 = 6; x1 = 3; x2 = -2. 3.4. Вычисление обратной матрицы методом Гаусса Пусть дана неособенная матрица A = [aij] (i,j = 1,2, ..., n). (3.19) Необходимо найти её обратную матрицу A-1 = [xij] (i,j = 1,2, ..., n). (3.20) Вспомним основное соотношение линейной алгебры: A·A-1 = E, (3.21) где Е – единичная матрица. Перемножая матрицы A и A-1, получаем n2 уравнений относительно n2 неизвестных xij: (i,j = 1, 2, ..., n), (3.22) где Таким образом, получим n систем линейных уравнений для j = 1, 2, ..., n, имеющих одну и ту же матрицу коэффициентов A и различные столбцы - свободные члены, которые можно одновременно решить методом Гаусса. Рассмотрим это подробнее, вычислив матрицу, обратную : Разделив все коэффициенты первой строки на a11 = 2, получим первую главную строку (обратите внимание, что с n столбцами свободных членов проводятся те же действия, что и с одним): 1.0 0.5 -0.05 0.5 0.5 0 0 0 1.0 13.4 -29 -0.6667 3.333 0 0 . Для проверки перемножим полученную обратную матрицу и исходную (должны получить единичную): . Благодаря округлению, убеждаемся, что обратная матрица вычислена неточно. В дальнейшем можно показать, как методом простой итерации можно уточнить A-1. 3.5. Вычисление определителей методом Гаусса Пусть дана исходная матрица . (3.23) Необходимо вычислить  = det A. Вспомним свойства определителей: • для того чтобы умножить (разделить) определитель на какое либо число, достаточно умножить (разделить) на это число строку или столбец: ; (3.24) • значение определителя не изменится, если его строку заменить суммой этой строки и другой, умноженной на произвольное число. Учитывая это свойство, умножая первую строку последовательно на a21, a31, ..., an1 и вычитая из второй, третьей и т.д., получим ; (3.25) • величина определителя равна сумме произведений элементов строки (столбца) на (-1)i+j |A|i j, где |A|i j – соответствующие миноры. Используя это свойство, представим определитель как сумму произведений элементов первого столбца на соответствующие миноры. При этом учтем, что за исключением первого элемента значения остальных элементов столбца равны нулю. Таким образом, мы понизили порядок определителя на 1. Применим к полученному определителю порядка n - 1 такие же преобразования. Выполняя n шагов, найдем определитель  как произведение ведущих элементов: (3.26) Пример 3.5. Замечания • При наличии решения, точные методы всегда дадут его через конечное число шагов. • В рамках точных методов вычислительная погрешность увеличивается с ростом размеров СЛАУ и не может быть уменьшена. 3.6. Метод простой итерации (метод Якоби) Рассмотрим систему A·x = f, (3.27) где матрица A = [aij] (i,j = 1, 2, …m) имеет обратную матрицу; x = (x1, x2, x3,…xm) – вектор неизвестных, f – вектор свободных членов. Преобразуем систему (3.27) к следующему виду: (i = 1, 2, …m), (3.28) где , , при этом предполагаем, что . Условимся, как обычно, считать значение суммы равным нулю, если верхний предел суммирования меньше нижнего. Тогда при i = 1 уравнение (3.28) имеет вид (3.29) В методе простой итерации (методе Якоби) исходят из записи системы в виде (3.28), итерации при этом определяют следующим образом: (3.30) Начальные значения – (i = 0, 1, … m) задаются произвольно. Окончание итерационного процесса определяют либо заданием максимального числа итераций n0, либо следующим условием: (3.31) где ε > 0. В качестве нулевого приближения в системе (3.30) примем . (3.32) Если последовательность приближений x1(0), x2(0), ..., xm(0), x1(1), x2(1), ..., xm(1), ..., x1(k), x2(k), ..., xm(k) имеет предел , (3.33) то этот предел является решением системы (3.28). Достаточным условием сходимости решения системы (3.27) является то, что матрица A является матрицей с преобладающими диагональными элементами, то есть (3.34) 3.7. Метод Зейделя Этот метод представляет собой некоторую модификацию метода простой итерации. Основная его идея заключается в том, что при вычислении (k+1)-го приближения неизвестной xi учитываются уже вычисленные ранее (k+1)-е приближения (x1 x2, ..., xi-1). Пусть дана приведенная линейная система: (i = 1, 2, …n). (3.35) Выберем произвольно начальные приближения корней , стараясь, конечно, чтобы они в какой-то мере соответствовали неизвестным x1, x2, x3, ..., xn. Предположим, что k-е приближение корней известно, тогда в соответствии с идеей метода будем строить (k+1) – е приближение по следующим формулам: (3.36) (k = 0, 1, 2,...). Обычно процесс Зейделя сходится быстрее, чем метод Якоби. Бывает, что процесс Зейделя сходится, когда простая итерация расходится и, т.п. Правда, бывает и наоборот. Во всяком случае, достаточные условия сходимости для метода Якоби достаточны и для сходимости метода Зейделя. Если выполняется достаточное условие сходимости для системы (3.35) – по строкам, то в методе Зейделя выгодно расположить уравнения (3.36) так, чтобы первое уравнение системы имело наименьшую сумму модулей коэффициентов: . (3.37) Пример 3.6. Для того чтобы обеспечить достаточные условия сходимости итерационного процесса (преобладающие значения диагональных элементов), преобразуем исходную систему и приведем к удобному виду. Чтобы дальнейшие преобразования были понятны, обозначим уравнения исходной системы буквами А, Б, В и Г соответственно: х1= -0.2х2 +0.1х3 – 0.2х4 – 0.4; (Г) х2 = -0.2х1 – 0.2х3 + 0.2; (А – Б) х3 = 0.2х1 – 0.4х2 + 0.2х4 – 0.4; (Б) х4 = 0.333х1 - 1.111. (2А – Б + 2В – Г) Преобразованную систему будем решать методом Зейделя, тогда, с учетом требования (3.37), окончательно получим: В качестве нулевого приближения (k = 0) возьмем . Зададим количество итераций k = 2 и все результаты вычислений сведем в табл. 3.1. Таблица 3.1 Итерация, k Значения неизвестных Невязки x1 x2 x3 x4 ε1 ε2 ε3 ε4 -0.4 0.2 -0.4 -1.111 -2.711 -1.911 0.444 -1.422 1 -0.263 0.36 -0.846 -1.244 -0.309 1.0 0.734 0.446 2 -0.329 0.422 -0.874 -1.199 0.095 -0.000 0.009 0.029 В приведенной таблице кроме значений неизвестных на каждом шаге оценивались невязки. Вспомним, что корнями уравнения называются такие значения неизвестных, которые превращают его в тождество. Так как мы используем итерационный (приближенный) метод, значения неизвестных вычисляем приближенно (три, четыре знака после десятичной точки), то, подставляя значения неизвестных в исходную систему, справа получим не ноль, а некоторые значения, называемые невязкой первого, второго, … уравнений на k –ом шаге. Анализ данных, приведенных в табл. 3.1, показывает, что итерационный процесс быстро сходится, о чем свидетельствуют как быстрое уменьшение невязок, так и уменьшение изменений неизвестных (см. формулу (3.31) метода Якоби). 3.8. Метод скорейшего спуска (градиента) для случая системы линейных алгебраических уравнений В рассматриваемом ниже итерационном методе вычислительный алгоритм строится таким образом, чтобы обеспечить минимальную погрешность на шаге (максимально приблизиться к корню). Представим систему линейных уравнений в следующем виде: (3.38) Запишем выражение (3.38) в операторной форме: (3.39) Здесь приняты следующие обозначения: ; ; . (3.40) В методе скорейшего спуска решение ищут в виде , (3.41) где и - векторы неизвестных на p и p+1 шагах итераций; вектор невязок на p-ом шаге определяется выражением , (3.42) а (3.43) В формуле (3.43) используется скалярное произведение двух векторов, которое определяется следующей формулой: (3.44) В формуле (3.43) - транспонированная матрица Якоби, вычисленная на p-ом шаге. Матрица Якоби вектор – функции f(x) определяется как (3.45) Нетрудно убедиться, что для системы (3.39) матрица Якоби равна (3.46) Как и для метода простой итерации, достаточным условием сходимости метода градиента является преобладание диагональных элементов. В качестве нулевого приближения можно взять . Замечания • Как видно из выражения (3.45), матрица Якоби не зависит от шага итерации. • Требования минимизации погрешности на каждом шаге обусловили то, что метод градиента более сложен (трудоемок), чем методы Якоби и Зейделя. • В методе градиента итерационный процесс естественно закончить при достижении , вектор невязок входит в вычислительную формулу. Обсуждение • В приближенных методах можно обеспечить практически любую погрешность, если итерационный процесс сходится. • Итерационный процесс можно прервать на любом k–ом шаге и продолжить позднее, введя в качестве нулевого шага значения x(k). • В качестве недостатка приближенных методов можно отметить то, что они часто расходятся, достаточные условия сходимости (преобладание диагональных элементов) можно обеспечить только для небольших систем из 3 – 6 уравнений. Пример 3.7. Методом скорейшего спуска решим систему уравнений Так как диагональные элементы матрицы являются преобладающими, то в качестве начального приближения выберем: Следовательно, вектор невязок на нулевом шаге равен Далее последовательно вычисляем Отсюда причем Аналогично находятся последующие приближения и оцениваются невязки. Что касается данного примера, можно отметить, что итерационный процесс сходится достаточно медленно (невязки уменьшаются). Вопросы для самопроверки • Назовите известные вам методы решения СЛАУ. • Чем точные методы отличаются от приближенных? • Что такое прямой и обратный ход в методе Гаусса? • Нужен ли обратный ход при вычислении методом Гаусса а) обратной матрицы; б) определителя? • Что такое невязка? • Сравните достоинства и недостатки точных и приближенных методов. • Что такое матрица Якоби? • Надо ли пересчитывать матрицу Якоби на каждом шаге итерации в методе градиента? • Исходная СЛАУ решается независимо тремя методами – методом Якоби, методом Зейделя и методом градиента. Будут ли равны значения а) начального приближения (нулевой итерации); б) первой итерации? • При решении СЛАУ (n > 100) итерационными методами решение расходится. Как найти начальное приближение? 4. Приближенное решение систем нелинейных уравнений 4.1. Метод Ньютона Рассмотрим нелинейную систему уравнений (4.1) с действительными левыми частями. Систему (4.1) можно представить в матричном виде (4.2) Здесь приняты следующие обозначения: - вектор аргументов, а - вектор – функция. Для решения системы (4.2) воспользуемся методом последовательных приближений. Предположим, что найдено р-ое приближение xp = (x1(p), x2(p) , ..., xn(p)) одного из изолированных корней x = (x1, x2, x3, ..., xn) векторного уравнения (4.2). Тогда точный корень уравнения (4.2) можно представить в виде (4.3) где - поправка (погрешность) корня на n – ом шаге. Подставив выражение (4.3) в (4.2), получим (4.4) Предположим, что функция f(x) - непрерывно дифференцируема в некоторой выпуклой области, содержащей x и x(p). Тогда левую часть уравнения (4.4) разложим в ряд Тейлора по степеням малого вектора ε(p), ограничиваясь линейными членами: , (4.5) или в развернутом виде: (4.6) Из анализа формул (4.5) и (4.6) следует, что под производной f(x) следует понимать матрицу Якоби системы функций f1 , f2, ..., fn, относительно переменных x1, x2, x3, ..., xn, то есть: . (4.7) Выражение (4.7) в краткой записи можно представить: (4.8) Выражение (4.6) представляет собой линейную систему относительно поправок (i = 1, 2, ..., n) с матрицей W(x), поэтому формула (4.5) может быть записана в следующем виде: (4.9) Отсюда, предполагая, что матрица W(x(p)) - неособенная, получим: (4.10) Теперь, подставив выражение (4.10) в формулу (4.3), окончательно получим: (4.11) Таким образом, получили вычислительную формулу (метод Ньютона), где в качестве нулевого приближения x(0) можно взять приближенное (грубое) значение искомого корня. Пример 4.1. Рассмотрим применение метода Ньютона на примере системы двух нелинейных уравнений (4.12) Прежде чем разбирать конкретные шаги по решению системы (4.12), распишем в общем виде якобиан для системы из двух уравнений Здесь A, B, C, D – функционалы от переменных x1, x2. Нас фактически интересует W-1. Пусть матрица W- неособенная, тогда обратная матрица вычисляется Теперь вернемся к системе (4.12). Графическое решение этой системы дает две точки пересечения: М1 (1.4; -1.5) и М2 (3.4; 2.2). Зададим начальное приближение: Используя формулу (4.11), получим: Аналогично получим: 4.2. Метод градиента (метод скорейшего спуска) Пусть имеется система нелинейных уравнений: (4.13) Систему (4.13) удобнее записать в матричном виде: (4.14) где - вектор – функция; - вектор – аргумент. Решение системы (4.14), как и для системы линейных уравнений (см. п. 3.8), будем искать в виде (4.15) Здесь и - векторы неизвестных на p и p+1 шагах итераций; вектор невязок на p-ом шаге – f(p) = f(x(p)); W'p – транспонированная матрица Якоби на p – ом шаге; ; . Пример 4.2. Методом градиента вычислим приближенно корни системы расположенные в окрестности начала координат. Имеем: Выберем начальное приближение: По вышеприведенным формулам найдем первое приближение: Аналогичным образом находим следующее приближение: Ограничимся двумя итерациями (шагами), и оценим невязку: Замечания • Как видно из примера, решение достаточно быстро сходится, невязка быстро убывает. • При решении системы нелинейных уравнений методом градиента матрицу Якоби необходимо пересчитывать на каждом шаге (итерации). Вопросы для самопроверки • Как найти начальное приближение: а) для метода Ньютона; б) для метода градиента? • В методе скорейшего спуска вычисляется Якобиан (матрица Якоби). Чем отличается Якобиан, вычисленный для СЛАУ, от Якобиана, вычисленного для нелинейной системы уравнений? • Каков критерий остановки итерационного процесса при решении системы нелинейных уравнений: а) методом Ньютона; б) методом скорейшего спуска? 5. ИНТЕРПОЛИРОВАНИЕ И ПРИБЛИЖЕНИЕ ФУНКЦИЙ 5.1. Задача интерполирования и аппроксимации функций Задача интерполирования состоит в том, чтобы по значениям функции f(x) в нескольких точках отрезка восстановить ее значения в остальных точках данного отрезка. Разумеется, такая постановка задачи допускает сколь угодно много решений. Задача интерполирования возникает, например, в том случае, когда известны результаты измерений yk = f(xk) некоторой физической величины f(x) в точках xk, k = 0, 1,…, n и требуется определить ее значение в других точках. Интерполирование используется также при необходимости сгущения таблиц, когда вычисление значений f(x) по точным формулам трудоемко. Иногда возникает необходимость приближенной замены (аппроксимации) данной функции (обычно заданной таблицей) другими функциями, которые легче вычислить. При обработке эмпирических (экспериментальных) зависимостей, результаты обычно представлены в табличном или графическом виде. Задача заключается в аналитическом представлении искомой функциональной зависимости, то есть в подборе формулы, корректно описывающей экспериментальные данные. 5.2. Интерполирование алгебраическими многочленами Пусть функциональная зависимость задана таблицей y0 = f(x0);…, y1= f(x1);…,yn = f(xn). Обычно задача интерполирования формулируется так: найти многочлен P(x) = Pn(x) степени не выше n, значения которого в точках xi (i = 0, 1 2,…, n) совпадают со значениями данной функции, то есть P(xi) = yi. Геометрически это означает, что нужно найти алгебраическую кривую вида (5.1) проходящую через заданную систему точек Мi(xi,yi) (см. рис. 5.1). Многочлен Р(х) называется интерполяционным многочленом. Точки xi (i = 0, 1, 2,…, n) называются узлами интерполяции. Рис. 5.1. Интерполирование алгебраическим многочленом Для любой непрерывной функции f(x) сформулированная задача имеет единственное решение. Действительно, для отыскания коэффициентов а0, а1, а2 ,…, аn получаем систему линейных уравнений (5.2) определитель которой (определитель Вандермонда) отличен от нуля, если среди точек xi (i = 0, 1, 2,…, n) нет совпадающих. Решение системы (5.2) можно записать различным образом. Однако наиболее употребительна запись интерполяционного многочлена в форме Лагранжа и в форме Ньютона. Запишем без вывода интерполяционный многочлен Лагранжа: (5.3) Нетрудно заметить, что старшая степень аргумента х в многочлене Лагранжа равна n. Кроме этого, несложно показать, что в узловых точках значение интерполяционного многочлена Лагранжа соответствует заданным значениям f(xi). 5.3. Интерполяционная формула Ньютона Интерполяционная формула Ньютона позволяет выразить интерполяционный многочлен Pn(x) через значение f(x) в одном из узлов и через разделенные разности функции f(x), построенные по узлам x0, x1,…, xn. Эта формула является разностным аналогом формулы Тейлора: (5.4) Прежде чем приводить формулу Ньютона, рассмотрим сведения о разделенных разностях. Пусть в узлах известны значения функции f(x). Предполагаем, что среди точек xk, k = 0, 1,…, n нет совпадающих. Тогда разделенными разностями первого порядка называются отношения (5.5) Будем рассматривать разделенные разности, составленные по соседним узлам, то есть выражения . По этим разделенным разностям первого порядка можно построить разделенные разности второго порядка: (5.6) Аналогично определяются разности более высокого порядка. То есть пусть известны разделенные разности k-го порядка тогда разделенная разность k+1-го порядка определяется как (5.7) Интерполяционным многочленом Ньютона называется многочлен (5.8) Показано, что интерполяционный многочлен Лагранжа (5.3) совпадает с интерполяционным многочленом Ньютона (5.8). Замечания • В формуле (5.8) не предполагалось, что узлы x0, x1,…, xn расположены в каком-то определенном порядке. Поэтому роль точки x0 в формуле (5.8) может играть любая из точек x0, x1,…, xn. Соответствующее множество интерполяционных формул можно получить из (5.8), перенумеровав узлы. Например, тот же самый многочлен Pn(x) можно представить в виде (5.9) • Если то (5.8) называется формулой интерполирования вперед, а (5.9) - формулой интерполирования назад. • Интерполяционную формулу Ньютона удобнее применять в том случае, когда интерполируется одна и та же функция f(x), но число узлов интерполяции постепенно увеличивается. Если узлы интерполяции фиксированы и интерполируется не одна, а несколько функций, то удобнее пользоваться формулой Лагранжа. 5.4. Сходимость интерполяционного процесса Обсудим следующий вопрос: будет ли стремиться к нулю погрешность интерполирования f(x) – Ln(x), если число узлов n неограниченно увеличивать: 1. Свойства сходимости или расходимости интерполяционного процесса зависят как от выбора последовательности сеток, так и от гладкости функции f(x). 2. Известны примеры несложных функций, для которых интерполяционный процесс расходится. Так последовательность интерполяционных многочленов, построенных для непрерывной функции по равноотстоящим узлам на отрезке [-1, 1], не сходится к функции ни в одной точке отрезка [-1, 1], кроме точек –1, 0, 1. На рис. 5.2 в качестве иллюстрации изображен график многочлена L9(x) при , построенного для функции по равноотстоящим узлам на отрезке [-1,1]. Рис. 5.2. Сходимость интерполяционных многочленов 3. Чтобы избежать этих некорректностей, в практике вычислений обычно избегают пользоваться интерполяционными многочленами высокой степени. 5.5. Задача обратного интерполирования Пусть функция y = f(x) задана таблично. Задача обратного интерполирования заключается в том, чтобы по заданному значению функции y определить соответствующее значение аргумента x. Для случая неравноотстоящих значений аргумента x0, x1,…, xn задача может быть непосредственно решена с помощью интерполяционного многочлена Лагранжа. В этом случае достаточно принять переменную y за независимую и написать формулу (аналог выражения (5.3)), выражающую х как функцию у: . (5.10) Можно также, считая у аргументом, использовать формулу Ньютона: (5.11) Замечание. Обратное интерполирование корректно только для взаимно однозначных функций. Пример 5.1. Исходная функция y = f(x) задана табл. 5.1: Таблица 5.1 x 10 15 17 20 y 3 7 11 17 Необходимо найти значение функции y при x = 12; найти значение x, для которого y = 10. Решение. В качестве примера задачу прямого интерполирования в начале таблицы с неравноотстоящими узлами решим по формулам Ньютона (5.8); для обратного интерполирования применим формулу Лагранжа (5.10). y(12) = f(x0) + (x –x0)f(x0, x1) + (x –x0) (x –x1) f(x0, x1, x2) + + (x –x0) (x –x1) (x –x2) f(x0, x1, x2, x3) = 3 + 2·0.8 + + 2·(-3)·0.02857 + 2·(-3)·(-5)·(-0.002857) = 4.3429. x(10) = 10·3·(-1)·(-7)/[(-4)(-8)(-14)] + 15·7·(-1)·(-7)/[4·(-4)(-10)] + + 17·7·3(-7)/[8·4·(-6)] + 20·7·4·1/[14·10·6] = 16.641. 5.6. Отыскание параметров эмпирических формул методом наименьших квадратов При эмпирическом (экспериментальном) изучении функциональной зависимости одной величины у от другой х производят ряд измерений величины у при различных значениях величины х. Полученные результаты можно представить в виде таблицы, графика: X x1 x2 … xn Y y1 y2 … yn Задача заключается в аналитическом представлении искомой функциональной зависимости, то есть в подборе функции, описывающей результаты эксперимента. Особенность задачи состоит в том, что наличие случайных ошибок измерений делает неразумным подбор такой формулы, которая точно описывала бы все опытные значения, то есть график искомой функции не должен проходить через все экспериментальные точки. Эмпирическую формулу обычно выбирают из формул определенного типа: (5.12) Таким образом, задача сводится к определению параметров a, b, c,… формулы, в то время как вид формулы известен заранее из каких-либо теоретических соображений или из соображения простоты аналитического представления эмпирического материала. Пусть выбранная эмпирическая зависимость имеет вид (5.13) с явным указанием всех параметров, подлежащих определению. Эти параметры а0, а1, а2,…, аn нельзя определить точно по эмпирическим значениям функции y0, y1, y2,…, yk, так как последние содержат случайные ошибки. Таким образом, речь может идти только о получении достаточно хороших оценок искомых параметров. Метод наименьших квадратов (МНК) позволяет получить несмещенные и состоятельные оценки всех параметров а0, а1, а2,…, аn. 5.7. Суть метода наименьших квадратов Дальнейшие рассуждения будем проводить в предположении, что все измерения значений функции y0, y1, y2,…, yn произведены с одинаковой точностью. Тогда оценка параметров а0, а1, а2,…, аn определяется из условия минимума суммы квадратов отклонений измеренных значений yk от расчетных f(xk; а0, а1, а2,…, аn): (5.14) Отыскание же значений параметров а0, а1, а2,…, аn, которые доставляют min значение функции (5.15) сводится к решению системы уравнений (5.16) Наиболее распространен способ выбора функции f(xk; а0, а1, а2,…, аn) в виде линейной комбинации: (5.17) Здесь базисные функции (известные); n << k; а0, а1, а2,…, аn – коэффициенты, определяемые методом наименьших квадратов. Запишем в явном виде условие (5.16), учитывая выражение (5.17): (5.18) Из системы линейных уравнений (5.18) определяются все коэффициенты ak. Система (5.18) называется системой нормальных уравнений, матрица которой имеет вид (5.19) Здесь (5.20) Матрица (5.19) называется матрицей Грама. Расширенную матрицу получим добавлением справа столбца свободных членов: (5.21), где (5.22) Основные свойства матрицы Грама 1. Матрица симметрична относительно главной диагонали, то есть . 2. Матрица является положительно определенной. Следовательно, при решении методом Гаусса можно воспользоваться схемой единственного деления. 3. Определитель матрицы будет отличен от нуля, если в качестве базиса выбраны линейно независимые функции ; в этом случае система (5.18) имеет единственное решение. В качестве базисных можно выбрать линейно независимые степенные функции (5.23) Следует учесть, что n << k. Тогда для этих функций расширенная матрица Грама примет вид (5.24) Если выбрать n = k, то на основании единственности интерполяционного полинома получим функцию , совпадающую с каноническим интерполяционным полиномом степени k. При этом аппроксимирующая кривая пройдет через все экспериментальные точки, и функция S будет равна нулю. Пример 5.2. Исходная функция y = f(x) задана в виде табл. 5.2: Таблица 5.2 x 10 15 17 20 y 3 7 11 17 Аппроксимируем экспериментальные данные линейной либо квадратичной функцией. Методом наименьших квадратов необходимо уточнить коэффициенты аппроксимирующего полинома. Решение 1. При линейной аппроксимации исходную зависимость представим в виде , где . Методом наименьших квадратов определим a0 и a1. Расширенная матрица Грама в нашем случае имеет вид  ; а1 = 1.3774; а0 =-11.8491. Таким образом, аппроксимирующая функция равна Оценим погрешность формулы, и результаты этой оценки сведем в табл. 5.3: Таблица 5.3 x y f y - f |y-f| / |y| 10 3 1.9249 1.0751 0.3584 15 7 8.8119 -1.8119 0.2588 17 11 11.5667 -0.5667 0.0515 20 17 15.6989 1.3011 0.07654 Для нашей линейной функции S1 = 6.4528. 2. Решим ту же задачу, аппроксимировав эмпирические данные полиномом второй степени: . Матрица Грама в этом случае имеет вид Все результаты сведены в табл. 5.4. Таблица 5.4 x y f y - f |y-f| / |y| 10 3 2.9511 0.0489 0.0163 15 7 7.3381 -0.3381 0.0483 17 11 10.6007 0.3993 0.0363 20 17 17.1101 -0.1101 0.0065 S2 = 0.2883. Обсуждение результатов 1. Аппроксимировав эмпирические результаты более простой функцией (линейной), мы получили погрешность в различных узловых точках, лежащую в пределах от 5 до 35 %. 2. Более сложная формула квадратичной интерполяции обеспечивает погрешность не более 5 %. 3. Косвенную оценку погрешности можно провести, сравнив значения S1 и S2. 4. Матрица Грама для полинома второй, третьей степени имеет простой вид и может быть решена, например, методом Гаусса. Вопросы для самопроверки • В чем заключается задача интерполирования и аппроксимации? • Запишите интерполяционные формулы Лагранжа и Ньютона. • Какие требования предъявляются а) к интерполяционным полиномам; б) к аппроксимационным полиномам? • Что такое разделенные разности? • В каких случаях применяются формулы Ньютона для интерполирования а) вперед, б) назад? • Что можно сказать о сходимости интерполяционных полиномов? • Что такое обратное интерполирование, при каких условиях оно возможно (корректно)? • В чем заключается идея метода наименьших квадратов? • Что такое матрица Грама, каковы ее свойства? Что такое базисные функции? Можно ли в качестве базисных функций выбрать а) линейно независимые функции; б) линейно зависимые функции? 6. ПРИБЛИЖЕННОЕ ДИФФЕРЕНЦИРОВАНИЕ И ИНТЕГРИРОВАНИЕ ФУНКЦИИ ОДНОЙ ПЕРЕМЕННОЙ 6.1. Численное интегрирование 6.1.1 Постановка задачи Задача численного интегрирования функции заключается в вычислении определенного интеграла на основании ряда значений подынтегральной функции. Численное вычисление однократного интеграла называется механической квадратурой. Мы будем рассматривать способы приближенного вычисления определенных интегралов , (6.1.1) основанные на замене интеграла конечной суммой: , (6.1.2) где Сk- числовые коэффициенты, а xk  [a, b], k = 0, 1, …, n. Приближенное равенство (6.1.3) называется квадратурной формулой, а xk – узлами квадратурной формулы. Погрешность квадратурной формулы определяется соотношением . (6.1.4) В общем случае погрешность квадратурной формулы (6.1.4) зависит как от выбора коэффициентов Ск , так и от расположения узлов хк. Введем на отрезке [a, b] равномерную сетку с шагом h, тогда xi = a + ih, где (i = 0, 1, ..., n; h·n = b-a). Теперь выражение (6.1.1) можно представить в виде суммы интегралов по частичным отрезкам: (6.1.5) Таким образом, для построения формулы численного интегрирования на отрезке [a, b] достаточно построить квадратурную формулу на частичном отрезке [xi-1, xi] и воспользоваться формулой (6.1.5). 6.1.2. Формула прямоугольников На частичном отрезке [xi-1, xi] заменим подынтегральную функцию полиномом Лагранжа нулевого порядка, построенным в одной точке. Естественно в качестве этой точки выбрать среднюю: xi-0.5 = xi - 0.5h. Тогда получим формулу . (6.1.6) Подставив (6.1.6) в (6.1.5), получим составную формулу средних прямоугольников: . (6.1.7) Графическая иллюстрация метода средних прямоугольников представлена на рис. 6.1. Рис. 6.1. Интегрирование методом средних прямоугольников Погрешность формулы (6.1.7) определяется выражением (6.1.8) Здесь . Таким образом, погрешность формулы (6.1.7) пропорциональна O(h2). Замечание. Формулу (6.1.7) можно представить в ином виде: . (6.1.9) Эти формулы в выражении (6.1.9) называются формулой левых и правых прямоугольников соответственно. Графически метод левых и правых прямоугольников представлен на рис. 6.2. а) б) Рис. 6.2. Метод левых (а) и правых (б) прямоугольников Однако из-за нарушения симметрии в формулах (6.1.9) их погрешность значительно больше, чем в методе средних прямоугольников и ~O(h). 6.1.3. Формула трапеций Если на частичном отрезке подынтегральную функцию заменить полиномом Лагранжа первой степени, то есть , (6.1.10) тогда искомый интеграл запишется следующим образом: (6.1.11) После подстановки выражения (6.1.11) в (6.1.5) составная формула трапеций примет вид (6.1.12) Графически метод трапеций представлен на рис. 6.3. Рис. 6.3. Метод трапеций Погрешность формулы (6.1.12) определяется выражением: (6.1.13) Таким образом, погрешность метода трапеций Ψ ~ O(h²), но она в два раза больше, чем для формулы средних прямоугольников. 6.1.4. Формула Симпсона В этом методе предлагается подынтегральную функцию на частичном отрезке аппроксимировать параболой, проходящей через точки (xj, f(xj)), где j = i-1; i-0.5; i, то есть подынтегральную функцию аппроксимируем интерполяционным многочленом Лагранжа второй степени: (6.1.14) Проведя интегрирование, получим: (6.1.15) Это и есть формула Симпсона или формула парабол. На отрезке [a, b] формула Симпсона примет вид (6.1.16) Графическое представление метода Симпсона показано на рис. 6.4. Рис. 6.4. Метод Симпсона Избавимся в выражении (6.1.16) от дробных индексов, переобозначив переменные: (6.1.17) Тогда формула Симпсона примет вид (6.1.18) Погрешность формулы (6.1.18) оценивается следующим выражением: , (6.1.19) где h·n = b - a, . Таким образом, погрешность формулы Симпсона пропорциональна O(h4). Замечание. Следует отметить, что в формуле Симпсона отрезок интегрирования обязательно разбивается на четное число интервалов. 6.1.5. Вычисление определенных интегралов методами Монте–Карло Рассматриваемые ранее методы называются детерминированными, то есть лишенными элемента случайности. Методы Монте–Карло (ММК) – это численные методы решения математических задач с помощью моделирования случайных величин. ММК позволяют успешно решать математические задачи, обусловленные вероятностными процессами. Более того, при решении задач, не связанных с какими-либо вероятностями, можно искусственно придумать вероятностную модель (и даже не одну), позволяющую решать эти задачи. Рассмотрим вычисление определенного интеграла (6.1.20) При вычислении этого интеграла по формуле прямоугольников интервал [a, b] разбиваем на N одинаковых интервалов, в серединах которых вычислялись значения подынтегральной функции. Вычисляя значения функции в случайных узлах, можно получить более точный результат: (6.1.21) (6.1.22) Здесь γi - случайное число, равномерно распределенное на интервале [0, 1]. Погрешность вычисления интеграла ММК ~ , что значительно больше, чем у ранее изученных детерминированных методов. На рис. 6.5 представлена графическая реализация метода Монте-Карло вычисления однократного интеграла со случайными узлами (6.1.21) и (6.1.22). Рис. 6.5. Интегрирование методом Монте-Карло (1-й случай) Однако при вычислении кратных интегралов детерминированными методами оценка погрешности перерастает в задачу порой более сложную, чем вычисление интеграла. В то же время погрешность вычисления кратных интегралов ММК слабо зависит от кратности и легко вычисляется в каждом конкретном случае практически без дополнительных затрат. Рассмотрим еще один метод Монте-Карло на примере вычисления однократного интеграла: (6.1.23) Рис. 6.6. Интегрирование методом Монте-Карло (2-й случай) Как видно на рис. 6.6, интегральная кривая лежит в единичном квадрате, и если мы сумеем получать пары случайных чисел, равномерно распределенных на интервале [0, 1], то полученные значения (γ1, γ2) можно интерпретировать как координаты точки в единичном квадрате. Тогда, если этих пар чисел получено достаточно много, можно приблизительно считать, что . Здесь S – число пар точек, попавших под кривую, а N – общее число пар чисел. Пример 6.1. Вычислить следующий интеграл: Поставленная задача была решена различными методами. Полученные результаты сведены в табл. 6.1. Таблица 6.1 Число интервалов (точек) Метод левых прямоугольников Метод средних прямоугольников Метод правых прямоугольников Метод трапеций Метод Симпсона Метод Монте-Карло 10 4.44112722 4.66882868 4.90820465 4.25683746 4.67077443 4.62289422 100 4.64745932 4.67075481 4.69416706 4.62903035 4.67077427 4.69812790 Замечание. Выбор табличного интеграла позволил нам сравнить погрешность каждого метода и выяснить влияние числа разбиений на точность вычислений. Вопросы для самопроверки • Сформулируйте задачу численного интегрирования. • Метод средних, левых и правых прямоугольников. Что можно сказать об их погрешности, трудоемкости? • Задача численного интегрирования решена методом трапеций. Предложите и обоснуйте пути повышения точности (уменьшения погрешности) расчетов. • Сравните метод трапеций и метод Симпсона. • Какие методы Монте–Карло численного интегрирования вы знаете? Сравните эти методы с любым детерминированным. • Необходимо вычислить интеграл методами трапеций, Симпсона и ММК, разбив область интегрирования на 77 интервалов (точек). Что можно сказать о точности и применимости этих методов? 6.2. Решение обыкновенных дифференциальных уравнений 6.2.1. Методы решения задачи Коши Среди задач, с которыми приходится иметь дело в вычислительной практике, значительную часть составляют различные задачи, сводящиеся к решению обыкновенных дифференциальных уравнений. Обычно приходится прибегать к помощи приближенных методов решения подобных задач. В случае обыкновенных дифференциальных уравнений в зависимости от того, ставятся ли дополнительные условия в одной или нескольких точках отрезка изменения независимой переменной, задачи обычно подразделяются на одноточечные (задачи с начальными условиями или задачи Коши) и многоточечные. Среди многоточечных задач наиболее часто в прикладных вопросах встречаются так называемые граничные задачи, когда дополнительные условия ставятся на концах рассматриваемого отрезка. В дальнейшем ограничимся рассмотрением численных методов решения задачи Коши. Для простоты изложения методов решения задачи будем рассматривать случай одного обыкновенного дифференциального уравнения первого порядка. Пусть на отрезке x0  x  b требуется найти решение y(x) дифференциального уравнения , (6.2.1) удовлетворяющее при x = x0 начальному условию (6.2.2) Будем считать, что условия существования и единственности решения поставленной задачи Коши выполнены. На практике найти общее либо частное решение задачи Коши удается крайне редко, поэтому приходится решать эту задачу приближенно. Отрезок [x0, b] накрывается сеткой (разбивается на интервалы) чаще всего с постоянным шагом h ( h = xn+1 - xn ), и по какому-то решающему правилу находится значение yn+1 = y(xn+1). Таким образом, в качестве решения задачи Коши численными методами мы получаем таблицу, состоящую из двух векторов: x = (x0 , x1 , …xn) – вектора аргументов и соответствующего ему вектора функции y = ( y0 , y1,… yn ). Численные методы (правила), в которых для нахождения значения функции в новой точке используется информация только об одной (предыдущей) точке, называются одношаговыми. Численные методы (правила), в которых для нахождения значения функции в новой точке используется информация о нескольких (предыдущих) точках, называются многошаговыми. Из общего курса обыкновенных дифференциальных уравнений широкое распространение получил аналитический метод, основанный на идее разложения в ряд решения рассматриваемой задачи Коши. Особенно часто для этих целей используется ряд Тейлора. В этом случае вычислительные правила строятся особенно просто. При этом приближенное решение ym(x) исходной задачи ищут в виде (6.2.3) Здесь а значения i = 2, 3,…m находят по формулам, полученным последовательным дифференцированием уравнения (6.2.1): (6.2.4) Для значений x, близких к x0, метод рядов (6.2.3) при достаточно большом значении m дает обычно хорошее приближение к точному решению y(x) задачи (6.2.1). Однако с ростом расстояния x - x0  погрешность приближенного равенства y(x)  ym(x), вообще говоря, возрастает по абсолютной величине, и правило (6.2.3) становится вовсе неприемлемым, когда x выходит из области сходимости соответствующего ряда (6.2.3) Тейлора. Если в выражении (6.2.3) ограничиться m = 1, то для вычисления новых значений y(x) нет необходимости пересчитывать значение производной, правда и точность решения будет невысока. Графическая интерпретация этого метода приведена на рис. 6.2.1. Рис. 6.2.1. Разложение функции в ряд Тейлора (m=1) 6.2.2. Метод рядов, не требующий вычисления производных правой части уравнения Естественно поставить задачу о таком усовершенствовании приведенного выше одношагового метода, которое сохраняло бы основные его достоинства, но не было бы связано с нахождением значений производных правой части уравнения (6.2.5) где xn+1 = xn + h. Чтобы выполнить это условие (последнее), производные y(i)(x), i = 2, 3,..., m, входящие в правую часть уравнения (6.2.5), можно заменить по формулам численного дифференцирования их приближенными выражениями через значение функции y' и учесть, что y'(x) = f [x, y(x)]. В случае m = 1 приближенное равенство (6.2.5) не требует вычисления производных правой части уравнения и позволяет с погрешностью порядка h2 находить значение y(xn+ h) решения этого уравнения по известному его значению y(xn). Соответствующее одношаговое правило можно записать в виде (6.2.6) Это правило (6.2.6) впервые было построено Эйлером и носит его имя. Иногда его называют также правилом ломаных или методом касательных. Метод Эйлера имеет простую геометрическую интерпретацию (см. рис. 6.2.2). Рис. 6.2.2. Нахождение решения методом Эйлера Замечание Метод Эйлера имеет порядок точности ~ h2 на одном шаге. Практическая оценка погрешности приближенного решения может быть получена по правилу Рунге. 6.2.3. Метод Рунге-Кутта Изложим идею метода на примере задачи Коши: (6.2.7) Интегрируя это уравнение в пределах от x до x + h (0 < h <1), получим равенство (6.2.8) которое посредством последнего интеграла связывает значения решения рассматриваемого уравнения в двух точках, удаленных друг от друга на расстояние шага h. Для удобства записи выражения (6.2.8) используем обозначение ∆y = y(x + h) – y(x) и замену переменной интегрирования t = x + h. Окончательно получим: (6.2.9) Указав эффективный метод приближенного вычисления интеграла в выражении (6.2.9), мы получим при этом одно из правил численного интегрирования уравнения (6.2.7). Постараемся составить линейную комбинацию величин i, i = 0, 1, ..., q, которая будет являться аналогом квадратурной суммы и позволит вычислить приближенное значение приращения y: (6.2.10) где Метод четвертого порядка для q = 3, являющийся аналогом широко известной в литературе четырехточечной квадратурной формулы "трех восьмых", имеет вид (6.2.11) где Особо широко известно другое вычислительное правило типа Рунге-Кутта четвертого порядка точности: (6.2.12) где Метод Рунге-Кутта имеет погрешность четвертого порядка (~ h4 ). Правило Рунге. Если приближенный метод имеет порядок погрешности m, то погрешность можно приближенно оценить по формуле (6.2.13) В формуле (6.2.13) O(xi) – главный член погрешности, и - приближенные решения в точке xi, найденные с шагом h и 2h соответственно. Пример 6.2.1. Решить дифференциальное уравнение на отрезке [0, 1] c начальным условием y(x=0) = 1. Найти первые три точки, приняв шаг h = 0.05. Решение. Поставленная задача была решена методом разложения в ряд Тейлора (6.2.3); методом Эйлера (6.2.6) и методом Рунге-Кутта (6.2.12). Для наглядности все полученные результаты сведем в табл. 6.2.1. Таблица 6.2.1 xi Ряд Тейлора (m=1) Метод Эйлера Метод Рунге-Кутта yi yi yi f(xi, yi) φ0 φ1 φ2 φ3 1 1 1 1 - - - - 0.05 1.05 1.05 1.0477 0.9089 0.05 0.0477 0.0476 0.0454 0.1 1.1 1.0931 1.0912 0.8321 0.0454 0.0435 0.0434 0.0416 0.15 1.15 1.1347 1.1311 0.7658 0.0416 0.0399 0.0399 0.0383 6.2.4. Многошаговые методы Ранее нами были рассмотрены одношаговые методы решения задачи Коши. Эти методы, обладая рядом удобных для практики вычислений особенностей, страдают одним существенным недостатком. При построении этих методов привлекается информация о решаемой задаче только на отрезке длиной в один шаг, поэтому подобная информация на каждом этапе процесса должна быть получена заново, что предопределяет большую трудоемкость соответствующих вычислительных правил. Если отказаться от требования одношаговости, можно вычислительные методы строить таким образом, чтобы часть получаемой информации могла быть использована повторно на нескольких следующих шагах вычислительного процесса. Такие методы, использующие информацию о решаемой задаче на отрезке длиной более одного шага, и называются многошаговыми. Будем, как и раньше рассматривать задачу Коши: (6.2.14) Ограничимся рассмотрением многошаговых методов с равномерной сеткой: xi = x0 + ih; i = 0, 1,..., n; n·h = b - x0. (6.2.15) Рассмотрим вычислительные правила вида (6.2.16) Среди вычислительных правил вида (6.2.16) особенно широко известны экстраполяционные (при s = 0) и интерполяционные (при s = 1, A-1  0). 6.2.5. Экстраполяционные методы Адамса Экстраполяционные формулы Адамса получаются из (6.2.16) при s = 0. Если же предположим при этом, что q = 0, то получим уже знакомый нам метод Эйлера: (6.2.17) При q = 3 из (6.2.16) получим следующий вид формулы Адамса: (6.2.18) Здесь приняты следующие обозначения: (6.2.19) Рекуррентная формула для определения конечных разностей j – го порядка имеет вид (6.2.20) Учитывая (6.2.20), получим: (6.2.21) 6.2.6. Интерполяционные методы Адамса При s = 1 формула (6.2.16) примет вид (6.2.22) Если q = 2, получим следующее вычислительное правило: (6.2.23) Обычно на практике используют экстраполяционную формулу (6.2.18), а затем корректируют полученное значение по формуле (6.2.23). И если результат уточненного значения не превышает допустимую погрешность расчета, то шаг h считается допустимым . Для расчетов на компьютере формулы (6.2.18) и (6.2.23) в конечно-разностном виде неудобны. С учетом (6.2.21) их можно представить в виде (6.2.24) Приведенные формулы имеют достаточно большую точность. Они дают погрешность порядка ~ О( h4 ), но сами формулы оценки погрешности достаточно сложны. Приближенно погрешность можно оценить по правилу Рунге. Пример 6.2.2. Решить дифференциальное уравнение на отрезке [0, 1] c начальным условием y(x=0) = 1. Найти решение методом Адамса (с коррекцией) в точке x4, решение в трех первых точках найти методом Рунге- Кутта, приняв шаг . Решение. Значения функции в четырех первых точках возьмем из табл. 6.2.1 (см. пример в предыдущем разделе). Теперь стало понятно, зачем мы сохраняли значения первой производной в этих точках (см. формулы (6.2.24)). x4 = x3 + h = 0.15 + 0.05 = 0.2; Для того чтобы скорректировать полученный результат, необходимо вычислить значение производной в этой точке: Теперь уточним значение по интерполяционной формуле (а можно этого и не делать, тогда погрешность метода будет больше): Так как в качестве нового значения функции принято скорректированное, то обязательно следует пересчитать значение производной. В нашем случае модуль разности экстраполяционной и интерполяционной формул меньше ε, что позволяет продолжить вычисления с тем же шагом. Вопросы для самопроверки • Сформулируйте задачу Коши для обыкновенных дифференциальных уравнений первого порядка. • Что является решением дифференциального уравнения: а) в высшей математике, б) в прикладной математике? • Какие методы решения дифференциальных уравнений называются одношаговыми, многошаговыми? Приведите примеры. • Сравните решения, полученные на первом, втором шаге методами Эйлера, Рунге-Кутта и разложением в ряд Тейлора (трудоемкость, погрешность…). • Как оценить погрешность применяемого метода? Как ее уменьшить? • Сравните одношаговые и многошаговые методы решения дифференциальных уравнений, указав достоинства и недостатки первых и вторых. • Что такое экстраполяционные и интерполяционные методы (формулы) Адамса? • Можно ли применять: а) только экстраполяционные методы Адамса, б) только интерполяционные? • Можно ли использовать: а) многошаговые методы без одношаговых; б) одношаговые методы без многошаговых? • При решении дифференциального уравнения методом Адамса на 27-м шаге необходимо сменить шаг. Как это сделать? 7. СТАТИСТИЧЕСКИЕ И ДИНАМИЧЕСКИЕ МОДЕЛИ. МЕТОД МОНТЕ КАРЛА • 7.1 Статистическое моделирование как метод прогнозирования долговечности технических систем На этапе проектирования технической системы следует осуществить оптимизацию структуры и подобрать оптимальные значения параметров, учитывая требования обеспечения надежности. Их подбор можно осуществить путем статистического моделирования. с этой целью следует построить статистическую модель системы, которая как можно лучше отображала бы ее структуру надежности, а также характер функционирования элементов при принятых внешних условиях. В результате проведенного статистического моделирования определяются величины характеристик (показателей) надежности отдельных конструкционных модификаций системы. решения, при которых прогнозируемые характеристики надежности обеспечивают ранее принятые требования, следует считать соответствующими (принятыми) решениями. Статистическое моделирование параметров долговечности и надежности базируется на методе монте карло. в этом методе находят применение генераторы случайных чисел с целью предусмотрения (прогнозирования) повреждений в процессе изнашивания в соответствии с принципами. В наиболее общем случае процедура статистического моделирования состоит в том, что определенные и генерированные ряды цифр, моделирующие отдельные случайные переменные х, вводятся в математическую модель функции t = f(xi), которая будет описывать долговечность триботехнической системы. Для выбранной модели t = f(xi) проводятся вычисления переменной t, а в дальнейшем осуществляется статистическая обработка результатов моделирования. сначала переменная группируется в параметрические группы (классы) по интервалам ее изменения – вариационный ряд, гистограмма. В дальнейшем вычисляются их дискретные статистические характеристики (среднее арифметическое, среднее квадратичное, коэффициент вариации, коэффициенты асимметрии и т.д.). Наконец, проводится определение функций отказа и надежности. Результаты прогнозируемого распределения долговечности Т при статистическом моделировании в дальнейшем проверяются (тестируются) на соответствие теоретическому распределению. 7.2 Метод Монте-Карло Специальный метод изучения поведения заданной статистики при проведении многократных повторных выборок, существенно использующий вычислительные возможности современных компьютеров. При проведении анализа по методу Монте-Карло компьютер использует процедуру генерации псевдослучайных чисел для имитации данных из изучаемой генеральной совокупности. Процедура анализа по методу Монте-Карло модуля Моделирование структурными уравнениями строит выборки из генеральной совокупности в соответствии с указаниями пользователя, а затем производит следующие действия: Для каждого повторения по методу Монте-Карло: 1. Имитирует случайную выборку из генеральной совокупности, 2. Проводит анализ выборки, 3. Сохраняет результаты. После большого числа повторений, сохраненные результаты хорошо имитирует реальное распределение выборочной статистики. Метод Монте-Карло позволяет получить информацию о выборочном распределении в случаях, когда обычная теория выборочных распределений оказывается бессильной. ЭВМ позволяют легко получать так называемые псевдослучайные числа (при решении задач их применяют вместо случайных чисел); это привело к широкому внедрению метода во многие области науки и техники (статистическая физика, теория массового обслуживания, теория игр и др.). Метод Монте-Карло используют для вычисления интегралов, в особенности многомерных, для решения систем алгебраических уравнений высокого порядка, для исследования различного рода сложных систем (автоматического управления, экономических, биологических и т.д.). Сущность метода Монте-Карло состоит в следующем: требуется найти значение а некоторой изучаемой величины. Для этого выбирают такую случайную величину X, математическое ожидание которой а: (1) Практически же поступают так: производят п испытаний; в результате которых получают п возможных значений X, вычисляют их среднее арифметическое (2) и принимают х в качестве оценки (приближенного значения) а * искомого числа а: (3) Поскольку метод Монте-Карло требует проведения большого числа испытаний, его часто называют методом статистических испытаний. Теория этого метода указывает, как наиболее целесообразно выбрать случайную величину X, как найти ее возможные значения. В частности, разрабатываются способы уменьшения дисперсии используемых случайных величин, в результате чего уменьшается ошибка, допускаемая при замене искомого математического ожидания а его оценкой а *. Отыскание возможных значений случайной величины Х (моделирование) называют «разыгрыванием случайной величины». Изложим лишь некоторые способы разыгрывания случайных величин и укажем, как оценить допускаемую при этом ошибку. Оценка погрешности метода Монте-Карло Пусть для получения оценки а * математического ожидания а случайной величины Х было произведено п независимых испытаний (разыграно п возможных значений Х) и по ним была найдена выборочная средняя , которая принята в качестве искомой оценки: а* = . Ясно, что если повторить опыт, то будут получены другие возможные значения X, следовательно, другая средняя, а значит, и другая оценка а*. Уже отсюда следует, что получить точную оценку математического ожидания невозможно. Естественно, возникает вопрос о величине допускаемой ошибки. Ограничимся отысканием лишь верхней границы допускаемой ошибки с заданной вероятностью (надежностью) : (4) Интересующая нас верхняя граница ошибки есть не что иное, как «точность оценки» математического ожидания по выборочной средней при помощи доверительных интервалов. Поэтому воспользуемся результатами, полученными ранее, и рассмотрим следующие три случая. 1. Случайная величина Х распределена нормально и ее среднее квадратическое отклонение - известно. В этом случае с надежностью у верхняя граница ошибки (5) где п — число испытаний (разыгранных значений X); t — значение аргумента функции Лапласа, при котором Ф(t) == /2, — известное среднее квадратическое отклонение X. 2. Случайная величина Х распределена нормально, причем ее среднее квадратическое отклонение неизвестно. В этом случае с надежностью верхняя граница ошибки (6) где п — число испытаний; s — «исправленное» среднее квадратическое отклонение, t находят по таблице значений ty == t{,n}. 3. Случайная величина Х распределена по закону, отличному от нормального. В этом случае при достаточно большом числе испытаний (n > 30) с надежностью, приближенно равной , верхняя граница ошибки может быть вычислена по формуле (5), если среднее квадратическое отклонение случайной величины Х известно; если же -неизвестно, то можно подставить в формулу (5) его оценку s — «исправленное» среднее квадратическое отклонение либо воспользоваться формулой (6). Заметим, что чем больше п, тем меньше различие между результатами, которые дают обе формулы. Это объясняется тем, что при п —> распределение Стьюдента стремится к нормальному. В частности, при п=--100, =0,95 верхняя граница ошибки равна 0,098 по формуле (5) и 0,099 по формуле (6). Как видим, результаты различаются незначительно. Замечание. Для того чтобы найти наименьшее число испытаний, которые обеспечат наперед заданную верхнюю границу ошибки , надо выразить п из формул (5) и (6): Литература 1. Демидович Б.П., Марон И.А. Основы вычислительной математики. - М.: Наука, 1970. – 664 с. 2. Копченова Н.В., Марон И.А. Вычислительная математика в примерах и задачах. - М.: Наука, 1972. – 308 с. 3. Самарский А.А., Гулин А.В. Численные методы: Учеб. пособие для вузов. – М.: Наука, 1989. – 432 с. 4. Мудров А.Е. Численные методы для ПЭВМ на языках Бейсик, Фортран и Паскаль. – Томск: МП "РАСКО", 1991. – 272 с. 5. Практикум по численным методам. / Л.Я. Егорова, Л.Л. Левин, Б.Г. Ослин и др. - Томск: Изд. ТГУ, 1979. – 212 с. 6. Press W.H., Teukolsky S.A., Vetterling W.T., Flannery B.P., Numerical Recipes in C. The Art of Scientific Computing. 2-nd ed. - Copyright © Cambridge University Press, 1992, - 966p.