Методы компьютерных вычислений и их приложение к физическим задачам 2 (стр. 10 из 12)
Даже если нулевое приближение лежит в области 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) α = ½. В этом случае
