Численное интегрирование нелинейного дифференциального уравнения

1. Теоретическое описание метода

1.1. Постановка задачи

Рассматривается колебательная система с одной степенью свободы. Тело массы \(m\) соединено с неподвижным основанием посредством упругой пружины. Движение тела происходит в среде с вязким сопротивлением под действием внешней силы (рис. 1.1).

Рис. 1.1. Расчетная схема колебательной системы

В начальный момент времени заданы начальное положение тела и его начальная скорость. Требуется определить перемещение, скорость и ускорение тела для каждого момента времени в интервале \([t_0, t_{\text{конеч}}]\).

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

Рис. 1.2. Характеристики элементов системы: (а) упругая сила; (б) сила вязкого сопротивления; (в) внешняя сила

1.2. Вывод дифференциального уравнения

Согласно второму закону Ньютона, движение системы описывается следующим соотношением:

\[ma = F^c + F^y + F^в \qquad (1.1)\]

или в эквивалентной форме:

\[F^с + F^y + F^в + F^u = 0, \qquad (1.2)\]

где с учетом выбранных направлений (рис. 1.3):

Рис. 1.3. Направления действующих сил

\[\begin{split}\begin{cases} F^c = Q \sin \dfrac{2\pi}{T} t \\[6pt] F^у = -kx \\[6pt] F^в = -\mu \nu \left| \nu \right| \\[6pt] F^u = -ma \end{cases} \qquad (1.3)\end{split}\]

Здесь \(x\), \(\nu\) и \(a\) обозначают соответственно перемещение, скорость и ускорение тела.

Подстановка соотношений (1.3) в уравнение (1.2) приводит к выражению:

\[-kx - \mu \nu \left| \nu \right| - ma + Q \sin \frac{2\pi}{T} t = 0 \qquad (1.4)\]

Принимая во внимание кинематические соотношения:

\[\nu = \frac{dx}{dt}, \qquad a = \frac{d^{2}x}{dt^{2}}, \qquad (1.5)\]

получаем дифференциальное уравнение движения тела:

\[kx + \mu \frac{dx}{dt} \left| \frac{dx}{dt} \right| + m \frac{d^{2}x}{dt^{2}} - Q \sin \frac{2\pi}{T} t = 0 \qquad (1.6)\]

1.3. Численное интегрирование

Применение численного подхода к интегрированию уравнения (1.6) предполагает нахождение приближенного решения для дискретного набора моментов времени. Временная ось представляется совокупностью точек \(t_0, t_1, t_2, \dots, t_i, t_{i+1}, \dots, t_n\) (рис. 1.4), в каждой из которых определяется приближенное решение.

Рис. 1.4. Дискретизация временной оси

Интегрирование осуществляется последовательно; выбор величины очередного шага \(\Delta t_i\) зависит как от требуемой точности, так и от результатов интегрирования на предшествующих шагах.

Таким образом, использование численного подхода позволяет осуществить переход от непрерывных функций \(x(t)\), \(\nu(t)\), \(a(t)\) на всем интервале \([t_0, t_{\text{конеч}}]\) к совокупности дискретных значений \(x_i\), \(\nu_i\), \(a_i\) для определенных моментов времени \(t_i\). При этом дифференциальные соотношения (1.5) заменяются алгебраическими формулами выбранного метода интегрирования.

Формулы неявного одношагового метода Штермера устанавливают следующие зависимости для переменных \(x_i\), \(\nu_i\) по известным с предыдущего шага значениям \(x_{i-1}\), \(\nu_{i-1}\) [1]:

\[\begin{split}\begin{cases} x_i = x_{i-1} + \nu_{i-1} \Delta t_i + a_i \dfrac{\Delta t_i^{2}}{2} \\[10pt] \nu_i = \nu_{i-1} + a_i \Delta t_i \end{cases} \qquad (1.7)\end{split}\]

где

\[\Delta t_i = t_i - t_{i-1}, \qquad i = \overline{1,n},\]

\(t_0\) — начальный момент времени,

\(x_0\), \(\nu_0\) — начальные значения перемещения и скорости,

\(t_n\) — конечный момент времени.

Для начального момента времени \(t_0\) значения \(x_0\) и \(\nu_0\) полагаются известными. Определим значения \(x_1\), \(\nu_1\), \(a_1\) для момента времени \(t_1 = t_0 + \Delta t_1\).

Уравнение (1.4) для момента времени \(t_i\) записывается в виде:

\[k x_i + \mu \nu_i \left| \nu_i \right| + m a_i - Q \sin \frac{2\pi}{T} t_i = 0 \qquad (1.8)\]

Дополняя это соотношение формулами метода интегрирования (1.7), получаем для момента времени \(t_1\) замкнутую систему уравнений:

\[\begin{split}\begin{cases} k x_1 + \mu \nu_1 \left| \nu_1 \right| + m a_1 - Q \sin \dfrac{2\pi}{T} t_1 = 0 \\[10pt] x_1 = x_0 + \nu_0 \Delta t_1 + a_1 \dfrac{\Delta t_1^{2}}{2} \\[10pt] \nu_1 = \nu_0 + a_1 \Delta t_1 \end{cases} \qquad (1.9)\end{split}\]

Приведем полученную систему к одному уравнению, выразив неизвестные \(x_1\) и \(a_1\) через \(\nu_1\):

\[\begin{split}\begin{cases} a_1 = \dfrac{\nu_1 - \nu_0}{\Delta t_1} \\[10pt] x_1 = x_0 + \dfrac{\nu_0 + \nu_1}{2} \Delta t_1 \end{cases} \qquad (1.10)\end{split}\]

В результате подстановки получаем:

\[k \left( x_0 + \frac{\nu_0 + \nu_1}{2} \Delta t_1 \right) + \mu \nu_1 \left| \nu_1 \right| + m \frac{\nu_1 - \nu_0}{\Delta t_1} - Q \sin \frac{2\pi}{T} t_1 = 0\]

