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

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

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

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

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

z (x0, x1, y) = [z(x0, y)-z(x1, y)]/(x0-x1)

z (x, y0, y1) = [z(x, y0)-z(x,y1)]/(y0-y1)

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

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

Постановка задачи аппроксимации по МНК. Условия наилучшего приближения.

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

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

. (1)

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

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

ц(x) = c0ц0(x) + c1ц1(x) + … + cmцm(x), (2)

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

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

. (3)

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

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

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

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

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

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

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

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

ц0(x) = x0 = 1; ц1(x) = x1 = x; цm(x) = xm, m < n.

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

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

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

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

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

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

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

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

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

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

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

T0(x) = 1;

T1(x) = x;

Tk+1(x) = 2xTk(x) - Tk-1(x).

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

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

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

L0(x) = 1;

L1(x) = x;

.
Сглаживание и линейная регрессия.

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

1) Линейная регрессия.

В случае линейного варианта МНК (линейная регрессия) ц(x) = a + bx можно сразу получить значения коэффициентов a и b по следующим формулам:

где , .

2) Линейное сглаживание по трём точкам.

3) Линейное сглаживание по пяти точкам.

10. Решение нелинейных уравнений с одним неизвестным.

Общие сведения о численном решении уравнений с одним неизвестным.

Пусть задана непрерывная функция f(x). Требуется найти корни уравнения f(x) = 0 численными методами - это и является постановкой задачи. Численное решение уравнения распадается на несколько подзадач:

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

· единственный корень;

· бесконечное множество решений;

· корней нет;

· имеется несколько решений, как действительных, так и мнимых (например, для полинома степени n). Корни четной кратности выявить сложно.

2) Локализация корней (разбиение на интервалы) и выбор начального приближения к каждому корню. В простейшем случае можно протабулировать функцию с заданным шагом.

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

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

Метод дихотомии (бисекций).

Иначе называется методом половинного деления. Пусть задан начальный интервал [x0, x1], на котором f(x0)f(x1) ? 0 (т.е. внутри имеется не менее чем один корень). Найдем x2 = ? (x0 + x1) и вычислим f(x2). Если f(x0)f(x2) ? 0, используем для дальнейшего деления отрезок [x0, x2], если > 0 - используем для дальнейшего деления отрезок [x1, x2], и продолжаем деление пополам.

Итерации продолжаются, пока длина отрезка не станет меньше 2о - заданной точности. Тогда середина последнего отрезка даст значение корня с требуемой точностью. В качестве иного критерия можно взять | f(x)| ? оy.

Скорость сходимости метода невелика, однако он прост и надежен. Метод неприменим к корням четной кратности. Если на отрезке несколько корней, то заранее неизвестно, к какому из них сойдется процесс.

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

Метод хорд.

Идея метода проиллюстрирована рисунком. Задается интервал [x0, x1], на котором f(x0)f(x1) ? 0, между точками x0 и x1 строится хорда, стягивающая f(x). Очередное приближение берется в точке x2, где хорда пересекает ось абсцисс. В качестве нового интервала для продолжения итерационного процесса выбирается тот, на концах которого функция имеет разные знаки. Условия выхода из итерационного цикла: или

| f(x)| ? оy.

Для вывода итерационной формулы процесса найдем точку пересечения хорды (описываемой уравнением прямой) с осью абсцисс: ax2 + b = 0, где ; b = f(x0) - ax0.

Отсюда легко выразить .

Метод хорд в большинстве случаев работает быстрее, чем метод дихотомии. Недостатки метода те же, что и в предыдущем случае.

Метод Ньютона (касательных).

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

Уравнение касательной, проведенной из точки (x0, f0): y(x) = f /(x0)(x-x0) + f(x0) дает для y(x1) = 0 следующее выражение:

, (1)

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

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

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

Метод секущих.

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

Для начала итерационного процесса необходимо задать x0 и x1, которые не обязательно ограничивают интервал, на котором функция должна менять знак; это могут быть любые две точки на кривой. Выход из итерационного процесса по условию .

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

От разболтки страхуются приемом Гарвика: выбирают некоторое оx и ведут итерации до выполнения условия . Затем продолжают расчет, пока убывает. Первое же возрастание может свидетельствовать о начале разболтки, а значит, расчет следует прекратить, а последнюю итерацию не использовать.

