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

Общий вид уравнения движения системы с одной степенью свободы и нелинейной характеристикой
\[
\ddot{x}=f(t, x, \dot{x}) .
\]

Решение можно начать с определения начального ускорения (в момент времени $t=0$ ) из уравнения (2.55), что дает
\[
\ddot{x}_{0}=f\left(0, x_{0}, \dot{x}_{0}\right) .
\]

Искомое решение уравнения (2.55) в любой последующий момент времени $t$ будем записывать в следующей символической форме:
\[
x=F(t) .
\]

На рис. 2.23 показан график движения, представляющий собой гладкую кривую в плоскости $x t$. Через $x_{1}, x_{2}, \ldots, x_{i}, \ldots$ обозначены значения перемещения $x$ в моменты времени $t_{1}, t_{2}, \ldots, t_{i}, \ldots$, отстоя-
Рис. 2.23

щие друг от друга на длину временных интервалов $\Delta t_{1}, \Delta t_{2}, \ldots$, $\Delta t_{i}, \ldots$. Эти временные интервалы обычно берутся постоянной длительности $\Delta t$, что не является обязательным правилом.

В подходе, который в дальнейшем будем называть как метод усреднения по ускорению, скорость $\dot{x}_{i}$ в момент времени определяется по приближенной формуле
\[
\dot{x}_{i}=\dot{x}_{i-1}+\frac{\ddot{x}_{i-1}+\bar{x}_{i}}{2} \Delta t_{i}
\]

где $\dot{x}_{i-1}$ — скорость на предыдущем временном интервале $t_{i-1}$. Данная формула, известная как правило трапеции, показывает, что этот подход в последние годы приобрел известность как метод Ньюмарка * с $\beta=1 / 4$ или метод с постолнным усреднением по ускорению. Авторы предпочитают называть его методом усреднения по ускорению. На каждом шаге по времени ускорение берем как среднее арифметическое от $\ddot{x}_{i-1}$ и $\ddot{x}_{i}$. Аналогично по правилу трапеции определяем и приближенное выражение для перемещения
\[
x_{i}=x_{i-1}+\frac{\dot{x}_{i-1}+\dot{x}_{i}}{2} \Delta t_{i},
\]

где скорость на каждом шаге по времени берем как среднее арифметическое от $\dot{x}_{i_{-1}}$ и $\dot{x}_{i}$. Подстановка выражения (2.58) в (2.59) дает
\[
x_{i}=x_{l-1}+\dot{x}_{i-1} \Delta t_{i}+\left(\ddot{x}_{i-1}+\ddot{x}_{i}\right)\left(\Delta t_{i}\right)^{2} / 4 .
\]

Для того чтобы определить погрешность, обусловленную непосредственным использованием выражения (2.60) для определения $x_{i}$, можно воспользоваться разложением в ряд Тейлора, что дает
\[
\begin{array}{c}
x_{i}=x_{i-1}+\Delta t_{i} \dot{x}_{i-1}+\frac{\left(\Delta t_{i}\right)^{2}}{2 !} \ddot{x}_{i-1}+\frac{\left(\Delta t_{i}\right)^{3}}{3 !} x_{i-1}^{3}+\cdots ; \\
\ddot{x}_{i}=\ddot{x}_{i-1}+\Delta t_{i} \ddot{x}_{i-1}+\frac{\left(\Delta t_{i}\right)^{2}}{2 !} x_{i-1}^{(4)}+\cdots
\end{array}
\]

Подставляя разложение (2.62) в выражение (2.60) и вычитая из результата разложение (2.61), найдем ошибку локального усечения рлда (или остаточный член)
\[
R_{x}=\frac{\left(\Delta t_{i}\right)^{3}}{12} \bar{x}_{i-1}+\frac{\left(\Delta t_{i}\right)^{4}}{12} x_{i-1}^{(4)}+\cdots=R_{x 1}+R_{x 2}+\cdots
\]

Первое слагаемое в этом остаточном члене ряда является наиболее существенным и называется главным членом ошибки локального усечения рлда.

В данном методе выражение (2.60) не использовалось непосредственно, а вместо него последовательно применялись выражения (2.58) и (2.59). Поскольку величина $\ddot{x}_{i}$ не известна заранее, решение необходимо отыскивать, прибегая на каждом шаге к итерациям;

ниже приведены рекуррентные соотношения для определения $j$-й итерации на $i$-м шаге:
\[
\begin{array}{c}
\left(\dot{x}_{i}\right)_{j}=A_{i-1}+\left(\ddot{x}_{i}\right)_{j-1} \Delta t_{i} / 2, j>1 ; \\
\left(x_{i}\right)_{j}=B_{i-1}+\left(\dot{x}_{i}\right)_{j} \Delta t_{i} / 2 ; \\
\left(\dot{x}_{i}\right)_{j}=f\left[t_{i},\left(x_{i}\right)_{j},\left(\dot{x}_{i}\right)_{j}\right]
\end{array}
\]

