Численные методы решения стохастических дифференциальных уравнений

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

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

Последовательные приближения определяются с помощью стохастических интегралов как

Здесь — номер приближения; в качестве нулевого приближения можно взять постоянную во времени (вообще говоря, случайную) величину .

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

(19.75)

При численном решении на ЦВМ значение первого интеграла в (74) может быть получено любым известным методом численного интегрирования [112]. Для вычисления второго интеграла в (74) на основании определения симметризованного стохастического интеграла имеем

(19.76)

где — шаг дискретизации по времени; — приращения винеровского процесса. Так как при отыскании значения случайного процесса известны, то при заданных значениях приращений винеровского процесса интеграл (76) легко вычисляется. Отметим, что последовательные приближения при . следует проводить при одних и тех же значениях , поскольку под точным решением (74) обычно понимается реализация случайного процесса при условии, что задана реализация процесса (нас интересует отклик динамической системы (1) на заданное входное воздействие ). Способ задания приращений винеровского процесса зависит от существа решаемой задачи. Так, например, при реализации на ЦВМ алгоритмов оптимальной фильтрации — выборки аддитивного шума [77].

При исследовании статистической динамики методом Монте-Карло значения согласно (14.9) могут быть получены по формуле

(19.77)

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

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

Разностные методы численного решения стохастического дифференциального уравнения (1) вытекают из рассмотренных определений стохастических интегралов [61-63, 144]. При этом следует иметь в виду, что при записи стохастических дифференциальных уравнений обязательно нужно указывать, в каком смысле они понимаются (см. выше).

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

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

(19.78)

Переходя к конечным разностям, из (78) и определения стохастического интеграла Ито (11) имеем

(19.79)

где приняты обозначения ,

Можно показать, что разностная схема (79) аппроксимирует стохастическое дифференциальное уравнение (1) с точностью .

Чтобы построить более точную разностную аппроксимацию уравнения (1), заменим в (3) интеграл его приближенным значением по формуле прямоугольников, а для стохастического интеграла воспользуемся определением Стратоновича. Тогда, вместо (79) получим

(19.80)

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

Раскладывая функции и в ряд Тейлора в окрестности точки , имеем

Соотношение (81) позволяет последовательно получить различные разностные аппроксимации стохастического дифференциального уравнения (1). Подставляя в (81) разностную аппроксимацию с точностью из (79) и учитывая члены порядка , получим

Разностная аппроксимация (82) имеет точность . В работе [63] экспериментально показано, что формула (82) действительно обладает более высокой точностью по сравнению с (79). Отметим, что если функция g не зависит от , то формулы (79), (82) совпадают и имеют одинаковую точность

Подставив в (81) значение из (82) и учитывая члены порядка , будем иметь

(19.83)

Соотношение (83) аппроксимирует стохастическое дифференциальное уравнение (1) с точностью . Аналогично из (81) и (83) можно получить разностную аппроксимацию (1) с точностью . Из-за громоздкости получаемой при этом формулы приведем ее вид для частного случая . Когда функция не зависит от , разностная аппроксимация (1) с точностью имеет вид

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

Для численного решения детерминированных дифференциальных уравнений широкое применение находят методы Рунге—Кутта [42], которые имеют точность до . Однако непосредственное применение этих методов для решения стохастического дифференциального уравнения (1) невозможно вследствие недифференцируемости процесса . Эту трудность можно обойти, если вместо (1) рассматривать уравнение

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

Для решения уравнения (85) при можно применять любые известные методы Рунге—Кутта [42]. Например, метод Рунге—Кутта четвертого порядка применительно к (85) дает следующий алгоритм:

(19-86)

Здесь — коррелированные выборки случайного процесса в соответствующие моменты времени.

При исследовании статистической динамики методом Монте— Карло в (86) для некоррелированных выборок процесса нужно положить , где задается выражением (77). Это соответствует ступенчатой аппроксимации на каждом отрезке , постоянной случайной величиной [144]. Сходимость (86) доказывается при помощи разложения коэффициентов в ряд с последующим переходом к пределу.

Сходимость рассмотренных численных методов в среднеквадратичном обеспечивает убывание среднеквадратичной погрешности (среднего квадрата ошибки) при уменьшении шага дискретизации, т. е. дает возможность воспроизвести отклик системы (1) на входное воздействие с заданной точностью. В ряде задач, когда требуется смоделировать случайный процесс с заданными статистическими характеристиками, а не воспроизвести отклик системы (1), можно воспользоваться другими определениями сходимости. Например, можно потребовать, чтобы первые моментов распределения численного решения отличались от истинных не более, чем на заданную величину. При этом точность рассмотренных численных методов, очевидно, возрастет (например, схема (79) будет иметь точность . В Приложении 1 рассмотрен численный метод решения стохастического дифференциального уравнения (1—40), который обеспечивает точное воспроизведение заданных статистических характеристик нормального процесса, но принципиально не позволяет воспроизвести отклик на входное воздействие (оно заменяется статистически эквивалентным).

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