Метод простых итераций.

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

f(x) = 0 (2)

к эквивалентному уравнению x = ц (x). Этот переход можно осуществить разными способами, в зависимости от вида f(x). Например, можно положить

ц (x) = x + bf(x), (3)

где b = const, при этом корни исходного уравнения (2) не изменятся.

Если известно начальное приближение к корню x0, то новое приближение x1 = ц (x0), т.е. общая схема итерационного процесса:

xk+1 = ц (xk). (4)

Наиболее простой критерий окончания процесса .

Критерий сходимости метода простых итераций: если вблизи корня |ц/(x)| < 1, то итерации сходятся. Если указанное условие справедливо для любого x, то итерации сходятся при любом начальном приближении. Исследуем выбор константы b в функции (3) с точки зрения обеспечения максимальной скорости сходимости. В соответствии с критерием сходимости наибольшая скорость сходимости обеспечивается при |ц/(x)| = 0. При этом, исходя из (3),

b = -1/f /(x), и итерационная формула (4) переходит в

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

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

Постановка задачи

Требуется решить систему нелинейных уравнений (1). В координатном виде эту задачу можно записать так: , где 1 ? k ? n.

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

Метод Ньютона.

Обозначим некоторое приближение к корню системы уравнений . Пусть малое . Вблизи каждое уравнение системы можно линеаризовать следующим образом:

, 1 ? k ? n. (2)

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

, 1 ? k ? n. (3)

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

, (4)

где в данном случае

- матрица Якоби, которая считается для каждого (s) приближения.

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

Блок-схема метода Ньютона для решения систем нелинейных уравнений.

Так как метод Ньютона отличается высокой скоростью сходимости при выполнении условий сходимости, на практике критерием работоспособности метода является число итераций: если оно оказывается большим (для большинства задач >100), то начальное приближение выбрано плохо.

Частным случаем решения (4) методом Ньютона системы из двух нелинейных уравнений

являются следующие легко программируемые формулы итерационного процесса:

, где ,

Метод простых итераций.

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

, (5)

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

Например, для исходной системы уравнений эквивалентная итерационная система (5) может быть представлена в следующем виде:

где множители = -0.15и = -0.1 подбираются из анализа условий сходимости.

Метод спуска.

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

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

Численные методы поиска минимума функции одной переменной.

Определения.

Функция f(x) имеет локальный минимум при некотором , если существует некоторая конечная о-окрестность этого элемента, в которой , . Требуется, чтобы на множестве X функция f(x) была по крайней мере кусочно-непрерывной.

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

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

Методы поиска минимума по нахождению корней уравнений.

Если функция f(x) аналитически дифференцируема, то решаем f /(x) = 0 методами, изложенными в предыдущих главах. При этом условие f //(x) > 0 в найденной точке указывает нам на минимум. Для использования этих методов необходимо знать либо аналитический вид первой и второй производных, либо рассчитать их численно, если это не приведет к потере точности.

Метод дробления.

Наиболее простой метод поиска минимума. Пусть дана начальная точка x0, а также величина и знак шага h, определяющие движение из этой точки в сторону предполагаемого минимума f(x). Метод заключается в последовательном дроблении исходного шага h с изменением его знака при выполнении условия f(xk+1) > f(xk), где k - порядковый номер вычисляемой точки. Например, как только очередное значение функции стало больше предыдущего, выполняется h = - h/3 и процесс продолжается до тех пор, пока

|xk+1 - xk| ? о . (1)

Данный метод является одним из самых медленных для поиска минимума. Основное достоинство данного алгоритма - возможность использования в программах управления экспериментальными исследованиями, когда значения функции f(x) последовательно измеряются с шагом h ? hmin.

Метод золотого сечения.

Пусть f(x) задана и кусочно-непрерывна на [xL, xR], и имеет на этом отрезке только один локальный минимум. Золотое сечение, о котором упоминал ещё Евклид, состоит в разбиении интервала [xL, xR] точкой x1 на две части таким образом, что отношение длины всего отрезка к его большей части равно отношению большей части к меньшей:

