Решение обыкновенных дифференциальных уравнений

Глава 4

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

Обычно ОДУ можно представить в следующем виде:

- правая часть является явной функцией времени. (4.1)

- правая часть зависит от переменной и времени. (4.2)

В технике - в правую часть входит входная переменная u, (4.3)

а переменная x является выходной.

Если уравнения записаны в виде системы уравнений первого порядка относительно производных, то говорят, что они записаны в форме Коши. Если утверждать, что уравнения записаны в векторной форме, то мы имеем дело с системой дифференциальных уравнений. Такая система может быть представлена одним дифференциальным уравнением n -го порядка.

Как правило, так записываются нелинейные дифференциальные уравнения. Если f(t) – имеет аналитическую форму записи, то уравнение (4.1) допускает такую форму записи:

.

Дифференциальное уравнение, таким образом, связано с интегралом. Следует помнить, что результат решения дифференциального уравнения – это функция времени x(t), полученная для заданных начальных условий. А определенный интеграл – это число, эквивалентное площади под кривой, стоящей под интегралом. Следовательно, нас интересует аналитическое значение интеграла. Кроме того, чтобы определить постоянную интегрирования, необходимо задать начальные условия, т.е. при t=t₀, обычно t₀=0, необходимо задать значение x(t₀)=x₀. Начальные условия необходимом задавать всегда при решении уравнений и аналитическими, и численными методами.

Исследование численных методов предполагает дискретизацию математической модели. Это означает, что искомая функция x(t) будет определена на сетке по времени на интервале [0, T], сетка может быть равномерной и неравномерной.

t

0 Δ tT

Алгоритмы на равномерной сетке проще. Δ t – шаг сетки, шаг интегрирования. Значения функции вне узлов не определены, но их можно найти интерполяцией. Вывод результатов в виде таблиц предполагает, как правило, равномерный шаг сетки. При построении графика можно использовать все имеющиеся узлы. В профессиональных программах чаще всего используется Δ t ~ var.

Большинство численных методов используют форму Коши. При этом в реальных задачах используются, чаще всего, два варианта формы записи:

1) Линейная или нелинейная системы в скалярной форме.

2) Линейная системы в векторно-матричной форме.

Изучение численных методов начнем с самого простого - метода Эйлера.

4.1. Метод Эйлера.

Для уравнения необходимо получить решение x(t) с учетом начального условия . Введем равномерную сетку по времени t = t+Δ t и представим x (t+Δ t) в виде ряда Тейлора: . Нелинейные члены разложения ряда не учитываем. Можно записать более кратко: , где значение производной в точке t. Таким образом, в методе Эйлера используется линейное соотношение, то есть «движение» происходит не по интегральной кривой, а по прямой, которая является касательной к кривой в точке t. Поэтому ясно, что решение получится приближенным и порядок точности его пропорционален . Ошибку ограничения можно оценить по первому отброшенному слагаемому ряда Тейлора. Его можно использовать для выбора переменного шага интегрирования с заданной погрешностью ε: . В алгоритме метода Эйлера ∆ t единственный параметр, который можно изменять. Запишем приближенное значение второй производной . При переходе к следующему моменту времени можно оценить и в зависимости от этого изменять его значение.

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

Уменьшение шага уменьшает ошибку, но проблему «ухода» от точной кривой решить невозможно. При этом, ошибка будет возрастать.

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

Проиллюстрируем иную ситуацию на следующих примерах:

Для звена первого порядка запишем аналитическое решение для свободного движения при различных начальных условиях и для возмущенного движения (u =1) при различных начальных условиях Графики показаны на рис 4.2 и 4.3. Используя метод Эйлера, для представленных интегральных кривых, можно получить решение. Процесс численного решения будет сходится к установившемуся режиму, так как производные с течением времени стремятся к нулю.

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

Рассмотрим звено 1-го порядка: . Аналитическое решение:

, (4.4)

где - корень характеристического уравнения. Начальное условие нулевое На входе звена единичный скачок .

Его структурная схема:

x₁ x₂

Пусть k = t = 1.

Рис. 4.4

