Краткие теоретические сведения

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

Например, для дифференциального уравнения первого порядка таблица решения будет представлять собой два столбца – x и y.

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

или , i = 1, …, N

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

· Задача Коши (начальная задача): Необходимо найти такое частное решение дифференциального уравнения, которое удовлетворяет определенным начальными условиям, заданным в одной точке:

то есть, задано определенное значение независимой переменной (х0), и значение функции и всех ее производных вплоть до порядка (n-1) в этой точке. Эта точка (х0) называется начальной. Например, если решается ДУ 1-го порядка, то начальные условия выражаются в виде пары чисел (x0, y0)

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

Численные методы решения задачи Коши ОДУ первого порядка

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

(1)

Необходимо найти значения функции y в заданных точках сетки , если известны начальные значения , где есть значение функции y(x) в начальной точке x0.

Преобразуем уравнение умножением на dx

И проинтегрируем левую и правую части между i-ым и i+1-ым узлами сетки.

(2)

Мы получили выражение для построения решения в i+1 узле интегрирования через значения x и y в i-ом узле сетки. Сложность, однако, заключается в том, что интеграл в правой части есть интеграл от неявно заданной функции, нахождение которого в аналитическом виде в общем случае невозможно. Численные методы решения ОДУ различным способом аппроксимируют (приближают) значение этого интеграла для построения формул численного интегрирования ОДУ.

Из множества разработанных для решения ОДУ первого порядка методов рассмотрим метод Эйлера, Он достаточно прост и даёт начальное представление о подходах к решению данной задачи в рамках численного решения задачи Коши.

Метод Эйлера

Исторически первым и наиболее простым способом численного решения задачи Коши для ОДУ первого порядка является метод Эйлера. В его основе лежит аппроксимация производной отношением конечных приращений зависимой (y) и независимой (x) переменных между узлами равномерной сетки:

где yi+1 это искомое значение функции в точке xi+1.

Если теперь преобразовать это уравнение, и учесть равномерность сетки интегрирования, то получится итерационная формула, по которой можно вычислить yi+1 , если известно yi в точке хi:

(3)

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

Графическая интерпретация метода Эйлера также не представляет затруднений (см. рисунок ниже). Действительно, исходя из вида решаемого уравнения (1) следует, что значение есть значение производной функции y(x) в точке x=xi - , и, таким образом, равно тангенсу угла наклона каcательной, проведенной к графику функции y(x) в точке x=xi.

Из прямоугольного треугольника на рисунке можно найти

,

откуда и получается формула Эйлера. Таким образом, суть метода Эйлера заключается в замене функции y(x) на отрезке интегрирования прямой линией, касательной к графику в точке x=xi. Если искомая функция сильно отличается от линейной на отрезке интегрирования, то погрешность вычисления будет значительной. Ошибка метода Эйлера прямо пропорциональна шагу интегрирования:

Ошибка ~ h

Процесс вычислений строится следующим образом. При заданных начальных условиях x0 и y0 можно вычислить

Таким образом, строится таблица значений функции y(x) с определенным шагом (h) по x на отрезке [x0, xN]. Ошибка в определении значения y(xi) при этом будет тем меньше, чем меньше выбрана длина шага h (что определяется точностью формулы интегрирования).

При больших h метод Эйлера весьма неточен. Он дает все более точное приближение при уменьшении шага интегрирования. Если отрезок [xi, xi+1] слишком велик, то каждый участок [xi, xi+1] разбивается на N отрезков интегрирования и к каждому их них применяется формула Эйлера с шагом , то есть шаг интегрирования h берется меньше шага сетки, на которой определяется решение.

Пример1:

Используя метод Эйлера, построить приближенное решение для следующей задачи Коши:

На сетке с шагом 0,1 в интервале [0, 1] (4)

Решение:

Данное уравнение уже записано в стандартном виде, резрешенном относительно производной искомой функции.

Поэтому, для решаемого уравнения имеем

Примем шаг интегрирования равным шагу сетки h = 0,1.

И т. д.

Полные результаты (с точностью до пятого знака после запятой) приведены в таблице 1 в третьей колонке - h=0,1 (N=1). Во второй колонке таблицы для сравнения приведены значения, вычисленные по аналитическому решению данного уравнения .

Во второй части таблицы приведена относительная погрешность полученных решений. Видно, что при h=0,1 погрешность весьма велика, достигая 100% для первого узла x=0,1.

Таблица 1 Решение уравнения (4) методом Эйлера (для колонок указан шаг интегрирования и число отрезков интегрирования N между узлами сетки)

Уменьшим шаг интегрирования вдвое, h = 0.05, в этом случае для каждого узла сетки вычисление будет проводиться за два шага (N=2). Так, для первого узла x=0,1 получим:

И так далее, до конца отрезка.

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

При шаге интегрирования h=0,025 для каждого узла сетки необходимо выполнить 4 вычисления по формуле Эйлера в промежуточных точках (N=4).

(Для других узлов значения приведены в таблице 1, колонка N=4)

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

Усовершенствованный метод Эйлера ( метод Эйлера с пересчетом).

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

(5)

Данная формула оказывается неявной относительно yi+1 (это значение есть и в левой и в правой части выражения), то есть является уравением относительно yi+1, решать которое можно, например, численно, применяя какой-либо итерационный метод (в таком виде его можно рассматривать как итерационную формула метода простой итерации). Однако, можно поступить иначи и приблизительно вычислить значение функции в узле i+1 с помощью обычной формулы Эйлера:

,

которое затем использовать при вычислении по (5).

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

(6)

Благодаря более точной формуле интегрирования, погрешность метода Эйлера с пересчетом пропорциональна уже квадрату шага интегрирования.

Ошибка ~ h2

Пример2

Проведем вычисления для уравнения (4) с помощью метода Эйлера с пересчетом .

При шаге интегрирования h=0,1 в первом узле сетки x1 получим:

Что намного точнее значения, полученного методом Эйлера при том же шаге интегрирования. В таблице 2 ниже приведены сравнительные результаты вычислений при h= 0,1 методов Эйлера и Эйлера с пересчетом

Таблица 2 Решение уравнения методами Эйлера и Эйлера с пересчетом

Отметим существенное увеличение точности вычислений метод Эйлера с пересчетом по сравнению с методом Эйлера. Так, для узла x=0,1 относительное отклонение значения функции, определенного метод Эйлера с пересчетом , оказывается в 30 (!) раз меньше. Такая же точность вычислений по формуле Эйлера достигается при числе отрезков интегрирования N примерно 30. Следовательно, при использовании метода Эйлера с пересчетом при одинаковой точности вычислений понадобится примерно в 15 раз меньше времени ЭВМ, чем при использовании метода Эйлера.

Проверка устойчивости решения

Решение ОДУ в некоторой точке xi называется устойчивым, если найденное в этой точке значение функции yi мало изменяется при уменьшении шага интегрирования. Для проверки устойчивости, таким образом, надо провести два расчета значения (yi) – с шагом интегрирования h и при уменьшенной (например, двое) величине шага

В качестве критерия устойчивости можно использовать малость относительного изменения полученного решения при уменьшении шага интегрирования (ε – наперед заданная малая величина)

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