. (2)

Таким образом, возьмем на отрезке две точки x1 и x2, симметрично относительно границ делящие исходный отрезок в отношении золотого сечения:

где коэффициент .

Если f(x1) < f(x2), мы должны сузить отрезок справа, т.е. новое значение xR = x2, в противном случае xL = x1. Оставшаяся внутри нового отрезка точка является первым приближением к минимуму и делит этот отрезок в отношении золотого сечения. Таким образом, на каждой итерации приближения к минимуму (см. рисунок) нам нужно ставить только одну точку (x1 или x2), в которой считать значение функции и сравнивать его с предыдущим. Условием выхода из итерационного процесса будет, подобно предыдущему случаю, условие |x2 - x1| ? о.

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

Метод парабол.

Пусть f(x) имеет первую и вторую производную. Разложим f(x) в ряд Тейлора в некоторой точке xk, ограничиваясь при этом тремя членами разложения:

. (3)

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

Обычно в практических реализациях данного метода не используют аналитический вид первой и второй производных f(x). Их заменяют конечно-разностными аппроксимациями. Наиболее часто берут симметричные разности с постоянным шагом h:

;

Это эквивалентно аппроксимации функции параболой, проходящей через три близкие точки xk+h, xk, xk-h. Окончательное выражение, по которому можно строить итерационный процесс, таково:

. (4)

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

1) знаменатель формулы (4) должен быть >0. Если это не так, нужно сделать шаг в обратном направлении, причем достаточно большой. Обычно в итерационном процессе полагают . Иногда ради упрощения расчетов полагают , однако это существенно уменьшает скорость сходимости.

2) . Если это не так, то от xk следует сделать шаг с ф = ?. Если и при этом условие убывания не выполнено, уменьшают ф и вновь делают шаг.

Численные методы поиска минимума функции нескольких переменных.

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

2) Истинный овраг

1) Котловинный (гладкая функция)

3) Разрешимый овраг	4) Неупорядоченный

Метод координатного спуска.

Пусть требуется найти минимум f(x, y). Выберем нулевое приближение (x0, y0). Рассмотрим функцию одной переменной f(x, y0) и найдем ее минимум, используя любой из рассмотренных выше способов. Пусть этот минимум оказался в точке (x1, y0). Теперь точно так же будем искать минимум функции одной переменной f(x1, y). Этот минимум окажется в точке (x1, y1). Одна итерация спусков завершена. Будем повторять циклы, постепенно приближаясь ко дну котловины, пока не выполнится условие .

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

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

Метод градиентного (наискорейшего) спуска.

В этом методе функция минимизируется по направлению, в котором она быстрее всего убывает, т.е. в направлении, обратном .

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

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

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

Метод оврагов.

Ставится задача найти минимум для овражной функции.

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

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

Проблемы поиска минимума в задачах с большим числом измерений.

Пусть в n-мерном векторном пространстве задана скалярная функция . Наложим дополнительные условия , 1 ? i ? m; , 1 ? j ? p. Условия типа равенств выделяют в пространстве некоторую (n - m)-мерную поверхность, а условия типа неравенств выделяют n-мерную область, ограниченную гиперповерхностями . Число таких условий может быть произвольным. Следовательно, задача inf есть поиск минимума функции n переменных в некоторой (n - m)-мерной области E. Функция может достигать минимального значения как внутри области, так и на ее границе. Однако перейти к (n - m)-мерной системе координат практически никогда не удается, поэтому спуск приходится вести во всем n-мерном пространстве.

Даже если нулевое приближение лежит в области E, естественная траектория спуска сразу выходит из этой области. Для принудительного возврата процесса в область E, например, используется метод штрафных функций: к прибавляются члены, равные нулю в E, и возрастающие при нарушении дополнительных условий (ограничений). Метод прост и универсален, однако считается недостаточно надежным. Более качественный результат дает использование симплекс-методов линейного программирования, однако данный вопрос в рамках настоящего курса не рассматривается.

13. Решение обыкновенных дифференциальных уравнений.

Постановка задачи. Типы задач для ОДУ.

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

