1) Создание модели, фиксирующей главные исследуемые факторы. Одновременно формулируются рамки применимости модели.

2) Предварительное исследование математической модели: поверка корректности постановки задачи, существования и единственности решения.

3) Разработка метода расчета сформулированной задачи, построение эффективных вычислительных алгоритмов.

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

5) Сравнение полученных результатов моделирования с тестовыми примерами и экспериментальными данными; решение вопроса о правильности практического моделирования (иначе повторяются пункты 3 и 4).

6) Решение вопроса о достоверности предложенной математической модели. Если модель не описывает экспериментальные данные, возврат на пункт 1.

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

Задача называется корректно поставленной, если для любых входных данных х из некоторого класса решение y существует, единственно и устойчиво по входным данным.

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

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

Отсутствие устойчивости обычно означает, что сравнительно небольшой погрешности дx соответствует весьма большое дy, а значит получаемое решение будет далеко от истинного. К такой задаче численные методы применять бессмысленно, ибо погрешности численного расчета будут катастрофически нарастать. Устойчивость задачи определяется (1) математической формулировкой, (2) используемым алгоритмом расчета. Пример неустойчивой задачи в первом случае:

Система имеет решение , однако

Система имеет решение ,

то есть разница в коэффициенте менее 1% приводит к изменению решения в 300%.

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

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

Если выполняется условие для норм (модулей), то задача устойчива. Однако если константа С очень велика, то фактически наблюдается слабая устойчивость. Такую задачу называют плохо обусловленной. Пример - дифференциальное уравнение с начальными условиями . Общее решение дифференциального уравнения есть . Начальные условия приводят к обнулению первого слагаемого, но если из-за погрешности начальных данных это будет не так, то при возрастании x влияние первого слагаемого будет катастрофически нарастать.

2. Точность вычислений, классификация погрешностей

Существуют четыре источника погрешности результата:

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

. То есть:

,

где - априорная погрешность метода на интервале интегрирования,

а r(x) - априорная погрешность метода на отдельном шаге интегрирования.

Обзор методов интегрирования.

Методы вычисления однократных интегралов называются квадратурными (для кратных интегралов - кубатурными).

1) Методы Ньютона-Котеса. Здесь ц(x) - полином различных степеней. Сюда относятся метод прямоугольников, трапеций, Симпсона.

2) Методы статистических испытаний (методы Монте-Карло). Здесь узлы сетки для квадратурного или кубатурного интегрирования выбираются с помощью датчика случайных чисел, ответ носит вероятностный характер. В основном применяются для вычисления кратных интегралов.

3) Сплайновые методы. Здесь ц(x) - кусочный полином с условиями связи между отдельными полиномами посредством системы коэффициентов.

4) Методы наивысшей алгебраической точности. Обеспечивают оптимальную расстановку узлов сетки интегрирования и выбор весовых коэффициентов с(x) в задаче. Сюда относится метод Гаусса-Кристоффеля (вычисление несобственных интегралов) и метод Маркова.

Метод прямоугольников.

Мы записывали априорные оценки главного члена погрешности в виде R₀ = Ah^p, (1) где A - коэффициент, зависящий от метода интегрирования и вида подинтегральной функции; h - шаг интегрирования, p - порядок метода. Такого сорта оценку можно применить не только к методам интегрирования, но и ко многим другим численным алгоритмам.

Первая формула Рунге.

Пусть w - точное значение, к которому должен прийти численный метод (мы его не знаем). Результат численного расчета дает нам величину w_h такую, что . (2)

Теперь вычислим ту же величину w с шагом kh, где константа k может быть как больше, так и меньше единицы. Коэффициент A будет одинаковый, так как вычисление осуществляется одним и тем же методом. Получаем . (3)

Приравняем правые части выражений (2) и (3) и пренебрежем бесконечно малыми величинами одинакового порядка малости.

.