Запишем дифференциальное уравнение в форме Коши: Пусть шаг интегрирования Dt = 1 с.

Счёт по шагам:

, так как x₀ = 0; .

, так как .

; так как и т.д.

Составим таблицу 4.1, в которую поместим результаты расчетов по точной формуле (4.4) и по алгоритму Эйлера.

Запишем алгоритм Эйлера, используя MATLAB:

function Eyler1

%Метод Эйлера

%Tau=1; k=1; u=1;

dt=1; x=0; t=0;

display(x)

display(t)

for m=1: 6

dx=1-x;

x=x+dx*dt

t=t+dt

end

Таблица 4.1

Для того, чтобы определить влияние шага интегрирования на характер протекания переходной функции, проведем ее расчет для того же звена методом Эйлера с шагом Dt=2 с, использовав старые значения k = t = 1; u = 1[t]. Получим

Переходная функция имеет характер гармонических колебаний с амплитудой равной 2 и периодом 4 с.

Совершенно очевидно, что от величины шага интегрирования зависит не только точность расчета переходной функции, но и характер переходного процесса. Но к действительным свойствам данного звена это не имеет ни малейшего отношения. Речь идет о чисто вычислительных проблемах выбора шага интегрирования. Имея программу расчета в виде.m файла Eyler1, легко убедиться, что при шаге интегрирования Dt< 1 с, будет обеспечен требуемый характер протекания переходной функции. При 1< Dt< 2 получим затухающие колебания, а при Dt> 2 с процесс будет неустойчивым. Для численной неустойчивости характерна “раскачка” колебаний на первых шагах счета.

Запишем однородное дифференциальное уравнение того же звена: .

Его общее решение называется свободным движением системы, поскольку отсутствует возмущение на входе звена, оно запишется: . Для того же начального условия t = 0, x(0)= 0, получим решение: , где l= -1/τ корень характеристического уравнения. Знак минус содержится в l.

Дискретизируем х₂, введя t = iΔ t, где i = 0, 1, …, n и разложим экспоненту в ряд Тейлора:

Отбросив нелинейные слагаемые ряда, получим:

.

Следовательно, метод Эйлера имеет порядок точности 0(Dt). Очевидно, что сходимость решения можно получить, если коэффициент |1+lΔ t | по модулю меньше единицы. Решая это неравенство для положительного значения 1+lΔ t, получим lΔ t < 0. Т.е. процесс вычислений сходится, если звено устойчиво l < 0. Для отрицательного значения 1+lΔ t, получим -2 < lΔ t. Следовательно, для обеспечения устойчивого счета следует выбирать значение шага интегрирования на основе неравенства -2 < lΔ t < 0. Проведя пробные расчеты, мы выяснили, что значение Δ t следует выбирать Δ t < τ. При таком выборе шага интегрирования сохраняется характер протекания переходной функции.

Таким образом, для системы уравнений величину шага интегрирования будут определять звенья с наименьшими постоянными времени. И проблема состоит в том, что, даже после затухания быстрых компонент, нельзя увеличивать шаг интегрирования без учета условий устойчивости. Приходится считать всю систему уравнений с малым Δ t и затраты времени при проведении расчетов возрастают. Для систем, у которых постоянные времени отличаются на порядки (жесткие системы), разрабатывают иные методы решения дифференциальных уравнений. Мы их не рассматриваем. При этом ясно, что более точные методы легче «приспособить» для решения жёстких систем обыкновенных дифференциальных уравнений (ОДУ). При обзоре методов интегрирования, используемых в MATLAB, об этих методах будет сказано.

Если рассмотреть систему однородных дифференциальных уравнений n порядка и записать ее детерминант det|Ep-A|, то приравняв его нулю получим характеристическое уравнение det|Ep-A|=0. Матрица А, в качестве собственных значений, имеет корни характеристического уравнения. Разложение матричной экспоненты в ряд Тейлора позволяет получить аналогичный результат и для системы уравнений / E+AΔ t/< 1. Поэтому все сказанное для одного дифференциального уравненияимеет смысл для системы в целом.

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

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

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