где
\[
\begin{array}{l}
A_{i-1}=\dot{x}_{i-1}+\ddot{x}_{i-1} \Delta t_{i} / 2 ; \\
B_{i-1}=x_{i-1}+\dot{x}_{i-1} \Delta t_{i} / 2 .
\end{array}
\]

Данный итерационный процесс является независящим от начальных условнй, поскольку требует использования специальной формулы для определения на каждом шаге по времени первого приближения. Определив $\ddot{x}_{0}$ из выражения (2.56), можно начать вычислять итерацию для первого шага, находя приближенное выражение для $\dot{x}_{1}$ с помощью экстраполлционной формуль Эйлера:
\[
\left(\dot{x}_{1}\right)_{1}=\dot{x}_{0}+\ddot{x}_{0} \Delta t_{1} .
\]

Тогда первые приближения для $x_{1}$ и $\ddot{x}_{1}$ находим соответственно по формулам (2.65) и (2.66). Все последующие итерации на первом шаге по времени состоят в повторном использовании формул (2.64), (2.65) и (2.66).

Для того чтобы начать итерационный процесс на $i$-м шаге по времени, можно снова воспользоваться формулой Эйлера и определить первое приближение для $\dot{x}_{i}$ :
\[
\left(\dot{x}_{i}\right)_{1}=\dot{x}_{i-1}+\ddot{x}_{i-1} \Delta t_{i} .
\]

В формулах (2.67) и (2.68) предполагается, что ускорение является постоянным внутри шага по времени. Подставляя $\dot{x}_{i}$ из последнего выражения в формулу (2.59), получим
\[
x_{i}=x_{i-1}+\dot{x}_{i-1} \Delta t_{i}+\ddot{x}_{i-1}\left(\Delta t_{i}\right)^{2} / 2 .
\]

Вычитая из выражения (в) разложение (2.61) в ряд Тейлора для $x_{1}$, получим остаточный член ряда
\[
R_{x}^{\prime}=-\frac{\left(\Delta t_{i}\right)^{3}}{6} \ddot{x}_{i-1}-\frac{\left(\Delta t_{i}\right)^{4}}{24} x_{i-1}^{(4)}-\cdots=R_{x 1}^{\prime}+R_{x 2}^{\prime}+\cdots
\]

В данном выражении абсолютная величина главного члена вдвое больше, чем в выражении (2.63). Для того чтобы повысить точность результатов, получаемых на первой итерации $i$-го шага, можно воспользоваться следующей несколько более сложной формулой, которая справедлива только при постоянном шаге по времени:
\[
\left(\dot{x}_{i}\right)_{1}=\dot{x}_{i-2}+2 \ddot{x}_{i-1} \Delta t .
\]

Это выражение охватывает два одинаковых шага по времени от $t_{i-2}$ до $t_{i}$ (см. рис. 2.23) и в нем используется значение ускорения средней точки в момент времени $t_{i-1}$. Подставляя выражение (2.69) в (2.59), найдем
\[
x_{i}=x_{i-1}+\left(\dot{x}_{i-2}+\dot{x}_{i-1}+2 \ddot{x}_{i-1} \Delta t\right) \Delta t / 2 .
\]

Скорость $\dot{x}_{i-2}$ можно разложить в ряд Тейлора:
\[
\dot{x}_{i_{-2}}=\dot{x}_{i-1}-\ddot{x}_{i-1} \Delta t+\ddot{x}_{i-1}(\Delta t)^{2} / 4-\cdots
\]

Подставляя ряд (е) в выражение (д) и вычитая из полученного результата выражение (2.61), найдем главный член ошибки локального усечения ряда
\[
R_{x}^{\prime \prime}=\frac{(\Delta t)^{3}}{12} \ddot{x}_{i-1},
\]

который совпадает [см. выражение (2.63)] с $R_{x_{1}}$.
Выражения (в) и (д) являются лвными (или открытыми) экстраполяционными формулами (или «предикпором»), с помощью которых приближенное значение $x_{i}$ в явном виде выражается через ранее найденные значения $x, \dot{x}$ и $\ddot{x}$. С другой стороны, выражение (2.60) называется нелвной (или скрытой) интерполяционной формулой (или «корректором»), которая позволяет находить более точные значения $x_{i}$, если найдено приближенное значение $\ddot{x}_{i}$. Метод усреднения по ускорению состоит в однократном использовании предиктора, после чего применяются итерации с корректором. Такой подход известен как метод предиктора-корректора.

