Математическое моделирование механических систем со многими степенями свободы Текст научной статьи по специальности «Физика»

УДК 681.5:681.3:519.6:531.1:517.91

Математическое моделирование механических систем со многими степенями свободы

© Ю.В. Журавлёв МГТУ им. Н.Э. Баумана, Москва, 105005, Россия

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

Ключевые слова: математическое моделирование, парашютная система, уравнения Лагранжа, численное интегрирование и дифференцирование.

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

Исходя из лагранжевых принципов теоретической механики, изучение систем с большим числом степеней свободы сводится к исследованию систем дифференциальных уравнений Лагранжа 2-го рода [1-4]. Неоднозначность выбора обобщенных координат позволяет изучать динамику систем в различных конфигурационных пространствах (фазовых пространствах, пространствах состояний — в другой терминологии). При использовании аппарата лагранжевой механики проблема моделирования и соответствующей адекватности решается существенно проще. Первоначальной задачей теоретической части исследования после выбора обобщенных координат является получение формулы кинетического потенциала (разности между кинетической и потенциальной энергиями механической системы), а итогом — формулы обобщенных сил (коэффициентов влияния вариаций обобщенных коорди-

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

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

независимых обобщенных координат ч = (, q2, ..., дп) . Движение системы с п степенями свободы описывается дифференциальным уравне-

Ш дТ дТ ди _ ■ {■ ■ - V нием Лагранжа 2-го рода — ----— - — + Ц, где qq2,..., <1п) —

ш дс[ дq дq

обобщенные скорости, Ц - (, Ц2, ..., Цп) — обобщенные силы непотенциального происхождения, и — силовая функция потенциальной природы, Т — кинетическая энергия системы. Будем считать известными Т - Т(ч,х), и - и(х), Ц - Ц(q,q,х), при этом

(

где А - А -

а

V ап1

Т - 2 <q, Л[ > + <Ь, >+Т0,

а

1 п

^ Г Ь1 ^

а

Ь -

пп у

V Ьп J

ач - аЧ (Ч Х) Ьг - Ьг (^ Х),

То- То (q, Х)< q, л[ > - 2 ал ч], <Ь, q > - 2 ьг[[г .

и ]-1

Путем введения в рассмотрение вектора обобщенных импульсов

дТ дТ

Р-—, Р-((, ..., Р„) , Рк(к -1,2, ..., п) функции Ла-

дq дЧк

гранжа, или иначе кинетического потенциала, Ь = Т + и, и градиента

дЬ дТ ди

кинетического потенциала —---1--, перейдем от одного урав-

дq дq дq

нения Лагранжа 2-го рода к Гамильтоновой системе двух уравнений

в канонической форме Якоби [1]: — - — + Ц, — - А_1(р - Ь), или , ШХ дq ШХ

$ -1 (У,Х),

аХ

где у =

{ рЛ ( / Л дТ

/ = ; ; / = — + б; Л = А"1 (р - я). Так как

IЛ) дя

кинетического потенциала — при вычислении вектора /, для чего

дТ дТ , ,

р = — = — = Ая + о, то вместо пары (я,я) можно оперировать экви-

дя дя

валентной парой (р, я) и говорить о задаче Коши для гамильтоновой

системы — = /(у,^), у= у0, ^ > у е Я2". Для решения последней

т

задачи Коши необходимо выбрать алгоритмы интегрирования [5-10]. Дополнительно укажем на необходимость вычисления градиента

дТ дя

потребуются регуляризующие алгоритмы [11-17]. Отдельного внимания заслуживает вопрос о вычислении вектора / из матричного уравнения А/11 = р - я, что должно опираться на метод квадратного корня [18]. В данной статье проблемы алгоритмизации моделирования многомерных механических систем с голономными связями изучены на примере математического моделирования системы груз— парашют (СГП).

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

Аналитический этап формирования модели пространственного движения системы груз—парашют. СГП является примером голо-номной системы. Груз сцеплен с коушем парашюта посредством идеального сферического шарнира — вертлюга. Груз и парашют приняты за абсолютно твердые осесимметричные тела вращения. Движение СГП рассмотрено в спокойной земной атмосфере в режиме потенциального обтекания при наполненном куполе парашюта. Груз именуется телом 1, парашют — телом 2. В обозначениях индексированных величин соответствующая величина относится к телу, номер которого указан индексом. На СГП действуют внешние силы потенциальной и непотенциальной природы. Потенциальные силы определяются плоскопараллельной моделью гравитационного поля Земли с постоянной напряженностью величины g0, и вектор напряженности направлен вниз. К непотенциальным относятся аэродинамические силы, для определения которых использована гипотеза о стационарном обтекании твердого тела в спокойной земной атмосфере. Аэродинамические

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

На рис. 1 показана СГП и абсолютная система О X У 2 .

^ ё ё ё ё

Рис. 1. Конфигурация СГП в абсолютной системе координат

Движение СГП исследуется в правой ортогональной земной системе координат О X У 2 при неподвижном полюсе О относительно

ё ё ё ё ё

Земли, причем ось У' направлена вверх по земной вертикали. С достаточной точностью система О X У 2 полагается абсолютной системой

ё ё ё ё

отсчета. С каждым телом (V = 1, 2) связана правая прямоугольная система осей ОХ У 2 , причем связанные системы имеют общий полюс

V V V' г

О в центре вертлюга. Ось ОХ,у совпадает с осью вращения v-го тела, и положительное направление оси ОХ2 ориентировано от купола пара-

шюта к коушу. Центры масс тел Су е ОХу имеют координаты хс1 = Ъ1 > 0, хс2 = Ъ2 < 0. Предполагаем, что расчет сил и моментов дает их проек-ци и на оси связанной системы координат. На рис. 1 показаны векторы: У0 — абсолютная скорость точки О; ш(у) — угловая скорость враще-

т? (V)

ния у-го тела; ; — главный вектор аэродинамических сил у-го тела, приложенный в точке Су; М(у) — главный аэродинамический момент, воздействующий со стороны среды на у-е тело относительно С . Поместим в точку О «подвижную» земную систему координат ОХУ2&, (оси Кёнига [4]), участвующую в поступательном перемещении относительно неподвижной земной системы О X У Z . Положение у-го

ё ё ё ё

тела относительно осей Кёнига определяется тремя эйлеровыми углами: рыскания уу, тангажа и крена фу. При этом угол уу — это угол между осью ОХу и продольной координатной плоскостью ОХУё; угол — угол между осью ОХу и горизонтальной координатной плоскостью ОХ2&; угол фу — угол между осью ОУу и продольной плоскостью ОХУСвязь между координатами вектора в связанной системе и в системе Кёнига задается матрицей направляющих косинусов:

Га 11 а12 а ^ 13

A = связ^зем а 21 а22 а23

Ча31 а32 а33 у

ап = сд-су, а12 = 5д, а13 = сд-5у, а21 = -сф - 5д - су + 5ф - 5у, а22 = сф - сд, а23 = сф - 5д - 5у + 5ф - су, а31 = 5ф- sd- су + сф - 5у, а32 = -5ф - сд, а33 = -5ф - sd- 5у + сф - су,

где 5у = sin у; су = cos у; sd = sin S; сд = cos д; 5ф = sin ф; сф = cos ф.

Так как положение связанных осей v-го тела характеризуется эйлеровыми углами у^ Sv и ф^ то на место у, д и ф в матрицу Асвяз^зем необходимо подставлять у^ Sv и фv соответственно.

Разложения векторов сил и скоростей могут приводиться к связанным осям

R (v) = R^i v + Rvj+ R<v)k v,

MM (v)= Mv)? v+ M;(v) j v + m ^k v,

c5 (v)=Q(5 v+ra(,v)Jv+®ív)k v,

где (1, ,ку) — ортонормированный правоориентированный репер у-й связанной системы осей ОХ У Z , либо к абсолютным осям

V V V'

к = я^ 1 ] ф к в

ы м = м£Ч я + ы( Ця + ы (

со м=ш(5,+ш(:)к,

где (1 г, ] г ,к^ — ортонормированный правоориентированным репер абсолютной системы координат 0ХУ2^.

Выбор обобщенных координат, определяющих пространственную конфигурацию СГП. СГП имеет 9 степеней свободы (6 у одного свободного тела и 3 у другого, как связанного по отношению поступательных перемещений). В качестве обобщенных координат примем 3 декартовых координаты точки О: х0, у0, г0, а также по 3 эйлеровых угла на каждое тело: и фу (V = 1, 2).

Выберем вектор обобщенных координат (, д4, д5, q7, д8, ) =

= (х0,у0,г0,у1,^1,ф1,у2,$2,ф2) и соответствующий вектор обобщенных

скоростей (1,с[2,с[з,с[4,с[5,с[6,С[7,сс8,с[9) = (х0,у0^0 ^^2^2,(Р2).

Необходимо получить аналитические выражения кинетического потенциала Ь = Т + и, а обобщенных сил Q1, Q2, ..., 09, как функций обобщенных координат q1, q2, ..., q9, обобщенных скоростей сс 1,сс г,...,^ 9 и времени X.

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

вращения тела и правую ортонормированную тройку ((^к^ связанной с телом системы осей ориентируем так, что орт 1 направлен из О по оси вращения и задает направление координатной оси ОХ, орт ] задает направление оси ОУ с привязкой к некоторому элементу тела, орт к задает ось 0Z. Пусть С — центр масс тела с известным положением гс' = ОС = Ы. Будем предполагать, что тензор инерции тела в связанных осях 0ХУ2 в точке О охарактеризован диагональной матрицей 00 = {©X,,0*|, что соответствует выбору координатных осей совмещенными с главными осями инерции тела.

Если У0 — абсолютная скорость полюса, С — абсолютная угловая скорость вращения тела, то кинетическая энергия тела массой т

2 (

T = -(mV2 + 2m(( хс5)rc' + cc -0О -ш).

Представим векторы в связанной системе координат

Vo =VX i +Vy j +VZ k, cc = c i j +шг k, r' = b i,

и заменим тензорное произведение на матричное шт0Ош, где ш = (шХ,шу,oz) ; 0О = diag{©X,©У,©*}; ©X,©У,©* — главные осевые моменты инерции тела в полюсе О, тогда

T = 2 ( + 2m (cz - Vo, с y )b + ©Хш2 +0 У шу +0* ш2).

Воспользуемся теоремой Штейнера, связывающей моменты инерции относительно параллельных осей ©X = Jx, ©У = Jy + mb2, ©* = Jz + mb2, где Jx, J Jz — главные осевые моменты инерции тела в центре масс, причем Jy = Jz.

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

T = 2(mVO + 2m(Cz -У0гшу)b + Jxc2x +(Jy + mb)( +ш2)).

Теперь за полюс О примем центр вертлюга СГП, а за связанные оси 0XYZ — координатные оси OXvYvZv v-го тела; mv — масса v-го тела,

r(v) r(v) j(v)í T-(v) r(v)\

Jx , Jy , Jy\Jy - Jz ¡ — главные осевые моменты инерции v-го тела

в точке Cv.

Заметим, что

Cxv = V vcSv + Ф v, Cyv = V v^A + & v SWv, Czv = "М7 vCWv c$v + ^ Cv,

(V ^ Oxv í (v) [ an (v) a 12 ( v)> a 13 (X ^ ло

V v Oyv = a (v) a 22 a 23 Уо

V V Ozv J a(v) V 31 a(v) 32 a(v) 33 J V z0 J

(V )

где а- — элементы матрицы поворота связанной системы координат v-го тела относительно подвижной земной системы координат ОХУХ^. При подстановке последних выражений в формулу кинетической энергии тела вращения, произвольно перемещающегося в пространстве, получим выражение кинетической энергии v-го тела СГП в следующем виде:

Т, =1 (41) Х2 + ^У2+03} ¿2 + 2 + а£$ 2+}Ф V +

+2а\1> хХо^к V + 2<) ^ у + 2«24; УоУ V + 2а2(>УоК +

+ 2а<4) ¿о\К V + 2аз((} ^ V + > V + 2а4б)у vФ V),

где

а (v) = а(у ) = а (v) = т "11 "22 "33 '"V'

^ = УV) К)2 + 2 (( + тЛ2 )(ФЛ)2

а(() = У (,v) + тЛ

7( V) =

г66

/v) = 14

,(V) -

г( V ) У

г( V )

2

.(V )

1 5

-(V) _

^21

6( V) 31

ас = mvbvк2V\

к(2Ч,) =

(V)

"а3 1)^>

а ^) 24 = тАк3"\

а (v) 2( = тКк4\

а (v) 34 = тАк("\

а (v) 3( = тКк6\

а (V) 4( = ((+т.

а (v) "46 =у: ^.

к^ = -( + а32) )сф сЗ„

к4")=а 22)сФ

(V)

"а3 2)^>

-(V) _

(

^2 3

33

^ = а 2V3)cФv

(V)

"а3 3) ^>

А2)

к

(V)

к

(V) _

= (,- CФv2CФv

Полная кинетическая энергия СГП, равная Т = ^ Т, примет вид

V=1

квадратичной формы Т = 1 цА¿¡. Матрица А — симметрическая и сильно разреженная:

А =

' а11 о о а14 а1( о а17 а18 о

а22 о а24 а2( о а27 а28 о

а33 а34 а3( о а37 а38 о

* а44 а4( а46 о о о

* а(( о о о о

* а66 о о о

* а77 а78 а79

а88 о

V а99 У

Ненулевые элементы матрицы следующие:

= £ а»\ ( = 1,6 2,

2

a =

V=1

4V) = a. ., (V = 1, 2), (i, j = 1,6);

iv = i при i = 1, 2, 3;

i = i + 3(v - 1) при i = 4, 5, 6. Обозначим ск = cos qk, sk = sin qk (k = 1,9).

a11 = a22 = a33 = m1 + m2,

av4v = Jv) ( )2 + 2 ( + m.b. ) C5V)

a = J(V) + m b2

5V5V t"'ví'V>

a = J (V)

U6V6V ^ X 5

a14V = mvbv^ K(V) =-(a21) +a3V) ) C5V ,

a15V = т.ь.к2 ), к2 )=a2 1^6.- a3 1)s6V,

a24V = mvbvк3 ), к3 ) =-(a22) +a32) )) c5v, a25V = mvbvк4 ), к4 ) = a2 2C6V - a3 2)S6V , a34V = m.b. к5ч,), ^ = - (a 23 + a3V) ) C5., a35V = mvbvк6 ), к6 ) = a2 3)c6V - a3 3)S6V , a4v5v = ( + m.bV), к™ = ( -C6v )c a = J (V)C

4V6V ^ X 5V •

Здесь

(V) = (V) = (V) =

а11 = С( С4 , а12 = , а13 = С( ,

а 2 1) = —С6., С4„ + 56„ , а2 2) = С6„ С(„ , а 2 3) = С6„ + Зб^, С4У ,

а3 1) = 36^, S(v С4У + С6.^ , а3 2) = — ^б,, С(г , а3 3) = — 36^, S(v S4V + С6.^ С4^! .

Силовая функция СГП. Силовая функция СГП и = -П, где П — потенциальная энергия СГП, равная сумме потенциальных энергий двух тел:

n=Z mv go Ус v •

Здесь ус — высота центра масс у-го тела в земной системе координат. Учитывая

Ус,, = Уо + A

(v )

О 1 ^связ^зем

Г ЪЛ о

v 0 у

где АЯ^зем = (а(/} ) = На?), найдем Ус, = Уо + Vv

s5 = sin ($v) • Итак, П = ( + m2) g0q2 + m1b1 sin q5 + m2b2 sin q8. Напомним, что Ъх > О, Ъ2 < 0.

Кинетический потенциал L = T + U = T - П можно, наконец, записать в виде

L = 21qAq - (m + m2) goq2- mibisin q5- mKsin q8.

Обобщенные силы непотенциальной природы, действующие на СГП со стороны внешней среды. Обобщенная сила Qk определяется как функция влияния вариации k-й обобщенной координаты 8qk на элементарную работу 5Wk = Qk8qk со стороны аэродинамических сил при отсутствии вариаций других обобщенных координат. Возде й-

т? (v)

ствие воздушной среды на v-е тело сводимо к главному вектору Rv ;

—(v)

и главному моменту M с известными разложениями в связанной v-й системе координат OXvYvZv. Пусть 5qx= 8x0 Ф 0 и одновременно 5qk =

2

= 0 V k Ф 1. Очевидно, что Q1 = ^ Rx v), где R(x v) — проекция вектора R(v)

v=1 g

на ось OXg. Так как R^) = R^Цj) + R(v^) + Rv)а(3), гдеR?},R(,v),Rv)— проекции вектора R(v) на связанные оси OXv , OYv , OZv тела, то