Группируя сомножители при одинаковых степенях неизвестного \(\nu_1\), соотношение может быть представлено в виде:

\[\alpha \nu_1 \left| \nu_1 \right| + \beta \nu_1 + \gamma = 0, \qquad (1.12)\]

где

\[\alpha = \mu\]

\[\beta = \frac{k \Delta t_1}{2} + \frac{m}{\Delta t_1} \qquad (1.12а)\]

\[\gamma = k x_0 + \frac{k \nu_0 \Delta t_1}{2} - \frac{m \nu_0}{\Delta t_1} - Q \sin \frac{2\pi}{T} t_1\]

Соотношение (1.12) сохраняет свою форму для любого момента времени \(t_i\) при соответствующей замене индексов (\(1 \to i\), \(0 \to i-1\)).

Таким образом, применение формул метода интегрирования позволяет преобразовать исходное дифференциальное уравнение (1.6) к нелинейному уравнению вида (1.12), которое решается на каждом временном шаге.

1.4. Решение нелинейного уравнения методом Ньютона

Для решения уравнения (1.12) применяется метод Ньютона. Рассмотрим последовательность действий при решении нелинейного уравнения.

Рассматривается уравнение вида:

\[f(z) = 0, \qquad (1.13)\]

где \(f(z)\) — нелинейная функция относительно неизвестной переменной \(z\).

Рис. 1.5. Геометрическая интерпретация метода Ньютона: (а) построение касательной; (б) проверка сходимости

Алгоритм численного решения включает следующие этапы:

1. Выбор начального приближения к решению — величины \(z^0\).