Для решений итерационного типа требуется использовать некоторый критерий для остановки процесса или изменения шага, а также необходимо задавать предельное число выполняемых итераций. Наиболее удобным критерием сходимости процесса на $i$-м шаге является сравнение разности двух последовательных значений $x_{i}$ и главного члена в $R_{x}$. Однако вычисление производных более высокого порядка, чем второй, не совсем удобно (причина состоит в том, что сами по себе ряды Тейлора являются не очень хорошей экстраполяционной формулой). Более удобный критерий состоит в контролировании числа значащих цифр в $x_{i}$ следующим образом:
\[
\left|\left(x_{i}\right)_{j}-\left(x_{i}\right)_{j-1}\right|<\varepsilon_{x}\left|\left(x_{i}\right)_{j}\right|,
\]

где $\varepsilon_{x}$ — заданная малая величина. Например, желая получить приближенно точность до четырех значащих цифр, можно взять $\varepsilon_{x}=0,0001$. Именно с такой точностью получены решения числовых примеров в данном параграфе.

Пример 1. Рассмотрим уже знакомое читателю линейное уравнение движения системы с одной степенью свободы
\[
m \ddot{x}+c \dot{x}+k x=Q(t) .
\]

Здесь задано: $m=1,79 \cdot 10^{2} \mathrm{H} \cdot \mathrm{c}^{2} / \mathrm{m} ; c=2,15 \cdot 10^{2} \mathrm{H} \cdot \mathrm{c} / \mathrm{m}, \quad k=1,61 \mathrm{H} / \mathrm{m}, Q(t)=$ $=Q_{1}=40,1 \mathrm{H}$ (ступенчатая функция). Тогда имеем $p=\sqrt{\frac{1,}{k / m}}=3 \mathrm{c}^{-1}, n=$ $=c /(2 \mathrm{~m})=0,6 \mathrm{c}^{-1}, \gamma=n / p=0,2$ и уравнение (з) принимает вид
\[
\ddot{x}+1,2 \dot{x}+9 x=9 .
\]

Точное выражение для динамических перемещений системы при демпфировании и действии возмущающей силы в виде ступенчатой функции $Q_{1}$ известно:
\[
\begin{array}{c}
x=\frac{Q_{1}}{k}\left[1-e^{-n t}\left(\cos p_{\text {д }} t+\frac{n}{p_{\text {Д }}} \sin p_{\text {д }} t\right)\right]=\frac{Q_{1}}{k}[1- \\
\left.-A e^{-n t} \cos \left(p_{\text {д }} t-\alpha_{\text {Д }}\right)\right]
\end{array}
\]

