DECgroup Inc


RADAU5

RADAU5 - алгоритм решения систем обыкновенных дифференциальных уравнений любого типа (в том числе и "жестких") основанный на неявном методе Рунге-Кутты пятого порядка точности с помощью квадратур Радау. Детальное описание алгоритма приведено в книге Хаирьера и Ваннера [1]. Ниже представлены некоторые резуьтаты тестирования алгоритма RADAU5, а также их сравнение с методом 'ode23s', основанного на модифицированной квадратуре Розенброка, из коммерческой системы Matlab® или открытой системы Octave.

Валидизация и верификация метода на примерах "жестких" (систем) обыкновенных дифференциальных уравнений

1. Одномерная модель распространения ламинарного пламени

При возгорании простой спички, поверхность пламени в виде шара начинает быстро расти пока не достигнет своего критического значения, после которого остается постоянной. С помощью интегральных соотношений можно показать, что количество кислорода, поглощаемого при горении внутри шара ограничено поверхностью самого шара. Простейшая математическая модель, описывающая данный процесс, следующая:
y(t) = y2 - y3
y(0) = δ; 0 < t < δ/2;
Скалярная переменная y(t) описывает радиус шара. Члены y2 и y3 представляют собой поверхность воспламенения и объем, соответственно. Критическим параметром в данной задачи является радиус воспламенения, δ, в момент времени t = 0, который считается "малым". Решение задачи находится на отрезке времени, обратно пропорциональному δ.

Данная система изначально не является "жесткой". Жесткость начинает проявляться по мере приближения системы к стационарному режиму. На рис. 1 показаны результаты интегрирования данной задачи для значений начального радиуса δ = 0.01 и δ = 0.0001.

Рис. 1. Результаты численного решения задачи о распространение пламени для δ = 0.01 и δ = 0.0001.

2. Уравнение Ван Дер Роля (van der Pol Equation)

Уравнение Ван Дер Роля имеет следующий вид: dy2/d2t = µ (1 - y2) dy/dt + y1, где µ - положительный скалярный параметр. Введя переменную dy/dt = y2, уравнение Ван Дер Роля сводится к следующий системе обыкновенных дифференциальных уравнений:
dy1/dt = y2
dy2/dt = µ (1 - y1) y2 - y1
y1(0) = 2; y2(0) = 0;

В зависимости от параметра µ, решение уравнение Ван Дер Роля может содержать сильные колебания на существенных интервалах интегрирования (т.е. система становится "жесткой"). На Рис. 2 приведены результаты интегрирования для µ = 1 и µ = 1000 с помощью алгоритмов ode23s и RADAU5.

Рис. 2. Результаты численного решения уравнения Ван Дер Роля для µ = 1 и µ = 1000.

Ссылки::

1. Hairer, E., and Wanner, G., Solving Ordinary Differential Equations II. Stiff and Differential-Algebraic Problems. Springer Series in Computational Mathematics 14, Springer-Verlag 1991, Second Edition (1996)