2. Организация итерационной последовательности, на каждой из которых производится уточнение значения \(z\) согласно схеме (рис. 1.5а):

   \[\begin{split}\begin{cases} z^{j} = z^{j-1} + \Delta z^{j} \\[6pt] \Delta z^{j} = -\dfrac{f(z^{j-1})}{f'(z^{j-1})} \end{cases} \qquad (1.14)\end{split}\]

   где \(f(z^{j-1})\) — значение функции \(f(z)\) при \(z = z^{j-1}\), \(f'(z^{j-1})\) — значение производной \(df/dz\) при \(z = z^{j-1}\).

3. Проверка на каждой итерации критерия прекращения итерационного процесса (рис. 1.5б):

   \[|z^{j} - z^{j-1}| = |\Delta z^{j}| \leq \delta_z \qquad (1.15)\]
   \[|f(z^{j})| \leq \delta_f\]

   где \(\delta_f\) — допустимая невязка (отклонение от нуля) правой части уравнения (1.13); \(\delta_z\) — допустимая величина изменения решения на двух соседних итерациях.

4. Контроль соблюдения ограничения на максимальное количество итераций:

   \[j \leq j_{\text{max}} \qquad (1.16)\]

Геометрическая интерпретация метода Ньютона заключается в отыскании абсциссы точки пересечения кривой \(f(z)\) с осью \(z\). На каждой \(j\)-й итерации решение аппроксимируется точкой пересечения касательной к кривой \(f(z)\) в точке \(z = z^{j-1}\) с осью \(z\).

2. Численная реализация метода интегрирования

2.1. Первый этап: численный расчет

Возвращаясь к численному решению уравнения (1.12), обозначим \(z = \nu_1\). Тогда:

\[\alpha z |z| + \beta z + \gamma = 0 \qquad (1.17)\]

или

\[f(z) = 0,\]

где

\[f(z) = \alpha z |z| + \beta z + \gamma \qquad (1.18)\]

Для решения уравнения (1.17) методом Ньютона требуется выражение для производной функции \(f(z)\):

\[f'(z) = 2\alpha |z| + \beta \qquad (1.19)\]

Зададимся следующими исходными данными:

• \(k = 20000\) Н/м,

• \(\mu = 1000\) Н·с²/м²,

• \(m = 0.1\) кг,

• \(Q = 1000\),

• \(T = 0.2\pi\),

• \(F = 1000 \sin(10t)\).

Начальные условия и шаг интегрирования:

\[x_0 = 0,\quad \nu_0 = 0,\quad \Delta t_1 = 0.001\ \text{с}\]

Тогда, согласно (1.12):

\[\alpha = 1000\]

\[\beta = \frac{20000 \cdot 0.001}{2} + \frac{0.1}{0.001} = 110\]

\[\gamma = 20000 \cdot 0 + 20000 \cdot \frac{0 \cdot 0.001}{2} - \frac{0.1 \cdot 0}{0.001} - 1000 \sin(0.01) = -10\]

Таким образом,

\[f(z) = 1000 z |z| + 110 z - 10 \qquad (1.20)\]

\[f'(z) = 2000 |z| + 110 \qquad (1.21)\]

Зададим значения допустимых погрешностей:

\[\delta_z = 0.001,\qquad \delta_f = 0.1 \qquad (1.22)\]

Максимально допустимое количество итераций примем равным \(5\).

Начальное приближение к решению выберем \(z^0 = 0\).

Первая итерация.

\[\Delta z^1 = -\frac{f(z^0)}{f'(z^0)} = -\frac{1000 \cdot 0 \cdot 0 + 110 \cdot 0 - 10}{2000 \cdot |0| + 110} = 0.09091\]

\[z^1 = 0 + 0.09091 = 0.09091\]

Проверка критериев завершения итераций:

\[|\Delta z^1| > \delta_z\]

\[f(z^1) = 1000 \cdot 0.09091 \cdot |0.09091| + 110 \cdot 0.09091 - 10 = 8.26\]

\[|f(z^1)| > \delta_f\]

Требуется продолжение итерационного процесса.

Вторая итерация.

\[\Delta z^2 = -\frac{f(z^1)}{f'(z^1)} = -\frac{8.26}{2000 \cdot 0.09091 + 110} = -0.02831\]

\[z^2 = 0.09091 - 0.02831 = 0.06260\]

Проверка:

\[|\Delta z^2| > \delta_z,\qquad |f(z^2)| = 0.80 > \delta_f\]

Третья итерация.

\[\Delta z^3 = -0.00341,\quad z^3 = 0.05918\]

\[|\Delta z^3| > \delta_z,\qquad |f(z^3)| = 0.012 < \delta_f\]

Четвертая итерация.

\[\Delta z^4 = -0.00005,\quad z^4 = 0.05913\]

\[|\Delta z^4| < \delta_z,\qquad |f(z^4)| = 0.0006 < \delta_f\]

Оба условия (1.15) выполнены, ограничение (1.16) не превышено. Решение достигнуто. Вычислено значение скорости для момента времени \(t_1\):

\[\nu_1 = 0.05913\ \text{м/с}\]

Используя формулы (1.10), определим значения ускорения и перемещения:

\[a_1 = \frac{0.05913 - 0}{0.001} = 59.13\ \text{м/с}^2\]

\[x_1 = 0 + \frac{0 + 0.05913}{2} \cdot 0.001 = 2.96 \times 10^{-5}\ \text{м}\]

2.2. Второй этап: выбор величины шага

Решение для момента времени \(t_1\) получено. Выполним еще один шаг по времени для иллюстрации процедуры выбора величины шага. Уравнения (1.9)–(1.12) справедливы для любого момента времени с учетом соответствующей замены индексов. Для второго шага имеем:

\[\begin{split}\begin{cases} k x_2 + \mu \nu_2 |\nu_2| + m a_2 - Q \sin \dfrac{2\pi}{T} t_2 = 0 \\[10pt] x_2 = x_1 + \nu_1 \Delta t_2 + a_2 \dfrac{\Delta t_2^2}{2} \\[10pt] \nu_2 = \nu_1 + a_2 \Delta t_2 \end{cases} \qquad (1.23)\end{split}\]

Приводя эту систему к одному уравнению относительно скорости, получаем:

\[\alpha \nu_2 |\nu_2| + \beta \nu_2 + \gamma = 0,\]

где

\[\alpha = \mu\]

\[\beta = \frac{k \Delta t_2}{2} + \frac{m}{\Delta t_2}\]

\[\gamma = k x_1 + \frac{k \nu_1 \Delta t_2}{2} - \frac{m \nu_1}{\Delta t_2} - Q \sin \frac{2\pi}{T} t_2 \qquad (1.24)\]

Предварительно величину шага \(\Delta t_2\) примем равной \(\Delta t_1\), т.е. \(0.001\) с. Тогда, с учетом исходных данных и полученного на первом шаге решения, вычислим коэффициенты:

\[\alpha = 1000\]

\[\beta = \frac{20000 \cdot 0.001}{2} + \frac{0.1}{0.001} = 110\]

\[\begin{split}\begin{align} \gamma = & 20000 \cdot 2.96 \times 10^{-5} + 20000 \cdot \frac{0.05913 \cdot 0.001}{2} \\ & - \frac{0.1 \cdot 0.05913}{0.001} - 1000 \sin\left(\frac{2\pi}{0.2\pi} \cdot 0.002\right) = -24.7 \end{align}\end{split}\]

Вновь получаем нелинейное уравнение:

\[1000 \nu_2 |\nu_2| + 110 \nu_2 - 24.7 = 0, \qquad (1.25)\]

которое решается методом Ньютона.

Особого внимания заслуживает выбор начального приближения в алгоритме метода Ньютона. В качестве начального приближения принимается значение скорости, которое тело имело бы в момент времени \(t_2\) при условии неизменности ускорения с момента времени \(t_1\):

\[\nu_2^0 = \nu_1 + a_1 \Delta t_2 \qquad (1.26)\]

Данная величина представляет собой так называемый явный шаг (прогноз). Полученное значение скорости используется не только как начальное приближение в методе Ньютона, но и при оценке выбранного временного шага.

Итак, начальное приближение (прогноз):

\[\nu_2^0 = 0.05913 + 59.13 \cdot 0.001 = 0.11826\]

Опуская промежуточные выкладки (аналогичные приведенным для первого шага), итерации метода Ньютона дают следующую последовательность:

• начальное приближение: \(\nu_2^0 = 0.11826\)

• первая итерация: \(\nu_2^1 = 0.11172\)

• вторая итерация: \(\nu_2^2 = 0.11159\) — решение достигнуто.

Таким образом, при величине шага \(\Delta t_2 = 0.001\) с скорость для момента времени \(t_2\) составляет:

\[\nu_2 = 0.11159\ \text{м/с}\]

2.3. Оценка локальной погрешности

Погрешность метода интегрирования на \(i\)-м шаге, именуемая локальной погрешностью, оценивается по формуле:

\[lp_i = \left| \frac{\nu_i^p - \nu_i^c}{2} \right| \qquad (1.27)\]

где \(\nu_i^p\) — явный прогноз величины скорости на \(i\)-м шаге, определяемый соотношением

\[\nu_i^p = \nu_{i-1} + a_{i-1} \Delta t_i \qquad (1.28)\]

\(\nu_i^c\) — значение скорости, полученное в результате итерационного решения с использованием неявной формулы

\[\nu_i^c = \nu_{i-1} + a_i \Delta t_i \qquad (1.29)\]

Рис. 1.6. Геометрическая интерпретация локальной погрешности

Соотношение (1.28) уже применялось при выборе начального приближения (1.26). Смысл оценки локальной погрешности поясняется на рис. 1.6.

В момент времени \(t_{i-1}\) решение находится в точке \(\nu_{i-1}\). Если для вычисления \(\nu_i\) воспользоваться явной формулой (1.28), то точка \(\nu_i = \nu_i^p\) будет лежать на касательной к кривой \(\nu(t)\) в точке \(t_{i-1}\), поскольку \(a_{i-1}\) представляет собой тангенс угла наклона этой касательной.

При вычислении \(\nu_i\) с использованием формулы (1.29) требуется значение \(a_i\), т.е. тангенса угла наклона касательной к кривой \(\nu(t)\) в точке \(t_i\). Поскольку в момент времени \(t_{i-1}\) поведение функции \(\nu(t)\) при \(t = t_i\) неизвестно, значение \(\nu_i = \nu_i^c\) определяется путем совместного решения системы уравнений (1.7), включающей соотношение (1.29). Данная процедура требует проведения нескольких ньютоновских итераций.

Как показано на рис. 1.6, явный прогноз \(\nu_i^p\) и скорректированное решение \(\nu_i^c\) располагаются по разные стороны от кривой \(\nu(t)\), проходящей через точку \(t_{i-1}\). Чем больше разность между \(\nu_i^c\) и \(\nu_i^p\), тем значительнее отклонение графика скорости от линейной зависимости на текущем шаге и тем выше погрешность интегрирования. Уменьшение величины шага \(\Delta t_i\) приводит к снижению локальной погрешности, оцениваемой по формуле (1.27).

2.4. Управление величиной шага

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

Задается значение предельно допустимой локальной погрешности на шаге интегрирования \(\delta_l\). По результатам очередного \(i\)-го шага производится сравнение значений \(\delta_l\) и фактически полученной локальной погрешности \(lp_i\).

• Если \(lp_i \leq \delta_l\), то шаг признается успешным. Осуществляется переход к следующему шагу; его величина для одношаговых методов интегрирования первого порядка точности (к которым относятся формулы (1.7)) выбирается по зависимости:

  \[\Delta t_{i+1} = c \Delta t_i \sqrt{\frac{\delta_l}{lp_i}} \qquad (1.30)\]

  где \(\Delta t_i\) — величина выполненного шага, \(\Delta t_{i+1}\) — рекомендуемое значение следующего шага, \(c\) — поправочный коэффициент (\(c < 1\)).

• Если \(lp_i > \delta_l\), то значение шага \(\Delta t_i\) является чрезмерно большим и не обеспечивает требуемой точности. В этом случае необходимо выполнить повторный расчет на \(i\)-м шаге с уменьшенным значением \(\Delta t_i\). Для выбора величины шага также используется формула (1.30), однако полученное значение применяется для повторного расчета на текущем шаге.

Вернемся к рассматриваемому примеру. Для второго шага интегрирования имеем:

\[\nu_2^p = 0.11826,\qquad \nu_2^c = 0.11159\]

Локальная погрешность:

\[lp_2 = \left| \frac{0.11826 - 0.11159}{2} \right| = 0.00333\]

Приняв допустимую погрешность \(\delta_l = 0.001\), констатируем, что \(lp_2 > \delta_l\). Выполненный шаг с \(\Delta t_2 = 0.001\) с не обеспечивает требуемой точности; необходимо повторить расчет на втором шаге с уменьшенным значением \(\Delta t_2\).

Рекомендуемое значение для повторного расчета определим по формуле (1.30) с коэффициентом \(c = 0.8\):

\[\Delta t_2 = 0.8 \cdot 0.001 \cdot \sqrt{\frac{0.001}{0.00333}} = 0.438 \times 10^{-3}\ \text{с}\]

Результаты повторного расчета с шагом \(\Delta t_2 = 0.438 \times 10^{-3}\) с:

\[\nu_2^p = 0.08505,\qquad \nu_2^c = 0.08509\]

\[lp_2 = \left| \frac{0.08505 - 0.08509}{2} \right| = 0.00002\]

Поскольку \(lp_2 < \delta_l\), второй шаг признается успешным. Дополним значение скорости \(\nu_2 = 0.08509\) м/с значениями ускорения и перемещения, используя формулы связи (1.7):

\[a_2 = \frac{\nu_2 - \nu_1}{\Delta t_2} = \frac{0.08509 - 0.05913}{0.438 \times 10^{-3}} = 59.21\ \text{м/с}^2\]

\begin{aligned} x_2 &= x_1 + \nu_1 \Delta t_2 + a_2 \frac{\Delta t_2^2}{2} \\ &= 2.96 \times 10^{-5} + 0.05913 \cdot 0.438 \times 10^{-3} + 59.21 \cdot \frac{(0.438 \times 10^{-3})^2}{2} \\ &= 6.12 \times 10^{-5}\ \text{м} \end{aligned}

2.5. Геометрическая интерпретация

В момент времени \(t_{i-1}\) решение находится в точке \(\nu_{i-1}\). Через эту точку проходит кривая \(\nu(t)\) — так называемая интегральная кривая для момента времени \(t_{i-1}\), представляющая собой точное решение уравнения (1.6) при начальном условии

\[\nu|_{t = t_{i-1}} = \nu_{i-1} \qquad (1.31)\]

Поскольку уравнение (1.6) решается приближенно, на каждом \(i\)-м шаге численного интегрирования происходит переход с одной интегральной кривой, удовлетворяющей начальному условию (1.31), на другую интегральную кривую, являющуюся точным решением уравнения (1.6) при начальном условии

\[\nu|_{t = t_i} = \nu_i \qquad (1.32)\]

На рис. 1.6 интегральная кривая для \(t = t_i\) обозначена пунктирной линией.

Результатом численного решения служит ломаная линия, проходящая через совокупность интегральных кривых, каждая из которых представляет собой точное решение уравнения (1.6) при начальных условиях, определяемых численным решением на текущем шаге (рис. 1.8).

Рис. 1.8. Переход между интегральными кривыми при численном интегрировании

2.6. Основные этапы численного решения

Обобщим ключевые аспекты, имеющие значение для численного анализа рассмотренного примера:

1. Формирование дифференциального уравнения, описывающего поведение системы:

   \[kx + \mu \frac{dx}{dt} \left| \frac{dx}{dt} \right| + m \frac{d^2 x}{dt^2} - Q \sin \frac{2\pi}{T} t = 0\]

   При формировании уравнения использован второй закон Ньютона.

2. Представление уравнения в форме, не содержащей явно дифференциальных соотношений:

   \[\begin{split}\begin{cases} kx + \mu \nu |\nu| + ma - Q \sin \dfrac{2\pi}{T} t = 0 \\[8pt] \nu = \dfrac{dx}{dt} \\[8pt] a = \dfrac{d\nu}{dt} = \dfrac{d^2 x}{dt^2} \end{cases}\end{split}\]

3. Замена дифференциальных связей алгебраическими уравнениями выбранного метода интегрирования:

   \[\begin{split}\begin{cases} k x_i + \mu \nu_i |\nu_i| + m a_i - Q \sin \dfrac{2\pi}{T} t_i = 0 \\[8pt] x_i = x_{i-1} + \nu_{i-1} \Delta t_i + a_i \dfrac{\Delta t_i^2}{2} \\[8pt] \nu_i = \nu_{i-1} + a_i \Delta t_i \end{cases}\end{split}\]

4. Приведение системы к одному уравнению относительно \(\nu_i\):

   \[\alpha \nu_i |\nu_i| + \beta \nu_i + \gamma = 0\]

5. Решение нелинейного уравнения методом Ньютона. Начальное приближение определяется формулой явного прогноза.

6. Оценка точности интегрирования путем контроля локальной погрешности. При неудовлетворительной величине локальной погрешности расчет на текущем шаге повторяется с уменьшенным шагом.

7. Вычисление ускорения и перемещения по уравнениям связи после успешного завершения шага.

8. Определение величины следующего шага на основе соотношения допустимой и фактической локальной погрешности.

2.7. Автоматическое формирование математической модели

На рассмотренном примере была обозначена схема численного решения, реализуемая в вычислительном ядре PRADIS. Однако формирование дифференциального уравнения производилось «вручную», а некоторые вопросы решались неформально (например, аналитическое определение производной \(df(z)/dz\)). В связи с этим целесообразно продолжить рассмотрение методов PRADIS с выяснения принципов автоматического формирования системы дифференциальных уравнений.

Вернемся к системе, представленной на рис. 1.1. Уравнение равновесия имеет вид:

\[F^c + F^y + F^в + F^u = 0 \qquad (2.1)\]

Для каждого \(i\)-го момента времени уравнение может быть записано как:

\[F_i^c + F_i^y + F_i^в + F_i^u = 0 \qquad (2.2)\]

где

\[\begin{split}\begin{cases} F_i^c = -Q \sin \dfrac{2\pi}{T} t \\[8pt] F_i^y = k x_i \\[8pt] F_i^в = \mu \nu_i |\nu_i| \\[8pt] F_i^u = m a_i \end{cases} \qquad (2.3)\end{split}\]

(Следует обратить внимание на отличие знаков от (1.3), обусловленное тем, что в PRADIS при рассмотрении условий равновесия суммируются усилия, действующие со стороны системы на элементы.)

Связь между \(x_i\), \(\nu_i\), \(a_i\) задается формулами метода интегрирования (1.7), которые воспроизводятся ниже:

\[x_i = x_{i-1} + \nu_{i-1} \Delta t_i + a_i \frac{\Delta t_i^2}{2} \qquad (2.4)\]

\[\nu_i = \nu_{i-1} + a_i \Delta t_i\]

Соотношение (2.2) представляет собой алгебраическое нелинейное уравнение, которое может быть решено методом Ньютона. Имеем:

\[f(z) = 0, \qquad (2.5)\]

где

\[f(z) = F_i^c + F_i^y + F_i^в + F_i^u\]

В качестве переменной \(z\) может быть выбрана любая из компонент \(x_i\), \(\nu_i\) или \(a_i\), поскольку они связаны соотношениями (2.4). Примем, как и ранее, \(z = \nu_i\).

На каждой итерации необходимо вычислить \(f(z)\) и \(df(z)/dz\).

Вычисление \(f(z)\) сводится к суммированию текущих значений сил при текущих значениях \(x_i\), \(\nu_i\), \(a_i\) (\(i\) — номер шага по времени, \(j\) — номер итерации по Ньютону). Вычисляемое значение \(f(z)\) представляет собой невязку условия равновесия, которую необходимо уменьшить до допустимых пределов.

Производная \(df(z)/dz\) определяется как:

\[\frac{df(z)}{dz} = \frac{dF_i^c}{dz} + \frac{dF_i^y}{dz} + \frac{dF_i^в}{dz} + \frac{dF_i^u}{dz} \qquad (2.6)\]

Каждую из производных следует представить как производную сложной функции. Для силы \(F_i^c\):

\[\frac{dF_i^c}{dz} = \frac{\partial F_i^c}{\partial x_i} \frac{dx_i}{dz} + \frac{\partial F_i^c}{\partial \nu_i} \frac{d\nu_i}{dz} + \frac{\partial F_i^c}{\partial a_i} \frac{da_i}{dz}\]

Поскольку \(z = \nu_i\):

\[\frac{dF_i^c}{dz} = \frac{\partial F_i^c}{\partial x_i} \frac{dx_i}{d\nu_i} + \frac{\partial F_i^c}{\partial \nu_i} \frac{d\nu_i}{d\nu_i} + \frac{\partial F_i^c}{\partial a_i} \frac{da_i}{d\nu_i} \qquad (2.7)\]

Аналогично для остальных сил.

Используем уравнения связи (2.4) для получения зависимостей \(a_i\) и \(x_i\) от \(\nu_i\):

\[a_i = \frac{\nu_i - \nu_{i-1}}{\Delta t_i}\]

\[x_i = x_{i-1} + \frac{\nu_i + \nu_{i-1}}{2} \Delta t_i \qquad (2.8)\]

Дифференцируя (2.8) по \(\nu_i\), получаем:

\[\frac{dx_i}{d\nu_i} = \frac{\Delta t_i}{2},\qquad \frac{d\nu_i}{d\nu_i} = 1,\qquad \frac{da_i}{d\nu_i} = \frac{1}{\Delta t_i} \qquad (2.9)\]

2.7.1. Вычисление частных производных

С использованием зависимостей (2.3) вычислим частные производные, входящие в выражения (2.7) и соответствующие для других сил.

Для силы \(F_i^c\):

\[\frac{\partial F_i^c}{\partial x_i} = 0,\quad \frac{\partial F_i^c}{\partial \nu_i} = 0,\quad \frac{\partial F_i^c}{\partial a_i} = 0,\]

поскольку \(F_i^c\) не зависит от перемещения, скорости и ускорения тела (2.10).

Для силы \(F_i^y\):

\[\frac{\partial F_i^y}{\partial x_i} = k,\quad \frac{\partial F_i^y}{\partial \nu_i} = 0,\quad \frac{\partial F_i^y}{\partial a_i} = 0 \qquad (2.11)\]

Для силы \(F_i^в\):

\[\frac{\partial F_i^в}{\partial x_i} = 0,\quad \frac{\partial F_i^в}{\partial \nu_i} = 2\mu |\nu_i|,\quad \frac{\partial F_i^в}{\partial a_i} = 0 \qquad (2.12)\]

Для силы \(F_i^u\):

\[\frac{\partial F_i^u}{\partial x_i} = 0,\quad \frac{\partial F_i^u}{\partial \nu_i} = 0,\quad \frac{\partial F_i^u}{\partial a_i} = m \qquad (2.13)\]

Подставляя найденные значения частных производных и коэффициенты из (2.9) в формулы (2.7) и аналогичные, получаем:

\[\frac{dF_i^c}{dz} = 0,\quad \frac{dF_i^y}{dz} = \frac{k \Delta t}{2},\quad \frac{dF_i^в}{dz} = 2\mu |\nu_i|,\quad \frac{dF_i^u}{dz} = \frac{m}{\Delta t_i} \qquad (2.14)\]

Суммируя слагаемые в формуле (2.6), приходим к выражению:

\[\frac{df(z)}{dz} = \frac{k \Delta t_i}{2} + 2\mu |\nu_i| + \frac{m}{\Delta t_i} \qquad (2.15)\]

Полученный результат структурно аналогичен ранее выведенной формуле (1.19), коэффициенты которой определяются соотношениями (1.12).

2.7.2. Преимущества автоматического подхода

Располагая возможностью вычислять \(f(z)\) и \(df(z)/dz\), продолжение расчета по описанному алгоритму не представляет затруднений. Какими же преимуществами окупаются дополнительные выкладки?

1. Отпадает необходимость в явном выписывании дифференциального уравнения движения.

2. Развернутая форма нелинейного уравнения вида (1.12) оказывается невостребованной.

3. Функциональная зависимость для производной \(df(z)/dz\) не требуется.

Рис. 2.1. Структурная схема формирования математической модели

Для формирования математической модели использована следующая информация (далее — «перечень»):

1. Сведения о стыковке элементов схемы (рис. 2.1).

2. Условие равновесия сил для \(i\)-го момента времени — нелинейное алгебраическое уравнение относительно \(x_i\), \(\nu_i\), \(a_i\):

   \[\sum_{k=1}^{4} F_i^{(k)} = 0 \qquad (2.16)\]

   (топологическое уравнение, определяемое структурой связей).

3. Компонентные уравнения — выражения, позволяющие определить усилия в каждом элементе как функции перемещения, скорости, ускорения и времени (2.3).

4. Частные производные усилий по перемещению, скорости, ускорению (2.10)–(2.13).

5. Алгебраические уравнения связи \(x\), \(\nu\), \(a\) — формулы метода интегрирования (2.4).

6. Выбор базисной переменной \(z\) (из \(x\), \(\nu\) или \(a\)) для решения нелинейного уравнения.

Данной информации достаточно для реализации машинных алгоритмов формирования математической модели объекта.

2.7.3. Модели элементов

Пусть имеется техническая система, процессы в которой подлежат анализу. Расчетная схема соответствует рис. 1.1.

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

Рис. 2.2. Типы элементов и их степени свободы

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

Для расчетной схемы (рис. 1.1) пользователь выбирает модели элементов из библиотеки и описывает структуру:

$ FRAGMENT: Пример
# BASE: 1
# STRUCTURE:
Пружина' K (1 2; Коэффициент жесткости)
Нелинейный демпфер ' MUNL (1 2; Коэффициент вязкости)
Масса ' M (2; Масса тела)
Воздействие ' FSIN (2 1; Q, T, начальная фаза)

Этим описанием пользователь сообщает программному комплексу PRADIS всю необходимую информацию по стыковке элементов схемы (пункт 1 перечня). Выбраны модели K, MUNL, M, FSIN; они соединены надлежащим образом; узел 1 описан как неподвижный.

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

Этап численного интегрирования представляет собой последовательность шагов по времени, каждый из которых сводится к решению нелинейного уравнения равновесия вида (2.16). Для этого необходима информация, приведенная в пунктах 3–6 перечня.

Входными данными для любой модели элемента являются:

• неизменный список параметров модели;

• текущие значения перемещений, скоростей и ускорений узлов, с которыми элемент соединен.

По этим данным модель элемента для текущего момента времени должна вычислить:

1. усилия, действующие со стороны системы на элементы (вектор усилий элемента) — пункт 3 перечня;

2. частные производные усилий по перемещениям, скоростям и ускорениям узлов элемента (матрицу Якоби элемента) — пункт 4 перечня.

Если элемент имеет \(N\) степеней свободы, длина вектора усилий элемента равна \(N\), а якобиан элемента имеет размерность \(N \times N \times 3\).

Рис. 2.3. Двухузловая модель упругой пружины

Например, разработчик двухузловой модели одномерной упругой пружины реализовал следующие зависимости (рис. 2.3):

\[F_1 = k(x_1 - x_2),\quad F_2 = k(x_2 - x_1) \qquad (2.17)\]

\[\frac{\partial F_1}{\partial x_1} = k,\; \frac{\partial F_1}{\partial x_2} = -k,\; \frac{\partial F_2}{\partial x_1} = -k,\; \frac{\partial F_2}{\partial x_2} = k \qquad (2.18)\]

\[\frac{\partial F_1}{\partial \nu_1} = \frac{\partial F_1}{\partial \nu_2} = \frac{\partial F_2}{\partial \nu_1} = \frac{\partial F_2}{\partial \nu_2} = 0 \qquad (2.19)\]

\[\frac{\partial F_1}{\partial a_1} = \frac{\partial F_1}{\partial a_2} = \frac{\partial F_2}{\partial a_1} = \frac{\partial F_2}{\partial a_2} = 0 \qquad (2.20)\]

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

Поскольку узел 1 пружины закреплен, востребуется только информация, связанная с незакрепленным вторым узлом:

\[F_2 = k(x_2 - x_1) \qquad (2.21)\]

\[\frac{\partial F_2}{\partial x_2} = k,\quad \frac{\partial F_2}{\partial \nu_2} = \frac{\partial F_2}{\partial a_2} = 0 \qquad (2.22)\]

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

2.8. Алгоритм выполнения шага интегрирования

В продолжение начатого интегрирования выполним третий шаг по времени с использованием формального алгоритма, базирующегося на формулах (2.5)–(2.15).

По результатам второго шага (\(t_2 = 1.438 \times 10^{-3}\) с):

\[x_2 = 6.12 \times 10^{-5}\ \text{м},\quad \nu_2 = 0.08509\ \text{м/с},\quad a_2 = 59.21\ \text{м/с}^2\]

Величина шага \(\Delta t_2 = 0.438 \times 10^{-3}\) с, локальная погрешность \(lp_2 = 0.000018\).

Рекомендуемое значение \(\Delta t_3\) определяем по формуле (1.30) с \(c = 0.8\) и \(\delta_l = 0.001\):

\[\Delta t_3 = 0.8 \cdot 0.438 \times 10^{-3} \cdot \sqrt{\frac{0.001}{0.000018}} = 2.63 \times 10^{-3}\ \text{с}\]

Действуем по схеме, представленной на рис. 2.4а–2.4в.

Рис. 2.4а. Алгоритм выполнения шага интегрирования (начало)

Рис. 2.4б. Алгоритм выполнения шага интегрирования (продолжение)

Рис. 2.4в. Алгоритм выполнения шага интегрирования (окончание)

Перед третьим шагом известны значения \(t_2\), \(x_2\), \(\nu_2\), \(a_2\), \(\Delta t_3\).

1. Определение коэффициентов приведения якобиана (2.9):

   \[\frac{dx_3}{d\nu_3} = \frac{\Delta t_3}{2} = \frac{2.63 \times 10^{-3}}{2} = 1.31 \times 10^{-3}\]
   \[\frac{d\nu_3}{d\nu_3} = 1,\quad \frac{da_3}{d\nu_3} = \frac{1}{\Delta t_3} = \frac{1}{2.63 \times 10^{-3}} = 380.2\]

2. Начальное приближение (явный прогноз по формуле (1.28)):

   \[\nu_3^0 = \nu_2 + a_2 \Delta t_3 = 0.08509 + 59.21 \cdot 2.63 \times 10^{-3} = 0.24081\]

   Значения \(x_3^0\) и \(a_3^0\), необходимые для расчета в моделях элементов:

   \[a_3^0 = a_2 = 59.21\]
   \[x_3^0 = x_2 + \nu_2 \Delta t_3 + a_3^0 \frac{\Delta t_3^2}{2} = 6.12 \times 10^{-5} + 0.08509 \cdot 2.63 \times 10^{-3} + 59.21 \cdot \frac{(2.63 \times 10^{-3})^2}{2} = 46.5 \times 10^{-5}\]

3. Первая итерация Ньютона (см. рис. 2.4в).

4. Обращение к моделям элементов. Вычисление вектора сил и якобиана каждого элемента по текущим значениям \(x_3^0\), \(\nu_3^0\), \(a_3^0\) для незакрепленного узла:

   Пружина:

   \[F^y = k x = 20000 \cdot 46.5 \times 10^{-5} = 9.3\]
   \[\frac{\partial F^y}{\partial x} = k = 20000,\quad \frac{\partial F^y}{\partial \nu} = \frac{\partial F^y}{\partial a} = 0\]

   Демпфер:

   \[F^в = \mu \nu |\nu| = 1000 \cdot 0.24081 \cdot |0.24081| = 58.0\]
   \[\frac{\partial F^в}{\partial x} = 0,\quad \frac{\partial F^в}{\partial \nu} = 2\mu |\nu| = 2 \cdot 1000 \cdot 0.24081 = 481.6,\quad \frac{\partial F^в}{\partial a} = 0\]

   Точечная масса:

   \[F^u = m a = 0.1 \cdot 59.21 = 5.9\]
   \[\frac{\partial F^u}{\partial x} = \frac{\partial F^u}{\partial \nu} = 0,\quad \frac{\partial F^u}{\partial a} = m = 0.1\]

   Прикладываемая сила:

   \[F^c = Q \sin \frac{2\pi}{T} t = -1000 \sin \frac{2\pi}{0.2\pi} (1.438 \times 10^{-3} + 2.63 \times 10^{-3}) = -40.6\]
   \[\frac{\partial F^c}{\partial x} = \frac{\partial F^c}{\partial \nu} = \frac{\partial F^c}{\partial a} = 0\]

5. Суммирование сил (значение \(f(z)\) на текущей итерации):

   \[\sum F = F^y + F^в + F^u + F^c = -40.6 + 9.3 + 58.0 + 5.9 = 32.6\]

6. Вычисление \(df(z)/dz\) по формулам (2.6), (2.7), (2.7а):

   \[\frac{dF^c}{dz} = 0 \cdot 1.31 \times 10^{-3} + 0 \cdot 1 + 0 \cdot 380.2 = 0\]
   \[\frac{dF^y}{dz} = 20000 \cdot 1.31 \times 10^{-3} + 0 \cdot 1 + 0 \cdot 380.2 = 26.2\]
   \[\frac{dF^в}{dz} = 0 \cdot 1.31 \times 10^{-3} + 481.6 \cdot 1 + 0 \cdot 380.2 = 481.6\]
   \[\frac{dF^u}{dz} = 0 \cdot 1.31 \times 10^{-3} + 0 \cdot 1 + 0.1 \cdot 380.2 = 38.0\]
   \[\frac{df(z)}{dz} = 0 + 26.2 + 481.6 + 38.0 = 545.8\]

7. Приращение \(\Delta z^1\):

   \[\Delta z^1 = -\frac{f(z^0)}{f'(z^0)} = -\frac{32.6}{545.8} = -0.05973\]

8. Очередное приближение:

   \[z^1 = z^0 + \Delta z^1 = 0.24081 - 0.05973 = 0.18108\]

9. Уточнение \(x_3^1\), \(\nu_3^1\), \(a_3^1\) по формулам (2.8):

   \[\nu_3^1 = z^1 = 0.18108\]
   \[a_3^1 = \frac{\nu_3^1 - \nu_2}{\Delta t_3} = \frac{0.18108 - 0.08509}{2.63 \times 10^{-3}} = 36.5\]
   \[x_3^1 = x_2 + \frac{\nu_2 + \nu_3^1}{2} \Delta t_3 = 6.12 \times 10^{-5} + \frac{0.08509 + 0.18108}{2} \cdot 2.63 \times 10^{-3} = 41.1 \times 10^{-5}\]

10. Проверка условий завершения итераций:

    \[|f(z^0)| = 32.6 > \delta_f = 0.1,\quad |\Delta z^1| = 0.05973 > \delta_z = 0.001\]

    Итерационный процесс следует продолжить.

11. Проверка лимита итераций:

    \[j = 1 < j_{\text{max}} = 5\]

12. Увеличение счетчика итераций: \(j = 2\).

13. Вторая итерация дает:

    \[f(z^1) = -3.6,\quad \Delta z^2 = -0.00858\]
    \[x_3^2 = 39.8 \times 10^{-5},\quad \nu_3^2 = 0.17250,\quad a_3^2 = 33.2\]

    Условия завершения не выполнены.

14. Третья итерация (последняя):

    \[|f(z^2)| = 0.07 < \delta_f,\quad |\Delta z^3| = 0.00018 < \delta_z\]
    \[x_3^3 = 39.8 \times 10^{-5},\quad \nu_3^3 = 0.17232,\quad a_3^3 = 33.2\]

15. Оценка локальной погрешности:

    \[lp_3 = \left| \frac{\nu_3^p - \nu_3^c}{2} \right| = \left| \frac{0.24081 - 0.17232}{2} \right| = 0.034\]

16. Корректировка правила выбора шага. Практика расчетов показывает, что формула (1.30) приемлема только при соотношении \(\delta_l/lp_i\) близком к единице. При значительных отклонениях рекомендуемое значение шага часто оказывается завышенным. Используется уточненное правило:

    \[\begin{split}\Delta t_{\text{рек}} = \begin{cases} c \cdot \Delta t_i \cdot \dfrac{\delta_l}{lp_i}, & \text{если } \dfrac{\delta_l}{lp_i} < 0.25 \\[12pt] c \cdot \Delta t_i \cdot \sqrt[4]{\dfrac{\delta_l}{lp_i}}, & \text{если } \dfrac{\delta_l}{lp_i} > 7 \\[12pt] c \cdot \Delta t_i \cdot \sqrt{\dfrac{\delta_l}{lp_i}}, & \text{если } 0.25 \leq \dfrac{\delta_l}{lp_i} \leq 7 \end{cases} \qquad (2.23)\end{split}\]

    Для рассматриваемого случая:

    \[\frac{\delta_l}{lp_3} = \frac{0.001}{0.034} = 0.03 < 0.25\]
    \[\Delta t_{\text{рек}} = 0.8 \cdot 2.63 \times 10^{-3} \cdot \frac{0.001}{0.034} = 0.061 \times 10^{-3}\ \text{с}\]

17. Повторный расчет с \(\Delta t_3 = 0.061 \times 10^{-3}\) с дает для момента времени \(t_3 = t_2 + \Delta t_3 = 1.499 \times 10^{-3}\) с:

    \[x_3 = 6.64 \times 10^{-5}\ \text{м},\quad \nu_3 = 0.08862\ \text{м/с},\quad a_3 = 58.1\ \text{м/с}^2\]

    Локальная погрешность находится в допустимых пределах. Рекомендуемое значение для следующего шага: \(\Delta t_{\text{рек}} = 0.264 \times 10^{-3}\) с.

Расчет на третьем шаге завершен.

Литература

1. Виттенбург Й. Динамика систем твердых тел. — М.: Мир, 1980.