u(p)(x) = f(x, u, u/, u//, … u(p+1)) заменой u(k)(x) = yk(x) можно свести к эквивалентной системе из p уравнений первого порядка, представленных в каноническом виде :

y/k(x) = yk+1(x) для 0 ? k ? p-2 (1)

y/p-1(x) = f(x, y0, y1, … yp-1), при этом y0(x) ? u(x). (2)

Покажем такое преобразование на примере уравнения Бесселя:

Предполагая тождественную замену y1(x) ? y(x) представим систему ОДУ в следующем виде:

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

Известно, что система p-го порядка имеет множество решений, которые в общем случае зависят от p параметров {C1, C2, … Cp}. Для определения значений этих параметров, т.е. для выделения единственного решения необходимо наложить p дополнительных условий на uk(x). По способу задания условий различают три вида задач, для которых доказано существование и единственность решения. Это

1) Задача Коши. Задается координата uk(x0) = uk0 начальной точки интегральной кривой в (p+1) - мерном пространстве (k = 1…p). Решение при этом требуется найти на некотором отрезке x0 ? x ? xmax .

2) Краевая задача. Это задача отыскания частного решения системы ОДУ на отрезке a ? x ? b, в которой дополнительные условия налагаются на значения функции uk(x) более чем в одной точке этого отрезка.

3) Задача на собственные значения. Кроме искомых функций и их производных в уравнение входят дополнительно m неизвестных параметров л1, л2, … лm, которые называются собственными значениями. Для единственности решения на некотором интервале необходимо задать p+m граничных условий.

В большинстве случаев необходимость численного решения систем ОДУ возникает в случае, когда аналитическое решение найти либо невозможно, либо нерационально, а приближенное решение (в виде набора интерполирующих функций) не дает требуемой точности. Численное решение системы ОДУ, в отличие от двух вышеприведенных случаев, никогда не покажет общего решения системы, так как даст только таблицу значений неизвестных функций, удовлетворяющих начальным условиям. По этим значениям можно построить графики данных функций или рассчитать для некоторого x > x0 соотвествующие uk(x), что обычно и требуется в физических или инженерных задачах. При этом требуется, чтобы соблюдались условия корректно поставленной задачи.

Метод Эйлера (метод ломаных).

Рассмотрим задачу Коши для уравнения первого порядка:

u/(x) = f(x, u(x)), x0 ? x ? xmax, u(x0) = u0. (3)

В окрестности точки x0 функцию u(x) разложим в ряд Тейлора:

(4)

Идея этого и последующих методов основывается на том, что если f(x,u) имеет q непрерывных производных, то в разложении можно удержать члены вплоть до O(hq+1), при этом стоящие в правой части производные можно найти, дифференцируя (3) требуемое число раз. В случае метода Эйлера ограничимся только двумя членами разложения.

Пусть h - малое приращение аргумента. Тогда (4) превратится в

Так как в соответствии с (3) u/(x0) = f(x0, u0), то .

Теперь приближенное решение в точке x1 = x0 + h можно вновь рассматривать как начальное условие, т.е. организуется расчет по следующей рекуррентной формуле:

(5),

где y0 = u0, а все yk - приближенные значения искомой функции (см. рисунок). В методе Эйлера происходит движение не по интегральной кривой, а по касательной к ней. На каждом шаге касательная находится уже для новой интегральной кривой (что и дало название методу - метод ломаных), таким образом ошибка будет возрастать с отдалением x от x0.

При приближенное решение сходится к точному равномерно с первым порядком точности. То есть, метод дает весьма низкую точность вычислений: погрешность на элементарном шаге h составляет ? h2 y//( ?(xk + xk+1)), а для всей интегральной кривой порядка h1. При h = const для оценки апостериорной погрешности может быть применена первая формула Рунге, хотя для работы метода обеспечивать равномерность шага в принципе не требуется.

Метод Эйлера легко обобщается для систем ОДУ. При этом общая схема процесса (5) может быть записана так:

(6),

где i = 1…m - число уравнений, k - номер предыдущей вычисленной точки.

Усовершенствованный метод Эйлера-Коши с уточнением.

Данный метод базируется на предыдущем, однако здесь апостериорная погрешность контролируется на каждом шаге вычисления. Как и в предыдущем случае, рассматриваем задачу Коши на сетке с постоянным шагом h. Грубое значение вычисляется по формуле (5) и затем итерационным циклом уточняется по формуле

