Численное решение дифференциальных уравнений

Моделирование линейных непрерывных систем

 

При цифровом моделировании непрерывных систем необходимо обеспечить близость процессов в моделируемой непрерывной системе и в ее цифровой модели. Несовпадение этих процессов связано с двумя причинами:

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

Математическая модель непрерывной системы представляет собой или нелинейное дифференциальное уравнение или совокупность соединенных между собой линейных и нелинейных блоков. В зависимости от принятой математической модели используются различные подходы к формированию цифровой модели.

Численное решение дифференциальных уравнений

 

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

du / dt = f (u,x,t). (1)

 

Здесь x = x (t) - независимая функция (входной процесс), u = u (t) - решение уравнения (выходной процесс).

Численное решение находится для дискретных значений аргумента t, отличающихся на шаг интегрирования D t. В одношаговых разностных методах для нахождения следующего значения u _к = u (t _к) требуется информация только об одном предыдущем шаге. Из одношаговых методов наибольшую известность получили методы Рунге-Кутта. В основу метода Рунге-Кутта первого порядка, называемого также явным или прямым методом Эйлера, положено разложение функции u (t) в ряд Тейлора в окрестности точки A (t _k-1,, u _k-1):

u (t) = S₀ + S₁ (t - t_{k - 1}) + S₂ (t - t_{k - 1}) ² + …, (5.2)

где S₀ = u (t_{k - 1}) = u_{k - 1,} S_i = ( 1/ i!) du (t) / dt при t = t_{k - 1}.

 

В методах Эйлера (и Рунге-Кутта тоже) ограничиваются только двумя первыми членами разложения в ряд. Запишем значение u_k = u (t_k), приняв в выражении (5.2) t = t_k и ограничившись двумя первыми членами ряда:

u_k = u_{k - 1} + S₁ (t_k - t_{k - 1}) = u_{k - 1} + S₁ Δ t

 

Учитывая, что производная du (t) / dt равна правой части дифференциального уравнения (1), имеем S₁ = f (u_{k - 1}, x_{k - 1}, t_{k - 1}) и окончательно получим:

u_k = u_{k - 1} + Δ t f (u_{k - 1}, x_{k - 1}, t_{k - 1}). (3)

 

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

На рис. 1 а) проиллюстрировано решение прямым методом Эйлера.

 

а)	б)
Рис.1

 

Видим, что при использовании этого метода используется линейная экстраполяция и тангенс угла наклона экстраполирующей прямой равен производной функции u (t) в точке А. Экстраполированное значение u_k отличается от точного на величину ошибки.

Неявный (обратный) метод Эйлера основан на разложении функции u (t) в ряд Тейлора в окрестности точки В (u _k,, t _k) (см. рис.1 б):

u (t) = u_k + S₁ (t - t_k) + S₂ (t - t_k) ² + …,

 

Приняв в этом выражении t = t_{k - 1} и ограничившись двумя первыми членами ряда, получим

u_{k - 1} = u_k - Δ t f (u_k, x_k, t_{k )}.

 

Откуда

u_k = u_{k - 1} + Δ t f (u_k, x_k, t_{k )}. (5.4)

 

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

Применим методы Эйлера для расчета переходной характеристики интегрирующей цепи. Передаточная функция интегрирующей цепи:

K (p) = 1/ (1 + pT).

 

Отсюда дифференциальное уравнение в операторной форме:

 

(pT + 1) y = x

и в канонической форме:

Tdy / dt + y = x.

 

Перепишем его в виде (1):

dy / dt = (1/ T) (x - y).

 

Запишем рекуррентную формулу для прямого метода Эйлера в соответствии с (5.3)

y_k = y_{k - 1} + (Δ t / T) (x_{k - 1} - y_{k - 1}), (5.5) или y_k = (1 - Δ t / T) y_{k - 1} + Δ t / T x_{k - 1}.

 

Формула для обратного метода Эйлера запишется в соответствии с (4)

 

y_k = y_{k - 1} + (Δ t / T) (x_k - y_k).

 

Так как уравнение линейное, то значение y_k вычисляется в явной форме:

y_k = (y_{k - 1} + (Δ t / T) x_k) / (1 + Δ t / T). (6)

 

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

Для вычислений по формулам (5.5) и (5.6) используем структуру Formula Node. Внутри этой структуры запишем точное выражение для переходной характеристики:

 

z = 1 - e^{- i Δ t / T},

 

и выражения для переходной характеристики, полученные прямым методом Эйлера:

y = y ₁ + (Δ t / T) (1 - y ₁)

 

и обратным методом Эйлера:

 

v = (v ₁ + (Δ t / T)) / (1 + Δ t / T)

 

при нулевых начальных условиях: y (0) = 0, v (0) = 0.

В этих выражениях использованы различные обозначения для выходных переменных и принято x = 1 (t) = 1, так как t > 0.

На рис.2 показана эта структура. В формулах Δ t обозначена как dt.

 

Рис.2	Рис.3

 

Напомним, что для образования входных и выходных терминалов нужно щелкнуть ПКМ на границе структуры в предполагаемом месте терминала и в раскрывшемся меню выбрать Add Input или Add Output.

Для формирования массивов выходных переменных структура Formula Node помещается внутрь структуры For Loop, при этом задержанные на интервал дискретизации отсчеты выходных переменных y ₁ и v ₁ получаются с помощью регистра сдвига (рис.3).

Прямой метод Эйлера при большом интервале дискретизации может дать неустойчивое решение. Это случится, если отклонение решения от входного процесса x_{k - 1} - y_{k - 1} (см формулу (5)) даст такое значение y_k. что отклонение на следующем шаге x_k - y_k будет той же величины, что и предыдущее, но обратным по знаку. Решение будет колебательным незатухающим.

 

К графическому индикатору

Рис.4

 

В предыдущих лабораторных работах развертка графического индикатора Graph осуществлялась автоматически в соответствии с типом данных, подаваемых на вход графического индикатора. В этой работе мы сформируем данные так, чтобы по горизонтальной оси откладывалось время. Для этого надо сформировать кластер, куда кроме массива данных будет входить информация о времени. Используем ВП Bundle (Объединить), который находится в подпалитре Cluster (Кластер). На его входы element подаются (см. рис.4): на верхний - время начала развертки - 0; на средний - интервал дискретизации - Δt; на нижний - массив данных