Исследование решателей переменной-шага с использованием жесткой модели

Этот пример показывает поведение решателей с переменным шагом в модели маятника Фуко. Решатели Simulink ® ode45, ode15s, ode23, и ode23t используются как тесты. Жесткие дифференциальные уравнения используются, чтобы решить эту задачу. Точное определение жесткости для уравнений отсутствует. Некоторые численные методы нестабильны, когда используются для решения жестких уравнений, и очень маленькие размеры шага требуются, чтобы получить численно стабильное решение жесткой задачи. Жесткая задача может иметь быстро изменяющийся компонент и медленно изменяющийся компонент.

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

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

Модель маятника Фуко

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

$$
\ddot{x} = 2\Omega \sin{\lambda} \dot{y} - \frac{g}{L} x
$$

$$
\ddot{y} = - 2\Omega\sin{\lambda} \dot{x} - \frac{g}{L} y
$$

$$\Omega = \mbox{ Earth's angular velocity of rotation around its axis}$$

$$\lambda = \mbox{ geographic latitude}$$

$$g = \mbox{ acceleration of gravity}$$

$$x \mbox{ and } y = \mbox{ coordinates of the pendulum bob}$$

Модель, показанная на фигуре 1, используется, чтобы решить дифференциальные уравнения, описывающие маятник Фуко. Тип sldemo_solvers в MATLAB ® Командное Окно, чтобы открыть модель. Пример описывает маятник Фуко в течение 86400 секунд. Константы и начальные условия сохраняются в рабочем пространстве модели.

Фигура 1: модель маятника Фуко, используемая верхней частью оценки эффективности решателя

Переменные решатели шага

Этот пример исследует эффективность ode45, ode15s, ode23, и ode23t решатели с переменным шагом. Чтобы прочитать больше о конкретном решателе, например ode45, тип help ode45 в Командном Окне MATLAB.

Фигура 2: Список решателей с переменным шагом, доступных в Simulink

Оценка эффективности решателя

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

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

Общее энергосбережение

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

Этот пример вычисляет нормированную общую энергию маятника на каждом временном шаге. Относительная погрешность в энергии равняется изменению общей энергии в течение симуляции. В идеале относительная погрешность в энергии должна быть нулем, потому что энергия сохранена. Общая энергия - это сумма потенциальной и кинетической энергии. Блок «NormalizeEnergy» вычисляет нормированную энергию маятника с помощью выражений, приведенных ниже.

$$ E = \frac{v_x^2 + v_y^2}{2} + g( L - \sqrt{L^2 - x^2 - y^2} ) $$

$$E_{normalized}(t) = \frac{E(t)}{E(0)}$$

$$E(0) = \mbox{ initial total energy}$$

Фигура 3: Нормированная энергия в зависимости от времени

Рисунок 3 показывает график нормированной энергии от времени, рассчитанный с помощью ode23 и ode23t. Понятно, что именно в этой проблеме ode23t намного точнее, чем ode23. В симуляции, которая использовала ode23нормированная энергия маятника уменьшилась более чем на 60%. Маятник с меньшей энергией имеет меньшую амплитуду колебаний. Вы можете увидеть этот эффект на фигуре 4, где амплитуда маятника вычисляется ode23 уменьшается при вращении плоскости колебаний.

Фигура 4: Положение маятника, рассчитанное с помощью ode23 и ode23t

Фигура иллюстрирует различие между жестким и нежестким решателем. Каждый график показывает положение маятника bob в течение дня (каждая 15-я точка данных изображена как синяя точка). Черная точка помечает начальное положение маятника, а черная линия - начальную плоскость колебаний маятника. Красная линия указывает плоскость колебаний через день. Желтая линия показывает плоскость колебаний в какую-то промежуточную точку времени. Обратите внимание, что плоскость колебаний маятника не завершает полное вращение в течение дня. Скорость вращения плоскости колебаний зависит от географической широты (см. детали в примере маятника Фуко). Заметьте, что амплитуда маятника на левом графике уменьшается, в то время как амплитуда на правом графике остается постоянной. Для той же относительной погрешности RelTol=1e-5, жесткий решатель дает более точный результат, но требует большего времени выполнения.

Рисунок 5 показывает более детальное исследование эффективности решателей Simulink. Четыре решателя были выбраны, чтобы проиллюстрировать, как относительные погрешности и время выполнения симуляции варьируются как функция относительной погрешности. Можно использовать sldemo_solvers_mcode.m скрипт для более обширных тестов решателя. Этот скрипт генерирует код С из модели, чтобы ускорить симуляции. Обратите внимание, что выполнение скрипта может занять несколько минут.

Building RSim executable for model..
Time taken: 16.718875 seconds.

Running generated RSim executable..
Time taken: 377.240099 seconds.

Фигура 5: Относительная погрешность и время выполнения для различных настроек решателя. AbsTol = '1e-7' для всех симуляций.

Обратите внимание, что в этом примере относительная погрешность не очень сильно уменьшается для относительных погрешностей ниже 1e-6. Это ограничение числового решателя, которое зависит от конкретной модели. Уменьшение относительной погрешности решателя не обязательно улучшает точность. Вам нужно оценить минимальную точность, которая требуется для вашей задачи, и выбрать решатель соответственно, чтобы минимизировать затраты на симуляцию. Для примера вы можете захотеть узнать положение маятника Фуко в пределах нескольких сантиметров. Нет необходимости вычислять положение маятника в пределах нескольких микрон, потому что вы не можете точно измерить положение.

Заключения

Числовые результаты симуляции могут вести себя по-разному в зависимости от настроек решателя. Это имеет решающее значение в случае жестких проблем. При работе с жесткими моделями выбирайте решатель, который даст точный результат при минимальной цене. Относительная погрешность (решатели переменного шага) или размер шага (решатели с фиксированным шагом) должны быть достаточно маленькими, чтобы результат был численно стабильным. Хотя некоторые решатели более эффективны, чем другие, жесткие решатели лучше подходят для решения жестких задач. Решатели переменного шага более устойчивы, чем решатели с фиксированным шагом.