Рис. 2.24
где $\quad p_{\text {д }}=p \sqrt{1-\gamma^{2}}=3 \sqrt{0,96} \mathrm{c}^{-1} ; \quad A=\sqrt{1+\left(n / p_{\text {д }}{ }^{2}\right.}=3 \sqrt{2} / 4 ; \quad \alpha_{\text {д }}=$ $=\operatorname{arctg}\left(n / p_{\text {j }}\right)=\operatorname{arctg}(0,2 / \sqrt{0,96})$. На рис. 2.24 представлен вид точного решения для данного случая.

Готовясь применять численный способ решения, запишем уравнение (и) в форме $(2.55)$ :
\[
\ddot{x}=9-9 x-1,2 \dot{x} .
\]

Если начальные условия взять в виде $x_{0}=0$ и $\dot{x}_{0}=0$, то начальное значение ускорения [см. уравнение (2.56)]
\[
\ddot{x}=9 .
\]

Будем использовать постоянный шаг по времени $\Delta t=0,1$ с и вести вычисления с точностью для перемещения до четырех значащих цифр [см. неравенство (2.70)].

Peшение. Первые приближения на первом шаге вычисляем по следующим формулам:
\[
\begin{array}{l}
(2.67):\left(\dot{x_{1}}\right)_{1}=0+9 \cdot 0,1=0,9 ; \\
(2.65):\left(x_{1}\right)_{1}=0+0,9 \cdot 0,1 / 2=0,045 ; \\
(2.66):\left(\ddot{x_{1}}\right)_{1}=9-9 \cdot 0,045-1,2 \cdot 0,9=7,515 .
\end{array}
\]

Второе приближение дает
\[
\begin{array}{l}
(2.64):\left(\dot{x_{1}}\right)_{2}=0,45+7,515 \cdot 0,1 / 2=0,82575 ; \\
(2.65):\left(x_{1}\right)_{2}=0+0,82575 \cdot 0,1 / 2=0,041288 ; \\
(2.66):\left(\ddot{x_{1}}\right)_{2}=9-9 \cdot 0,041288-1,2 \cdot 0,82575=7,6375 .
\end{array}
\]

Третье приближение дает
\[
\begin{array}{l}
(2.64):\left(\dot{x}_{1}\right)_{3}=0,45+7,6375 \cdot 0,1 / 2=0,83188 ; \\
(2.65):\left(x_{1}\right)_{3}=0+0,83188 \cdot 0,1 / 2=0,041594 ; \\
(2.66):\left(\ddot{x}_{1}\right)_{3}=9-9 \cdot 0,041594-1,2 \cdot 0,83188=7,6274 .
\end{array}
\]

Четвертое приближение дает
\[
\begin{array}{l}
(2.64): \dot{x_{1}}=0,45+7,6274 \cdot 0,1 / 2=0,83137 \\
(2.65) .: x_{1}=0+0,83137 \cdot 0,1 / 2=0,041569 \\
(2.66): \ddot{x_{1}}=9-9 \cdot 0,041569-1,2 \cdot 0,83137=7,6282 .
\end{array}
\]

Пятое приближение дает
\[
\begin{array}{l}
(2.64): \dot{x}_{1}=0,45+0,76282 \cdot 0,1 / 2=0,83141 ; \\
(2.65): x_{1}=0+0,83141 \cdot 0,1 / 2=0,041570 ; \\
(2.66): \ddot{x}_{1}=9-9 \cdot 0,041570-1,2 \cdot 0,83141=7,6282 .
\end{array}
\]

Видно, что решение для первого шага получено с точностью до четырех значащих цифр. За пять циклов итераций первое приближение для второго шага определяем следующим образом:
\[
\begin{array}{l}
(2.64):\left(x_{2}\right)_{1}=0+2 \cdot 7,6282 \cdot 0,1=1,5256 \\
(2.65):\left(x_{2}\right)_{1}=0,08314+1,5256 \cdot 0,1 / 2=0,15942 ; \\
(2.66):\left(\ddot{x_{2}}\right)_{1}=9-9 \cdot 0,15942-1,2 \cdot 1,5256=5,7345 .
\end{array}
\]

Второе приближение дает
\[
\begin{array}{l}
(2.64):\left(\dot{x}_{2}\right)_{2}=1,2128+5,7345 \cdot 0,1 / 2=1,4995 ; \\
(2.65):\left(x_{2}\right)_{2}=0,08314+1,4995 \cdot 0,1 / 2=0,15812 ; \\
(2.66):\left(\ddot{x}_{2}\right)_{2}=9-9 \cdot 0,15812-1,2 \cdot 1,4995=5,7775 .
\end{array}
\]

Третье приближение дает
\[
\begin{array}{l}
(2.64):\left(\dot{x}_{2}\right)_{3}=1,2128+5,7775 \cdot 0,1 / 2=1,5017 \\
(2.65):\left(x_{2}\right)_{3}=0,08314+1,5017 \cdot 0,1 / 2=0,15823 ; \\
(2.66):\left(\ddot{x}_{2}\right)_{3}=9-9 \cdot 0,15823-1,2 \cdot 1,5017=5,7739 .
\end{array}
\]

Четвертое приближение дает
\[
\begin{array}{l}
(2.64):\left(\dot{x}_{2}\right)_{4}=1,2128+5,7739 \cdot 0,1 / 2=1,5015 ; \\
(2.65):\left(x_{2}\right)_{4}=0,08314+1,5015 \cdot 0,1 / 2=0,15822 ; \\
(2.66):\left(\ddot{x_{2}}\right)_{4}=9-9 \cdot 0,15822-1,2 \cdot 1,5015=5,7742 .
\end{array}
\]

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

Из табл. 2.1а, где приведены результаты для 20 шагов по времени, можно видеть, что приближенные значения $x$ совпадают с точными [полученными с помощью выражения (к)] вплоть до трех значащих цифр. Таким образом, показанные на рис. 2.24 точки практически совпадают с соответствующими точками на кривой, представляющей точное решение.
2.1a. Решения примера 1 методом осреднения по ускорениям
Другой способ получения приближенного решения уравнения (2.55) известен как метод линейного ускорения. Как следует из

eго названия, в этом способе принято предположение, что ускорение изменяется по линейному закону и длине шага по времени. В соответствии с этим допущением выражение для $\ddot{x}$ на шаге по времени $\Delta t_{i}$ (см. рис. 2.23) можно записать в виде
\[
\ddot{x}\left(t^{\prime}\right)=\ddot{x}_{i-1}+\left(\ddot{x}_{i}-\ddot{x}_{i-1}\right) t^{\prime} / \Delta t_{i},
\]

где время $t^{\prime}$ отсчитывается от начала шага. Если ускорение изменяется по линейному закону, то соответствующая_ему скорость будет изменяться во времени по параболическому закону, а перемещение — по кубическому, и тогда имеем
\[
\begin{array}{c}
\dot{x}\left(t^{\prime}\right)=\dot{x}_{i-1}+\ddot{x}_{i-1} t^{\prime}+\left(\ddot{x}_{i}-\ddot{x}_{i-1}\right)\left(t^{\prime}\right)^{2} / 2 \Delta t_{i} ; \\
x\left(t^{\prime}\right)=x_{i-1}+\dot{x}_{i-1} t^{\prime}+\ddot{x}_{i-1}\left(t^{\prime}\right)^{2} / 2+\left(\ddot{x}_{i}-\ddot{x}_{i-1}\right)\left(t^{\prime}\right)^{3} /\left(6 \Delta t_{i}\right) .
\end{array}
\]

В конце шага по времени скорость и перемещение имеют вид
\[
\begin{array}{c}
\dot{x}_{i}=\dot{x}_{i-1}+\left(\ddot{x}_{i-1}+\ddot{x}_{i}\right) \Delta t_{i} / 2 \\
x_{i}=x_{i-1}+\dot{x}_{i-1} \Delta t_{i}+\left(2 \ddot{x}_{i-1}+\ddot{x}_{i}\right)\left(\Delta t_{i}\right)^{2} / 6 .
\end{array}
\]

Выражение (2.71) совпадает с выражением (2.58) метода усреднения по ускорению, а выражение (2.72) несколько отличается от соответствующего ему выражения (2.60). Если выражение (2.62) для $\tilde{x}_{i}$ подставить в выражение (2.72) и из последнего вычесть разложение в ряд Тейлора для $x_{i}$ [выражение (2.61)], то получим главный член остатка для этого случая:
\[
R_{x 1}^{\prime \prime \prime}=-\frac{\left(\Delta t_{i}\right)^{4}}{24} x_{i-1}^{(4)} .
\]

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

Применим метод линейных ускорений в духе, аналогичном описанному выше для подхода с усредненными ускорениями. Поскольку выражение (2.71) аналогично (2.58), то и рекуррентное выражение $j$-й итерации для $\dot{x}_{i}$ совпадает с выражением (2.64). Для получения соотношения, непосредственно связывающего $x_{i}$ и $\dot{x}_{i}$, найдем $\ddot{x}_{i}$ из выражения (2.71) и подставим его в (2.72). Тогда получим
\[
x_{i}=x_{i-1}+\left(2 \dot{x}_{i-1}+\dot{x}_{i}\right) \Delta t_{i} / 3+\ddot{x}_{i-1}\left(\Delta t_{i}\right)^{2} / 6 .
\]

В результате приходим к следующим рекуррентным формулам $j$-й итерации для $x_{i}$ :

где
\[
\left(x_{i}\right)_{j}=B_{i-1}^{*}+\left(x_{i}\right)_{j} \Delta t_{i} / 3,
\]
\[
B_{i-1}^{*}=x_{i-1}+2 \dot{x}_{i-1} \Delta t_{i} / 3+\bar{x}_{i-1}\left(\Delta t_{i}\right)^{2} / 6 .
\]

Полученные выше формулы (2.67), а также (2.68) или (2.69) можно вновь использовать для получения на каждом шаге первоначальных значений в итерационном процессе.

Когда приведенную в примере 1 задачу решали методом линейных ускорений, были получены результаты, приведенные в табл. 2.1б. В данном случае большинство приближенных значений $x$ ближе к точному решению, чем приведенные в табл. 2.1а, и были получены методом усредненных ускорений. Проверив эти методы на линейной задаче, применим их теперь для исследования примеров нелинейных задач.
2.16. Решения примера 1 методом линейных ускорений
Пример 2. Движение простого маятника, показанного на рис. 2.3 (см. п. 2.1), описывается нелинейным уравнением (2.4a)
\[
\ddot{\varphi}+p^{2} \cdot \sin \varphi=0,
\]

где $p^{2}=g / L$. Если положить длину $L$ численно равной ускорению $g$, получим $p^{2}=1$. Начальные условия возьмем в виде $\varphi_{0}=\pi / 2, \varphi_{0}=0$. В п. 2.2 было получено точное выражение (2.12) для периода колебания маятника. Используя начальное условие вида $\varphi_{0}=\varphi_{\mathrm{M}}=\pi / 2$ по таблицам эллиптических интегралов, найдем $k=F(k, \pi / 2) \rightleftharpoons 1,8541$. Отсюда следует, что четверть периода $\tau / 4 \doteq k / p=$ $=1,8541$ с. С другой стороны, если взять $p=k=1,8541$, то получим $\tau / 4=1 \mathrm{c}$. Именно это значение в силу его простоты и будет использовано ниже.
Таким образом, уравнение, которое надо решать численно, имеет вид
\[
\ddot{\varphi}=-p^{2} \sin \varphi=-3,437687 ;
\]

при этом начальные условия таковы:
\[
\ddot{\varphi}_{0}=-p^{2} \sin \pi / 2=-3,437687 .
\]

В табл. 2.2 приведены результаты для 20 шагов по времени (с шагом $\Delta t=$ $=0,1$ c) как методом усредненных ускорений, так и методом линейных ускорений. Величина угла $\varphi$ в момент времени $t_{10}$ должна равняться нулю, и метод линейных ускорений дает меньшее из получаемых обоими методами приближенное значение $\varphi_{10}$. Однако оба метода дают правильное конечное значение угла $\dot{\varphi}_{20}=-1,5708$ рад. На рис. 2.25 показан график приближенных значений угла $\varphi$ в зависимости от времени.

Пример 3. В качестве второго примера с нелинейной задачей рассмотрим следующее уравнение движения системы, в которой имеется пружина с возрастающей жесткостью:

или
\[
m \ddot{x}+c \dot{x}+k\left(x+\alpha x^{3}\right)=Q(t)
\]
\[
\ddot{x}+2 n \dot{x}+p^{2}\left(x+\alpha x^{3}\right)=q(t) .
\]

Из задачи 2.21 (см. п. 2.2) возьмем следующие данные: $\ddot{m}=1,79 \mathrm{H} \cdot \mathrm{m} / \mathrm{c}, k=$ $=715 \mathrm{H} / \mathrm{m}, \quad p^{2}=4 \mathrm{c}, \alpha=1,55 \cdot 10 \mathrm{~m}, c=Q(t)=0$. Для указанных значений параметров уравнение, которое необходимо решать численно, имеет вид
\[
\ddot{x}=-p^{2}\left(x+\alpha x^{3}\right)=-4\left(x+2 x^{2}\right) .
\]

Подставляя начальные условия $x_{0}=0$ и $\dot{x}_{0}=0$ в уравнение (ц), получим
\[
\ddot{x}_{0}=-4(0+0)=0 .
\]

Результаты, полученные обоимн методами с использованием 20 шагов по времени с шагом $\Delta t=0,025 \mathrm{c}$, приведены в табл. 2.3. На рис. 2.26 показан график приближенных значений перемещения $x$ в зависимости от времени. Максимальное значение, приближенно равное $0,05 \mathrm{~m}$, появляется, как это и должно быть, вблизи момента времени $t_{12}=0,30 \mathrm{c}$.

Хотя примеры с итерационными расчетами из этого параграфа можңо решать с помощью калькуляторов, тем не менее, вычисли-
Рис. 2.25
Рис. 2.26

тельные операции становятся довольно трудоемкими, поэтому лучше воспользоваться ЭВМ. Для расчетов примеров 1,2 и 3 методом осредненных ускорений были использованы специальные программы для ЭВМ под названием соответственно AVAC1A, AVAC2A, AVAC3A. Программы были составлены на языке БЕЙСИК, тексты их приведены в приложении. Для использования метода линейных ускорений эти программы легко переделать в программы LINAC1A, LINAC2A, LINAC3A путем изменения всего нескольких строк в каждой из них. Кроме того, большинство задач, приведенных в конце этого параграфа, можно решить, слегка изменив эти же программы.

Если уравнение движения (2.55) является линейным, то при численном решении можно избежать использования итераций с неявными формулами. Прямая формула линейной экстраполяции для метода усредненных ускорений может быть получена подстановкой выражений (2.58) и (2.60) в уравнение (3), решив которое относительно $\ddot{x}_{i}$ получим
\[
\ddot{x}_{i}=Q_{i}^{*} / m^{*},
\]

где
\[
\begin{array}{c}
m^{*}=m+c \Delta t_{i} / 2+k\left(\Delta t_{i}\right)^{2} / 4 ; \\
Q_{i}^{*}=Q_{i}-c\left(\dot{x}_{i-1}+\ddot{x}_{i-1} \Delta t_{i} / 2\right)-k\left[x_{i-1}+\dot{x}_{i-1} \Delta t_{i}+\ddot{x}_{i-1}\left(\Delta t_{i}\right)^{2} / 4\right] .
\end{array}
\]

В качестве нагрузки $Q_{i}$ обычно берут его постоянное среднее значение на шаге $\Delta t_{i}$ по времени. Уравнение (2.75a) представляет собой явную формулу для ускорения $\ddot{x}_{i}$, выраженного через известные к началу шага значения $x_{i-1}, \dot{x}_{i-1}$ и $\ddot{x}_{i-1}$. Это выражение вместе с выражениями (2.58) и (2.59) позволяет вычислять $x_{i-1}, \dot{x}_{i-1}$ и $\ddot{x}_{i-1}$ на каждом шаге. Таким образом, имеем прием прямой экстраполяции, не зависящий от начального приближения и не требующий итераций.

Қак метод линейных ускорений, так и аналогичные ему методы можно также применять в сочетании с процедурой прямой экстраполяции. Более того, уравнение движения для нелинейной системы может быть линеаризовано на малом шаге по времени, будучи записанным через приращения следующим образом:
\[
m_{i-1} \Delta \ddot{x}_{i}+c_{i-1} \Delta \dot{x}_{i}+k_{i-1} \Delta x_{i}=\Delta Q_{i} .
\]

Здесь $\Delta x_{i}, \Delta \dot{x}_{i}, \Delta \ddot{x}_{i}$ и $\Delta Q_{i}$ — приращения соответственно перемещения, скорости, ускорения и нагрузки на $i$-м шаге; $m_{i-1}, c_{i-1}$ и $k_{i-1}$ — значения массы, постоянной демпфирования и жесткости в начале каждого шага. Формулы явной схемы, аналогичные выражению (2.75a), можно записать относительно приращений ускорения $\Delta \ddot{x}_{i}$, скорости $\Delta \dot{x}_{i}$ и \»перемещения $\Delta x_{i}^{*}$.

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

Более точные, чем описанные здесь, методы можно найти в литературе по численному анализу *. Наиболее часто используемые подходы основываются либо на разложении искомой функции $x=$ $=F(t)$ и ее производных в ряды Тейлора **, либо на использовании формул интегрирования для полиномиальных интерполирующих функций ***. Они ориентируются на дифференциальные уравнения первого порядка и отражают ту точку зрения, что любое выражение вида уравнения (2.55) с производными второго порядка можно представить в форме двух уравнений первого порядка. Последнюю форму получаем введением вспомогательной зависимой переменной
\[
\dot{x}=y .
\]

Подставляя переменную $y$ вместо $\dot{x}$ в уравнение (2.55), получаем
\[
\dot{y}=f(t, x, y) .
\]

Уравнения (2.77a) и (2.77б) представляют систему двух уравнений первого порядка, которую можно интегрировать численно, используя параллельные одинаковые выражения для экстраполяции неизвестных $x$ и $y$. Хотя этот прием и прост, здесь не удалось обойти тот факт, что уравнение (2.77a) имеет более специфическую форму, чем уравнение (2.77б). Более эффективный подход состоит в том, чтобы вместо использования параллельной одинаковой экстраполяции выразить $y$, а затем и $x$ в виде ряда. Если поступить согласно сказанному, дальнейшая методология сводится к тому же численному решению уравнения второго порядка с использованием двух последовательных экстраполирующих формул для $\dot{x}$ и $x$, как и описано в данном параграфе.

ЗАДАЧИ

2.6.1. Рассмотреть пример 1 , взяв вместо ступенчатой функции линейную функцию вида $Q(t)=9 t$. Все остальные искомые данные взять из примера 1 .
Oтвет: $x_{10} \approx 1,93 \cdot 10^{-2} \mathrm{M}$.
2.6.2. Для рассмотренной в примере 1 линейной системы взять $Q(t)=0$, а начальные условия принять в виде $x_{0}=2,54 \cdot 10^{-2} \mathrm{~m} ; \dot{x}_{0}=-2,54 \cdot 10^{-2} \mathrm{~m} / \mathrm{c}$. Bсе остальные данные взять те же, что и в примере 1 .
Ответ: $x_{10} \approx 1,4 \cdot 10^{-2} \mathrm{M}$.
2.6.3. Вновь рассмотреть пример 2 , изменив начальное условие: вместо $\varphi_{0}=0$ взять $\varphi_{0}=2,618$ рад/с. Все остальные данные взять те же, что и в указанном примере.
Ответ: $\varphi_{10} \approx \pi / 2$ рад.
2.6.4. Для системы, рассмотренной в задаче 2.1.1 (см. п. 2.1), с помощью метода усредненных ускорений определить динамические перемещения при свобод-
ных колебаниях, взяв в качестве начальных условий следующие: $\dot{x}_{0}=l / 4 ; \dot{x}_{0}=$ $=y_{0}=\dot{y}_{0}=0$. Использовать такое уравнение движения и следующие значения параметров системы: $k=1,79 \cdot 10^{3} \mathrm{H} / \mathrm{m} ; l=0,102 \mathrm{~m} ; m=1,43 \cdot 10^{3} \mathrm{H} \cdot \mathrm{c}^{2} / \mathrm{m}$. Взять постоянный шаг по времени $\Delta t=0,1$ с, число шагов принять равным десяти и построить график полученного решения.
Oтвет: $\dot{x}_{\max } \approx 4,06 \cdot 10^{-2} \mathrm{~m} / \mathrm{c}$.
2.6.5. Для системы, рассмотренной в задаче 2.1.6 (см. п. 2.1), с помощью метода линейных ускорений определить динамические перемещения при свободных колебаниях с начальными условиями вида: $\varphi_{0}=0, \varphi_{0}=10,64$ рад/с. Использовать такое уравнение движения и следующие значения параметров системы: $W=22,7 \mathrm{H}$; $l=0,254 \mathrm{~m} ; \quad k_{\mathrm{K}}=1,16 \cdot 10^{2} \mathrm{H} \cdot \mathrm{M} /$ рад. Взять постоянный шаг по времени $\Delta t=$ $=0,025 \mathrm{c}$, число шагов принять равным 20 , построить график полученного решения.
Omвeт: $\varphi_{\max } \approx \pi / 2$ рад.,
2.6.6. Вновь рассмотреть пример 3 , но характеристики пружины взять не из задачи 2.2.1, а из задачи 2.2.2. С помощью метода усредненных ускорений определить максимальное перемещение и время, когда оно возникает.
Ответ: $x_{\mathrm{M}} \approx 0,121 \mathrm{~m} ; t_{\mathrm{M}} \approx 0,725$ с.
2.6.7. Вновь рассмотреть пример 3 , но характеристики пружины взять не из задачи 2.2.1, а из задачи 2.2.4. С помощью метода линейных ускорений определить максимальное перемещение и время, когда оно возникает.
Опвет: $x_{\mathrm{M}} \approx 7,95 \cdot 10^{-2} \mathrm{M} ; t_{\mathrm{M}} \approx 0,50 \mathrm{c}$.
2.6.8. Решить задачу 2.5 .2 (см. п. 2.5) методом усредненных ускорений *, используя следующие значения параметров системы: $m=1,79 \cdot 10^{2} \mathrm{H} \cdot \mathrm{c}^{2} / \mathrm{m} ; k_{1}=$ $=1,79 \cdot 10^{2} \pi^{2} \mathrm{H} / \mathrm{M} ; \quad k_{2}=4 k_{1} ; \quad x_{1}=2,54 \cdot 10^{-2} \mathrm{M}$.
Omвem: $x_{\mathrm{M}} \approx 3,81 \cdot 10^{-2} \mathrm{M} ; t_{\mathrm{M}}=0,75 \mathrm{c}$.
2.6.9. Решить задачу 2.5.4 (см. п. 2.5) методом линейных ускорений *, используя следующие значения параметров системы: $m=1,79 \cdot 10^{2} \mathrm{H} \cdot \mathrm{c}^{2} / \mathrm{m} ; k=1,79 \times$ $\times 10^{2} \pi^{2} \mathrm{H} / \mathrm{M} ; x_{1}=2,54 \cdot 10^{-2} \mathrm{M}$.
Ответ: $x_{\mathrm{M}}=3,81 \cdot 10^{-2} \mathrm{~m} ; t_{\mathrm{M}} \approx 1,22 \mathrm{c}$.
2.6.10. Решить задачу 2.5.11 (см. п. 2.5) методом линейных ускорений *, используя следующие значения параметров системы: $m=1,79 \cdot 10^{2} \mathrm{H} \cdot \mathrm{c}^{2} / \mathrm{m} ; k=$ $=1,79 \cdot 10^{2} \pi^{2} \mathrm{H} / \mathrm{M} ; x_{1}=2,54 \cdot 10^{-2} \mathrm{M}$.
Omвeт: $x_{\mathrm{M}}=8,23 \cdot 10^{-2} \mathrm{~m} ; t_{\mathrm{M}} \approx 1,62 \mathrm{c}$.
* Задачи $2.6 .8,2.6 .9$ и 2.6.10 составлены для систем с кусочно-линейными характеристиками восстанавливающей силы, и на каждом этапе исследования величина $k$ изменяется скачком от единого постоянного значения до другого.