Методы компьютерных вычислений и их приложение к физическим задачам

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. Для практических вычислений используются следующие соотношения, обобщенные в данном случае на решение системы ОДУ:

 Скачать реферат