Отсюда, учитывая (1), получим . (4) Эта формула, выражающая апостериорную оценку главного члена погрешности величины w путем двойного просчета с разным шагом, носит название первой формулы Рунге. При уменьшении шага главный член погрешности будет стремиться к полной погрешности R.

Вторая формула Рунге.

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

Численное дифференцирование - некорректная задача, так как отсутствует устойчивость решения. При численном дифференцировании приходится вычитать друг из друга близкие значения функции. Это приводит к уничтожению первых значащих цифр, т.е. к потере части достоверных знаков числа. А так как значения функции обычно известны с определенной погрешностью, то все значащие цифры могут быть потеряны. На графике кривая (1) соответствует уменьшению погрешности дифференцирования при уменьшении шага; кривая (2) представляет собой неограниченно возрастающий (осциллирующий) вклад неустранимой погрешности исходных данных - значений функции y(x). Критерий выхода за оптимальный шаг при его уменьшении - «разболтка» решения: зависимость результатов вычислений становится нерегулярно зависящей от величины шага.

Пусть введена как интерполяционный многочлен Ньютона. В этом случае для произвольной неравномерной сетки:

, для i = 0,1…n-1, интерполяция полиномом первой степени.

, интерполяция полиномом второй степени.

В общем случае

.

Минимальное число узлов, необходимое для вычисления k-й производной, равно k + 1.

Оценка погрешности при численном дифференцировании может быть осуществлена по формуле

,

где n - число узлов функции, k - порядок производной.

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

, где h = x₁ - x₀ = const.

.

Исходя из общего вида интерполяционного полинома можно вывести формулы для более высокого порядка точности или для более высоких производных.

7. Численное решение систем линейных уравнений

1) Прямой ход.

Идея метода состоит в последовательном исключении неизвестных из системы n линейных уравнений. На примере первого уравнения системы (2) рассмотрим выражение для x₁:

.

Подставим выражение для x₁ во второе и все остальные уравнения системы:

.

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

Общие формулы прямого хода:

, (3)

k = 1…n, j = 1…n+1. Звездочкой отмечены элементы k-й строки с измененными значениями, которые будут подставлены в следующую формулу. Для определенности будем считать первый индекс - по строкам, второй - по столбцам.

, (4)

i = k +1…n, j = 1…n+1, k фиксировано в уравнении (3). Для уменьшения количества действий достаточно изменять значения элементов, находящихся выше главной диагонали.

2) Обратный ход.

Второй этап решения СЛАУ методом Гаусса называется обратным ходом и состоит в последовательном определении x_k, начиная с x_n, так как для последнего решение фактически получено. Общая формула:

. (5)

Таким образом, вычисление корней происходит за 2/3 n³ арифметических действий.

3) Выбор главного элемента.

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

4) Погрешность метода. Расчет невязок.

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

, где k = 1…n. (6)

Специально отметим, что подставлять найденные значения следует в исходную (не преобразованную к верхнетреугольному виду) систему.

Преимущества и недостатки метода.

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

Блок-схема алгоритма метода Гаусса без выбора главного элемента.

1

Итерационные методы решения систем линейных уравнений.

Простейшим итерационным методом решения СЛАУ является метод простой итерации. При этом система уравнений (1) преобразуется к виду (2), а ее решение находится как предел последовательности (3), где {n} - номер итерации. Утверждается, что всякая система (2), эквивалентная (1), записывается в виде .

Теорема о достаточном условии сходимости метода простой итерации утверждает, что если норма матрицы (), то система уравнений (2) имеет единственное решение и итерационный процесс (3) сходится к решению со скоростью геометрической прогрессии.

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

На практике для обеспечения сходимости итерационных методов необходимо, чтобы значения диагональных элементов матрицы СЛАУ были преобладающими по абсолютной величине по сравнению с другими элементами.

Представим СЛАУ в следующей форме, удовлетворяющей (3):

(4)

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

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

Критерий сходимости метода Зейделя: Пусть - вещественная симметричная положительно определенная матрица. Тогда метод Зейделя сходится.

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

8)Интерполирование функций.