Q1 = ¿ ((<1° + Ryv)а(2) + Rv)а(3^. Аналогично найдем

v=1

Q2 = ¿((21 + R^ + RZv)а£), Q3 = + + R^?)•

v=1 v=1

Рассмотрим ненулевое рыскание первого тела 5q4 = 8^ Ф 0 при Sqk = 0 V k Ф 4. При этом векторная проекция силы R(1) на горизонталь-

v=1

ную земную плоскость совершает работу (R^^b sin ф1 - cos ф1) , а проекция момента M(1) на горизонтальную земную плоскость совершает работу cos ф1 - Mf} sin ф1

В итоге Q4 = (Mf - R^b1 )cos ф1 - ((b + Mf} )sin ф1. Аналогично получим

Q5 = (Mf + R(;]b1)cosф, + (M{;) - Rfb,)sinфр Q6 = M™,

Q7 = (M^2) - Rf»b2)cosф2 -(Rfb2 + MfVnф2, Q8 = (Mf + R^cos ф2 + (My2) - Rf^sin ф2, Q9 = Mf.

Таким образом, получены аналитические выражения кинетического потенциала L = T- П и обобщенных сил Q1,Q2,...,Qn как функций

обобщенных координат q1,q2,...,qn обобщенных скоростей qx,q2,...,qn и времени t. Аналитический этап математического моделирования динамической системы СГП можно считать завершенным.

Решение начальной задачи Коши для гамильтоновой системы. Наличие проблемы численного нахождения градиента кинетического потенциала заставило привлечь многошаговый метод численного интегрирования. Численное решение начальной задачи Коши — = f (y, t),

dt

строится в виде сеточной функции

{y^..^yn,Уn+l,...,yN} на равномерной сетке {t0,tx,...,tn,tn+x,...,tN},

где t0 < t1 <... < tn < tn+1 <... < tN < tf с шагом h = tn+1 - tn. Дифференциальному уравнению на n-м шаге соответствует равносильное инте-

tn+1

гральное y (tn+1) = y (tn) + J f (y (t), t) dt. Заменив подынтегральную

tn

вектор-функцию f (y (t), t) полиномиальным интерполянтом Pkn (t) степени k, построенным на (k + 1) точечном фронтовом отрезке данных

{- k, fn - k ), (tn - k+1, fn - k+1),..., (tn, fn Й, где fn = f (yn, О и Уп « y (tn), получим явную k-шаговую формулу Адамса:

ln+1

Уп+1 = Уп + { Pk,n (t) dt.

При k = 5 это явная 5-шаговая экстраполяционная формула Адамса—Башфорта

n

y-=y»+h (/»+2 V/» +v 2 /»+3 v f + fV 4 /» + HV ^

которая имеет локальную погрешность O'(h7), причем

Vfn = /» -/„-, (n = 1, 2, 3, ...), V2fn = f -2/»-j + f»-2, V3/» = /» -3/_i + 3/»-2 - /п_з, V4/» = /n - 4/»-i + 6/_2 - 4/»-3 + /„-4, V5/n = /» - 5/»-1 + 10/n-2 - 10/n-3 + 5/»-4 - /»-5.

Для расчета фронта (первых 6-ти узлов) сетки, t., i = 0, 1, 2, 3, 4, 5, можно использовать некоторый одношаговый метод.

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

d4k

Отдельно отметим двухточечную формулу центральных разностей. Пусть y = /(х) — функция и /'(0) — ее производная в точке х = 0. Пусть yk = f(kh\ к = -2, -1, 0, 1, 2 — отсчеты заданной функции.

Двухточечная формула центральных разностей у'0 = _ . Так как

2h

h2

/'(0) = у'0--< < х1? то методическая погрешность двух-

6

точечной формулы центральных разностей Аметод =— тах /"(£,)! •

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

Ц_<Ш ПС ш-15

лах будет округляться с погрешностью Двыч (х) < | (при p = 16 и t = 14). Поэтому

- р

t~ 1

х • 0,5-10"

д (у, = АВЬ1Ч (у,)+ АВЫЧ (у_,) <\у, | + \у_, |) 0 10-15 „ к 011 п-BbI4W 2h 2h 2h

\Уо

15

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

Из принципа равного вклада методической и вычислительной погрешностей Лметод = ЛВыч(>о) имеем

h = з^ =

' C

3 УоИ'

15

■¡«1 У" (?)|

1,44 • 10

-5

з з

\>y

max| y3)

с

Но из условия минимума (h* )2 C1 +—2 = min получим

h*

h* = зз|

f2C1

а,14 • 10

-5

W

max У(3)! Р0|

что на 20 % меньше шага, найденного из принципа равного вклада.

Почему двухточечная? А.Н. Тихонов предложил регуляризацию решения некорректных задач, к которым относится и численное дифференцирование [11]. С.Б. Стечкин решил экстремальную задачу [12]

sup

f"(*) ^m

f (*)- b (*)

где 5 — ограниченный однородный аддитивный оператор 5 е [С ^ С].

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

а =

R0 = R0f (х ) =

[25

f (х + a)-f (х-а)

5 >|| f 5 (х

2а

)-f (х

m >

(

Шаг в двухточечной формуле центральных разностей зависит от уровня вычислительной погрешности и от степени кривизны функции. Проблеме численного дифференцирования всегда уделялось особое внимание в отечественной и зарубежной математической периодике [13-17].

Заключение. Изложены теоретические основы методологии автоматизации моделирования механических систем со многими степенями свободы с использованием аппарата уравнений Лагранжа второго рода. Проведен аналитический этап построения математической модели СГП. Выделены такие ключевые этапы, как получение уравне-

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

ЛИТЕРАТУРА

[1] Гантмахер Ф.Р. Лекции по аналитической механике. Е.С. Пятницкий, ред. 3-е изд. Москва, Физматлит, 2005, 264 с.

[2] Лурье А.И. Аналитическая механика. Москва, Наука, 1961, 824 с.

[3] Дубошин Г.Н. Небесная механика. Аналитические и качественные методы. 2-е изд. Москва, Наука, 1978, 457 с.

[4] Голубев Ю.Ф. Основы теоретической механики, 2-е изд. Москва, Изд-во МГУ, 2000, 719 с.

[5] Shampine L.F., Gordon M.K. Computer solution of ordinary differential equations; the initial value problem. San Francisco, Freeman and Company, 1975, 318 p.

[6] Gear C.W. Numerical initial value problems in ordinary differential equations. Prentice-Hall, 1971, 253 p.

[7] Холл Дж., Уатт Дж. Современные численные методы решения обыкновенных дифференциальных уравнений. Москва, Мир, 1979, 312 с.

[8] Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. Москва, Наука, 1987, 600 с.

[9] R. van Wyk. Variable Mesh Multistep Methods for Ordinatory Differential Equations. J. Comput. Phys., 1970, vol. 5 (2), pp. 244-264.

[10] Хайрер Э., Нёрсетт С., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Нежесткие задачи. Москва, Мир, 1990, 512 с.

[11] Тихонов А.Н. О решении некорректно поставленных задач и методе регуляризации. ДАН СССР, 1963, т. 151, № 3, с. 501-504.

[12] Стечкин С.Б. Наилучшее приближение линейных операторов. Мат. заметки, 1967, т. 1, № 2, с. 137-148.

[13] Васин В.В. Об устойчивом вычислении производной в С(-да, да). ЖВММФ, 1973, т. 13, № 6, с. 1383-1389.

[14] Cullum J. SIAM J. Num. Anal., 1967, vol. 8 (2), pp. 254-262.

[15] Curtis A., Reid J.K. J. Inst. Math. Appl., 1974, vol. 13, pp. 121-126.

[16] Stepleman R.S., Winarsky N.D. Mathematics of Computation, 1979, vol. 33 (148), pp. 1257-1264.

[17] Dumontet J., Vignes J. Determination du pas Optimal Daus le Calcul des Derivees sur Ordinatoner R.A.I.R.O. Analyse Numerique. Numerical Analysis, 1977, vol. 11, № 1, pp. 13-25.

[18] Уилкинсон, Райнш. Справочник алгоритмов на языке АЛГОЛ. Линейная алгебра. Москва, Машиностроение, 1976, 600 c.

Статья поступила в редакцию 27.06.2013

Ссылку на эту статью просим оформлять следующим образом:

Журавлёв Ю.В. Математическое моделирование механических систем со многими степенями свободы. Инженерный журнал: наука и инновации, 2013, вып. 9. URL: http://engjournal.ru/catalog/mathmodel/technic/1117.html

Журавлёв Юрий Васильевич родился в 1947 г., окончил Московский авиационный институт в 1974 г., МГУ им. М.В. Ломоносова в 1981 г. Старший преподаватель кафедры «Вычислительная математика и математическая физика» МГТУ им. Н.Э. Баумана. Автор 20 научных публикаций в области моделирования, идентификации, управления динамическими системами. e-mail: zhurjurwas270747@yandex.ru