, (7)

где m - номер итерации. Итерационный цикл повторяется до тех пор, пока . Данная формула также легко обобщается на решение систем ОДУ. Априорная погрешность метода при m = 1 на каждом шаге порядка h3.

Блок-схема метода Эйлера-Коши с уточнением.

Метод Рунге-Кутты II порядка.

Увеличение точности решения ОДУ из предыдущей задачи при заданном шаге h может быть достигнуто учетом большего количества членов разложения функции в ряд Тейлора. Для метода Рунге-Кутты второго порядка следует взять три первых коэффициента, т.е. обеспечить:

. (8)

Переходя к приближенному решению y ? u и заменяя производные в (8) конечными разностями, получаем в итоге следующее выражение:

, (9)

где 0 ? б ? 1- свободный параметр. Можно показать, что если f(x,u) непрерывна и ограничена вместе со своими вторыми производными, то решение, полученное по данной схеме, равномерно сходится к точному решению с погрешностью порядка h2.

Для параметра б наиболее часто используют следующие значения:

1) б = 1. В этом случае

Графически это уточнение можно интерпретировать так: сначала по схеме ломаных делается шаг h/2, и находится значение . В найденной точке определяется наклон касательной к интегральной кривой, который и будет определять приращение функции для целого шага, т.е. отрезок [AB] (см. рисунок) будет параллелен касательной, проведенной в точке (xk + h/2, y(xk + h/2) ) к интегральной кривой.

2) б = Ѕ. В этом случае

Можно представить, что в этом случае по методу Эйлера сначала вычисляется значение функции и наклон касательной к интегральной кривой в точке xk+1. Затем находится среднее положение касательной из сравнения соответствующих наклонов в точках xk и xk+1, которое и будет использоваться для расчета точки yk+1.

Метод Рунге-Кутты IV порядка.

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

где i = 1…p, p - число уравнений в системе.

;

k - номер точки, для которой осуществляется расчет;

;

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

Все представленные выше схемы допускают расчет с переменным шагом. Например, шаг можно уменьшить там, где функция быстро изменяется, и увеличить в обратном случае. Так, метод Рунге-Кутты-Мерсона позволяет оценивать погрешность на каждом шаге и, в зависимости от полученной оценки принимать решение об изменении шага. Автоматический выбор шага позволяет значительно сократить время вычислений.

Метод Рунге - Кутта - Мерсона.

Этот метод отличается от метода Рунге - Кутта четвертого порядка возможностью оценивать погрешность на каждом шаге и в зависимости от этого принимать решение об изменении шага. Один из вариантов формул:

;

Rn+1 = 0.2k4 - 0.3k3 - 0.1k5 - погрешность на каждом шаге.

Пусть задана максимальна погрешность . Если , h = h/2 , и n+1 цикл расчета повторяется (с точки xn, yn) c новым шагом. Если же

, h = 2h

Автоматический выбор шага позволяет значительно сократить время решения ОДУ.

Схема РКМ обобщается на системы ОДУ аналогично классической схеме Рунге - Кутта.

Метод Адамса.

Метод основан на аппроксимации интерполяционными полиномами правых частей ОДУ.

Пусть с помощью любого из методов, рассмотренных выше, вычислено решение заданного дифференциального уравнения в точках x1, x2, x3 (а в точке x0 решение и так известно - поставлена задача Коши). Полученные значения функции обозначим как y0, y1, y2, y3, а значения правой части дифференциального уравнения как f0, f1, f2, f3, где fk = f(xk, yk). Начиная с четвертой точки, на каждом шаге интегрирования дифференциального уравнения вычисления осуществляются по схеме

P(EC){m}E

где P - прогноз решения; Е - вычисление f(x,y); С - коррекция решения; m - количество итераций коррекции. Схемы такого типа называют «прогноз-коррекция»: это подразумевает сначала приблизительное вычисление решение по формуле низкого порядка, а затем уточнение с учетом полученной информации о поведении интегральной кривой.

Прогноз осуществляется по экстраполяционной формуле Адамса:

. (10)

Коррекция осуществляется по интерполяционной формуле Адамса:

. (11)

Вычисление осуществляется по формуле:

Количество итераций m ? p, где p - порядок используемого метода. В ходе каждой итерации решается нелинейное уравнение (11) относительно неизвестной y4 (обычно методом простых итераций).

Иногда в методе Адамса используется схеме PECE на каждом шаге процесса интегрирования, т.е. осуществляется только одна коррекция. В силу сложности вычислений метод используется только в мощных программных пакетах численного анализа. Формулы метода также легко переносятся на решение систем ОДУ первого порядка.

Постановка краевой задачи. Метод стрельбы.

Краевая задача - задача отыскания частного решения системы

(12),

1 ? k ? p, на отрезке a ? x ? b, на котором дополнительные условия налагаются на значения функции uk(x) более чем в одной точке этого отрезка. В качестве примера моно привести задачу нахождения статического прогиба u(x) нагруженной струны с закрепленными концами:

, a ? x ? b, u(a) = u(b) = 0.

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

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

Рассмотрим систему двух дифференциальных уравнений первого порядка

(13)

с краевыми условиями

(14)

Сначала выберем некоторое значение , так что . Это уравнение мы можем решить и определить . Таким образом у нас появились два числа и в, которые будут определять задачу Коши , для системы (13), удовлетворять левому краевому условию, но не удовлетворять второму уравнению (14). Тем не менее, решим систему ОДУ с задачей Коши каким-либо из известных нам численных методов. Получим решение вида

(15), зависящее от как от параметра.

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

Если это осуществилось, пара таких значений определяет интервал , который можно обработать методом дихотомии до получения корня уравнения (16). Нахождение каждого нового значения функции (16) требует нового численного интегрирования для решения ОДУ, что делает этот метод достаточно медленным.

Решение уравнений в частных производных.

К уравнениям в частных производных приводят задачи газодинамики, теплопроводности, переноса излучения, электромагнитных полей, процессов переноса в газах, и др. Независимыми переменными в физических задачах обычно являются время t, координаты , скорости частиц . Пример - уравнение теплопроводности

, (17)

где U - температура, - теплоемкость, - коэффициент теплопроводности, q - плотность источников тепла.

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

Рассмотрим в качестве примера одномерную задачу, близкую по смыслу к (17):

. (18)

Здесь 0 ? x ? a, 0 ? t ? T.

Граничные условия:

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

Введем равномерную прямоугольную сетку по x и t с шагом h и ф соответственно и заменить производные соответствующими разностными отношениями. Тогда

. (19)

Здесь 1 ? k ? N-1 - число точек по координате x; 0 ? m ? M - число точек по координате t. Число неизвестных в (19) больше числа уравнений, недостающие уравнения выводятся из начальных и граничных условий:

; 0 ? k ? N.

; ; 0 ? m ? M.

Схема (19) содержит в каждом уравнении несколько неизвестных значений функции. Такие схемы называют неявными. Для фактического вычисления решения следует переписать схему так:

, где 1 ? k ? N-1.

; . (20)

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

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

Приложение 1.

Случайные величины в компьютерном моделировании.

Параметры случайных величин.

Величину о называют случайной с плотностью распределения с(x), если вероятность того, что величина примет значения между x1 и x2, равна

Псевдослучайные распределения. Иллюстрация возможностей Mathcad.

Моделирование броуновского движения в системе MathLab.

Литература

1. Н. Бахвалов, Н. Жидков, Г. Кобельков. Численные методы. М., 2002, 632 с.

2. Н. Калиткин. Численные методы. М., 1972,

3. А. Мудров. Численные методы для ПЭВМ на языках Бейсик, Фортран и Паскаль. Томск, 1992, 270с.

4. А. Самарский. Введение в численные методы. М., , 270с.

5. Ю. Тарасевич. Численные методы на Mathcad'е. Астрахань, 2000, 70с.

6. Г. Коткин, В. Черкасский. Компьютерное моделирование физических процессов с использованием MathLab. Новосибирск, 2001.

7. М. Лапчик, М. Рагулина, Е. Хеннер. Численные методы.М., 2004, 384с.

Страницы: 1, 2