1) Последовательная интерполяция на прямоугольной сетке. Пусть заданы z _{i j} = z(x_i, y_j) требуется найти z(x, y). Сначала при фиксированных y_j0 найдем значение z(x, y_j0),

Затем по полученному набору значений найдем z(x, y).

В случае интерполяции полиномом Лагранжа общая формула имеет вид

где k и m - количество узлов по сторонам прямоугольной сетки.

2) Треугольная конфигурация узлов.

z (x₀, x₁, y) = [z(x₀, y)-z(x₁, y)]/(x₀-x₁)

z (x, y₀, y₁) = [z(x, y₀)-z(x,y₁)]/(y₀-y₁)

Многочлен Лагранжева типа в этом случае имеет вид

9. Аппроксимация функций методом наименьших квадратов

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

Введем непрерывную функцию ц(x) для аппроксимации дискретной зависимости f(x_i), i = 0…n. Будем считать, что ц(x) построена по условию наилучшего квадратичного приближения, если

. (1)

Весу с для i-й точки придают смысл точности измерения данного значения: чем больше с, тем ближе аппроксимирующая кривая «притягивается» к данной точке. В дальнейшем будем по умолчанию полагать с = 1 для всех точек.

Рассмотрим случай линейной аппроксимации:

ц(x) = c₀ц₀(x) + c₁ц₁(x) + … + c_mц_m(x), (2)

где ц₀…ц_m - произвольные базисные функции, c₀…c_m - неизвестные коэффициенты, m < n. Если число коэффициентов аппроксимации взять равным числу узлов, то среднеквадратичная аппроксимация совпадет с интерполяцией Лагранжа, при этом, если не учитывать вычислительную погрешность, Q = 0.

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

. (3)

Иными словами, если , число коэффициентов аппроксимации недостаточно для правильного воспроизведения графика экспериментальной зависимости. Если , многие коэффициенты в (2) не будут иметь физического смысла.

Для решения задачи линейной аппроксимации в общем случае следует найти условия минимума суммы квадратов отклонений для (2). Задачу на поиск минимума можно свести к задаче поиска корня системы уравнений , k = 0…m. (4).

Подстановка (2) в (1), а затем расчет (4) приведет в итоге к следующей системе линейных алгебраических уравнений:

Далее следует решить полученную СЛАУ относительно коэффициентов c₀…c_m. Для решения СЛАУ обычно составляется расширенная матрица коэффициентов, которую называют матрицей Грама, элементами которой являются скалярные произведения базисных функций и столбец свободных коэффициентов:

,

где , , j = 0…m, k = 0…m.

После того как с помощью, например, метода Гаусса найдены коэффициенты c₀…c_m, можно построить аппроксимирующую кривую или вычислить координаты заданной точки. Таким образом, задача аппроксимации решена.

Аппроксимация каноническим полиномом.

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

ц₀(x) = x⁰ = 1; ц₁(x) = x¹ = x; ц_m(x) = x^m, m < n.

Расширенная матрица Грама для степенного базиса будет выглядеть следующим образом:

.

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

Выбор базисных функций в виде степеней x не является оптимальным с точки зрения достижения наименьшей погрешности. Это является следствием неортогональности выбранных базисных функций. Свойство ортогональности заключается в том, что для каждого типа полинома существует отрезок [x₀, x_n], на котором обращаются в нуль скалярные произведения полиномов разного порядка:

, j ? k, с - некоторая весовая функция.

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

Блок-схема алгоритма формирования матрицы Грама и аппроксимации полиномом.

1

Аппроксимация ортогональными классическими полиномами.

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

1) Полиномы Чебышева.

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

Строятся следующим образом (рекуррентная формула):

T₀(x) = 1;

T₁(x) = x;

T_k+1(x) = 2xT_k(x) - T_k-₁(x).

2) Полиномы Лежандра.

Определены и ортогональны на [-1, 1] с весом .

Строятся следующим образом (рекуррентная формула):

L₀(x) = 1;

L₁(x) = x;

.