Проект компьютеризации расчета траектории неуправляемого летательного аппарата

Проект компьютеризации расчета траектории неуправляемого летательного аппарата

Таблица. Исходные данные и полученные результаты:

Таблица. Дополнительные параметры

Таблица. Движение по направляющим

Таблица. Активный участок траектории ЛА

Таблица. Пассивный участок траектории ЛА.

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

.1 Постановка задач. Анализ и представление исходных данных

компьютерный летательный аппарат

В этой главе рассмотрены следующие задачи:

Разработка компьютерной модели расчета траектории неуправляемого ЛА, используя методы численного интегрирования Эйлера и Рунге-Кутта.

Расчет траектории пассивного участка полета неуправляемого ЛА с использованием параболической теории.

Траектория полета неуправляемого ЛА (снаряда РСЗО) состоит из трех характерных участков:

движение по направляющим;

активный участок;

пассивный участок.

Расчет дополнительных параметров:

. Секундный массовый расход (q) рассчитан по формуле:

. Тяга двигателя (R) рассчитана по формуле:

=Ue*q

. Коэффициент лобового сопротивления рассчитан по формуле:=Cx(M)*i , где Cx=Cx(M)-коэффициент лобового сопротивления эталонного ЛА. Принимаем Cx =0.3.

. Площадь миделя ЛА рассчитана по формуле:

. Закон изменения массовой плотности воздуха от высоты полета принимаем в виде, где ro= 1.23- плотность воздуха на уровне моря; высота y в формуле берется в метрах.

. Закон изменения массы принимаем в виде m= m0-q*t.

1.2 Модель динамики неуправляемого ЛА

.2.1 Движение ЛА по направляющим

При расчете движения по направляющем принимается во внимание действие на ЛА следующих сил: R- тяга, G- вес ЛА, F- реакция направляющих, которая состоит из силы трения Fтр и нормальной составляющей N. Аэродинамическими силами при расчете движения по направляющим обычно пренебрегают.

Кроме того, в идеальном случае направляющие принимаются жесткими и недеформированными. Расчетная схема движения представлена на рисунке 1.1.

Рис 1.1.

Уравнение движения ЛА на направляющих в проекции на направление скорости запишется в виде

где f- коэффициент трения скольжения, m=m0-Q*t - принимаемый в работе закон изменения массы ЛА, G=m*g. Записывая уравнение движения в нормальной форме Коши, получаем

Последовательно интегрируя уравнение с учетом соотношения ds/dt=v, находим

Полагая в формулах t=td, получаем выражение для дульной скорости Vd и длины направляющих:

Если задана длина направляющих Sh, то, решая численно уравнение, определим td. Для этого уравнение запишем в виде

Если пренебречь изменением массы при движении ЛА по направляющим, то, заменяя массу ЛА средним значением m= mср,

получаем

По этим формулам можно приближенно найти vd и td.

Примечание: 1. Коэффициент трения скольжения ЛА об направляющие f равен 0.15 (f=0.15).

. Изменением массы при движении ЛА по направляющим пренебрегаем.

1.2.2 Активный участок полета ЛА

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

. угол атаки равен нулю.

. тяга R постоянна.

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

Здесь X=0.5*Cx*Ro*S*V - сила лобового сопротивления. Остальные элементы системы были указаны выше.

1.2.3 Пассивный участок полета ЛА

При расчете пассивного участка полета ЛА также используется система (1.1). При этом:

тяга R равна нулю (R=0).

поверхность Земли представляется в виде бесконечной плоскости.

масса ЛА постоянна (m= m-mt=const).

1.3 Алгоритмы численного интегрирования на основе методов Эйлера и Рунге-Кутта

Рассмотрим метод Эйлера - линейное приближение, использующее первые два члена ряда Тейлора. Здесь искомая интегральная кривая аппроксимируется ломаной линией. Если шаг h достаточно мал, то интеграл в формуле  можно вычислить используя теорему о среднем, т.е. вынося подынтегральную функцию из- под знака интеграла средним значением.

В методе Эйлера подынтегральная функция выносится при нижнем пределе интегрирования: . Это приближение геометрически соответствует движению от точки x к точке х+h по касательной к кривой y(x) в точке х. Запишем расчетные формулы метода Эйлера:

+1=yk+f(xk,yk)h, xk=xk-1+h(x0)=y0, yk=y(xk)

В усовершенствованном методе Эйлера- Коши в первом приближении полагается:

а во втором

Погрешность метода Эйлера определяется остаточным членом ряда Тейлора

т.е. R~h2 на каждом шаге вычислений. Для обеспечения сходимости шаг h следует выбирать достаточно малым. Для метода Эйлера- Коши погрешность имеет порядок h2.

Рассмотрим метод Рунге и Кутта. В основе получения вычислительных схем этого метода лежит разложение функции y(x) в ряд Тейлора с последующим преобразованием отрезка ряда к виду, не содержащему производных. На шаге h производная dy/dx=f(x,y) аппроксимируется параболой второго порядка. Здесь функция (x,h) определяется формулой парабол Симпсона (формула Ньютона - Котеса для трех узлов):

Рассмотрим дифференциальное уравнение при начальном условии (хА,уА). Выполним следующие операции:

По известным начальным условиям (хА,уА) определим значение производной в начальной точке А: .

и отметим значение ее ординаты в середине шага интегрирования h (точка В с координатами

).

Рис 1.2.

) Найдем значение производной по формуле  в точке В: и проведем из точки А прямую . Отметим значение ординаты этой прямой в середине шага интегрирования h (точка С с координатами ).

) По уравнению найдем значение производной в точке С: и проведем из точки А прямую . Отметим значение ординаты этой прямой в конце шага интегрирования h (точка D с координатами ).

) По уравнению  найдем значение производной в точке D: .

В результате построений найдем значение производных  в точках А, В,С и D. Отложим эти значения на графике рис 1.3. Как видно из графика, в точке с абсциссой  получены два значения производной вместо одного. Это следствие приближенности метода. Примем в этой точке среднее значение производной: . Отложив на графике (рис 1.3.) ординату , получим точку М.

Рис 1.3.

Будем считать, что кривая, изображающая зависимость  должна проходить через точки A,M и D. Проведем через эти три точки параболу, уравнение которой:

.

Значения коэффициентов a,b и с выбираются из условия прохождения параболы через точки А, М, и D. Коэффициент . Из уравнения параболы  имеем систему:

Решив эти уравнения, найдем:

Проинтегрируем теперь уравнение параболы  в пределах от x=xA до x=xA+h. Значение этого интеграла является приращением искомой функции y при изменении х на величину h. Таким образом  Подставив сюда полученные выше выражения для a,b,c, после приведения подобных членов для общего случая () получим:

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

Запишем расчетные формулы метода Рунге- Кутта:

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

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

1.4 Решение контрольной задачи для одного шага интегрирования методами Эйлера и Рунге-Кутта (Проведение расчетов без ПК)

.4.1 Начальные условия

Для решения тестового примера на 1 шаг интегрирования системы (1.1) методами Эйлера и Рунге- Кутта задаем следующие начальные условия:

.Шаг интегрирования h принимаем равным 0.05 секунды (h=0.05 сек.);

. V=Vg, где Vg- скорость ЛА в момент схода с направляющих, рассчитываемая по формуле:

, где

q_g - угол наклона направляющих; q_g =450=0,785 рад;дульное время, рассчитываемые по формуле:

=  и tg =  , где

==

. q=q_g - угол наклона направляющих; g=450=0,785 рад.

. х- начальная абсцисса, рассчитываемая по формуле:

. у- начальная ордината, рассчитываемая по формуле:

1.4.2 Расчет системы уравнений методом Эйлера (1 шаг)

1. Нахождение коэффициентов К1 для каждого уравнения.

По определению коэффициенты К1- это правые части уравнений, вычисленные при начальных условиях.

а. плотность при начальной высоте :

б. лобовое сопротивление при начальной (дульной) скорости:

=0,5*0.26*0.01767*1.22967*(72.9421)2=15.0288

в. Масса при начальном (дульном) времени:

=60-20*0.1095=57.81

г. К1(1)= (40000 - 9.81*57.81* 0.785 - 15.0288) / 57.81 = 683.96 м/с2

К1(2)=-9,81*cos450/ 22,708= -0.7761 рад/с

К1(3)=22,708*cos450= 3.51304 м/с

К1(4)=22,708*sin450= 3.4857 м/с

2. Результаты одного шага интегрирования по методу Эйлера:

1.4.3 Расчет системы уравнений методом Рунге-Кутта (1 шаг)

. Нахождение коэффициентов К1 для каждого уравнения.

По определению коэффициенты К1- это правые части уравнений, вычисленные при начальных условиях.

а. плотность при начальной высоте :

б. лобовое сопротивление при начальной (дульной) скорости:

=0,5*0.26*0.01767*1.22967*(72.9421)2=15.0288

в. Масса при начальном (дульном) времени:

=60-20*0.1095=57.81

г. К1(1)= 8,1224487289

К1(2)= -0,0112446838

К1(3)= 1,091383734

К1(4)= 1,0905149817

2. Нахождение коэффициентов К2 для каждого уравнения.

По определению коэффициенты К2- это правые части уравнений, вычисленные при полученных новых значениях параметров.

К2(1)= 8,1419713966

К2(2)= -0,009992507

К2(4)= 1,2270805529

. Нахождение коэффициентов К3 для каждого уравнения

По определению К3- это правые части уравнений, вычисленные при полученных новых значениях параметров.

К3(1)= 8,1417519898

К3(2)= -0,0099835329

К3(3)= 1,2415239848

К3(4)= 1,2282011286

. Нахождение коэффициента К4 для всех уравнений

По определению К4- это правые части уравнений, вычисленные в конце шага

К4(1)= 8,1610849832

К4(2)= -0,0089854278

К4(3)= 1,3930436092

К4(4)= 1,364415478

. Результаты одного шага по методу Рунге- Кутта.

1.4.4 Сравнение решений тестового примера на 1 шаг, полученных вручную (п.1.4.1 и 1.4.2) и с использованием программного комплекса

Сравнение результатов производится в соответствии с данными, представленными в таблице 1.3.:

Таблица 1.3

Примечание: Результаты компьютерного решения тестового примера взяты из результатов работы программного комплекса.

1.5 Результаты расчета траекторных параметров неуправляемого ЛА

.5.1 Движение по направляющим

Результатами расчета траекторных параметров данного участка траектории являются значения дульной скорости (vd) и дульного времени (td).

Результаты расчета приведены в таблице 1.4.:

Таблица 1.4

1.5.2 Активный участок траектории ЛА

Расчет траекторных параметров на конце активного участка полета ЛА произведен с помощью программного комплекса с шагом интегрирования равным 0,05

Таблица 1.5

1.5.3 Пассивный участок траектории ЛА

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

Таблица 1.6.

1.5.4 Результаты расчета траекторных параметров неуправляемого ЛА для характерных точек траектории

Результаты расчета траекторных параметров неуправляемого ЛА в данном пункте приводятся для характерных точек траектории, показанных на рисунке 1.11.

Рис 1.11.

где A- точка начала движения по направляющим;точка конца движения по направляющим (начало активного участка);точка конца активного участка (начало пассивного участка);наивысшая точка траектории;точка контакта снаряда с поверхностью;

Результаты расчета представлены в таблице 1.7.

Таблица 1.7

Примечание: Результаты расчета получены при помощи программного комплекса

1.6 Расчет траектории пассивного участка полета ЛА с использованием параболической теории

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

Основные формулы параболической теории:

Здесь t0=0, y0=yc=228,664 м, x0=xc=270,014 м,

-параметры конца активного участка рассчитанные по данным п 1.5.2

Ниже производится расчет пассивного участка траектории ЛА с использованием параболической теории для характерных точек параболической траектории.

Для вершины траектории (точка D) время вычисляется по формуле:

Подставляя tD в формулы 3 и 4, получим X(tD)=6626.495 м Y(tD)=2804.056 м

Для точки Е (конца траектории) имеем Y(tE)=0. Тогда

Результаты расчетов представлены в сравнении с расчетами реальной траектории в таблице 1.8.

Таблица 1.8

Координата Х, м (числ/ П.Т.)

Координата Y, м (числ/ П.Т.)

C

270,014/270,014

228,664/228,664

D

5520,904/6626.495

2475,650/2804.056

E

10534,668/19615.69

0/0

Приложение 1

Графики

1.      V(t) методами численного интегрирования и по параболической теории.

2.      θ(t) методами численного интегрирования и по параболической теории.

3.      y(x) (Метод Эйлера и параболическая теория).

4.      y(x) (Метод Рунге-Кутта и параболическая теория).

Рис.

Рис.

Рис.

Рис.

Рис.

Рис.

Приложение 2

Текст программы, разработанной для интегрирования методами Эйлера и Рунге-Кутта. (Программа выполнена на Visual Basic 6.0 Professional Edition)

Рис.

Base 1r0 = 1.23r As Doubley As Doublex As DoubleV As DoubleQ As DoubleXx As Doublem As DoubleYg As DoubleXg As DoubleCx As DoubleVg As DoubleS As DoubleMn As DoubleRashod As DoubleTg As DoubleRs As DoubleT As Doubleh As DoubleQg As Doublei As DoubleK(4, 4) As Doublestr As Stringj As DoubleSub Command1_Click().ForeColor = &HFF&.Caption = "идет расчет"= CDbl(Text4.Text) * Sin(CDbl(Text3.Text))= CDbl(Text4.Text) * Cos(CDbl(Text3.Text))= CDbl(Text5.Text)= CDbl(Text1.Text)= CDbl(Text6.Text)= CDbl(Text17.Text)= CDbl(Text18.Text)= CDbl(Text2.Text)= CDbl(Text19.Text)= CDbl(Text16.Text)= CDbl(Text7.Text)= CDbl(Text3.Text).Text = "".Text = "".Text = "".Text = "".Text = ""

'----------------------------------

'Активный участок траектории полета

'Эйлер= Xg

y = Yg

V = Vg

Q = Qg

j = Tg

For i = Tg To T Step h= r0 * Exp(-y / 7800)= 0.5 * Cx * r * S * V ^ 2= Mn - Rashod * i= V + h * (Rs - 9.81 * m * Sin(Q) - Xx) / m= Q - h * 9.81 * Cos(Q) / V= x + h * V * Cos(Q)= y + h * V * Sin(Q)Option1.Value = True And Round(i, 2) = Round(j, 2) Then= j + 1.Text = Text33.Text & Round(i, 2) & vbCrLf.Text = Text34.Text & Round(V, 1) & vbCrLf.Text = Text36.Text & Round(Q, 3) & vbCrLf.Text = Text40.Text & Round(x, 0) & vbCrLf.Text = Text38.Text & Round(y, 0) & vbCrLfIfi = Tg Then = "Значения коэфициентов" & vbCrLf & "на первом шаге Э" & vbCrLf

str = str & "V = " & Round(V, CLng(Text21.Text)) & vbCrLf= str & "Q = " & Round(Q, CLng(Text21.Text)) & vbCrLf= str & "X = " & Round(x, CLng(Text21.Text)) & vbCrLf= str & "Y = " & Round(y, CLng(Text21.Text)) & vbCrLf & vbCrLfIfi.Text = Round(V, CLng(Text21.Text)).Text = Round(Q, CLng(Text21.Text)).Text = Round(x, CLng(Text21.Text)).Text = Round(y, CLng(Text21.Text))

'-----------------------------

'Рунге-Кутт= Xg= Yg= Vg= Qg= Tgi = Tg To T Step h= r0 * Exp(-y / 7800)= 0.5 * Cx * r * S * V ^ 2= Mn - Rashod * i(1, 1) = h * (Rs - 9.81 * m * Sin(Q) - Xx) / m(1, 2) = -9.81 * h * Cos(Q) / V(1, 3) = h * V * Cos(Q)(1, 4) = h * V * Sin(Q)= r0 * Exp(-(y + 0.5 * K(1, 4)) / 7800)= 0.5 * Cx * r * S * (V + 0.5 * K(1, 1)) ^ 2= Mn - Rashod * (i + h / 2)(2, 1) = h * (Rs - 9.81 * m * Sin(Q + K(1, 2) / 2) - Xx) / m(2, 2) = -9.81 * h * Cos(Q + K(1, 2) / 2) / (V + K(1, 1) / 2)(2, 3) = h * (V + K(1, 1) / 2) * Cos(Q + K(1, 2) / 2)(2, 4) = h * (V + K(1, 1) / 2) * Sin(Q + K(1, 2) / 2)= r0 * Exp(-(y + 0.5 * K(2, 4)) / 7800)= 0.5 * Cx * r * S * (V + 0.5 * K(2, 1)) ^ 2= Mn - Rashod * (i + h / 2)(3, 1) = h * (Rs - 9.81 * m * Sin(Q + K(2, 2) / 2) - Xx) / m(3, 2) = -9.81 * h * Cos(Q + K(2, 2) / 2) / (V + K(2, 1) / 2)(3, 3) = h * (V + K(2, 1) / 2) * Cos(Q + K(2, 2) / 2)(3, 4) = h * (V + K(2, 1) / 2) * Sin(Q + K(2, 2) / 2)= r0 * Exp(-(y + K(3, 4)) / 7800)= 0.5 * Cx * r * S * (V + K(3, 1)) ^ 2= Mn - Rashod * (i + h)(4, 1) = h * (Rs - 9.81 * m * Sin(Q + K(3, 2)) - Xx) / m(4, 2) = -9.81 * h * Cos(Q + K(3, 2)) / (V + K(3, 1))(4, 3) = h * (V + K(3, 1)) * Cos(Q + K(3, 2))(4, 4) = h * (V + K(3, 1)) * Sin(Q + K(3, 2))= V + (K(1, 1) + K(4, 1) + 2 * (K(2, 1) + K(3, 1))) / 6= Q + (K(1, 2) + K(4, 2) + 2 * (K(2, 2) + K(3, 2))) / 6= x + (K(1, 3) + K(4, 3) + 2 * (K(2, 3) + K(3, 3))) / 6= y + (K(1, 4) + K(4, 4) + 2 * (K(2, 4) + K(3, 4))) / 6Option2.Value = True And Round(i, 2) = Round(j, 2) Then= j + 1.Text = Text33.Text & Round(i, 2) & vbCrLf.Text = Text34.Text & Round(V, 1) & vbCrLf.Text = Text36.Text & Round(Q, 3) & vbCrLf.Text = Text40.Text & Round(x, 0) & vbCrLf.Text = Text38.Text & Round(y, 0) & vbCrLfIfi = Tg Then

str = str & "Значения коэффициентов" & vbCrLf & "на первом шаге Р-К" & vbCrLf

str = str & "K(1,1) = " & Round(K(1, 1), CLng(Text21.Text)) & vbCrLf= str & "K(1,2) = " & Round(K(1, 2), CLng(Text21.Text)) & vbCrLf= str & "K(1,3) = " & Round(K(1, 3), CLng(Text21.Text)) & vbCrLf= str & "K(1,4) = " & Round(K(1, 4), CLng(Text21.Text)) & vbCrLf= str & "K(2,1) = " & Round(K(2, 1), CLng(Text21.Text)) & vbCrLf= str & "K(2,2) = " & Round(K(2, 2), CLng(Text21.Text)) & vbCrLf= str & "K(2,3) = " & Round(K(2, 3), CLng(Text21.Text)) & vbCrLf= str & "K(2,4) = " & Round(K(2, 4), CLng(Text21.Text)) & vbCrLf= str & "K(3,1) = " & Round(K(3, 1), CLng(Text21.Text)) & vbCrLf= str & "K(3,2) = " & Round(K(3, 2), CLng(Text21.Text)) & vbCrLf= str & "K(3,3) = " & Round(K(3, 3), CLng(Text21.Text)) & vbCrLf= str & "K(3,4) = " & Round(K(3, 4), CLng(Text21.Text)) & vbCrLf= str & "K(4,1) = " & Round(K(4, 1), CLng(Text21.Text)) & vbCrLf= str & "K(4,2) = " & Round(K(4, 2), CLng(Text21.Text)) & vbCrLf= str & "K(4,3) = " & Round(K(4, 3), CLng(Text21.Text)) & vbCrLf= str & "K(4,4) = " & Round(K(4, 4), CLng(Text21.Text)) & vbCrLf= str & "V = " & Round(V, CLng(Text21.Text)) & vbCrLf= str & "Q = " & Round(Q, CLng(Text21.Text)) & vbCrLf= str & "X = " & Round(x, CLng(Text21.Text)) & vbCrLf= str & "Y = " & Round(y, CLng(Text21.Text)).Text = strIfi.Text = Round(V, CLng(Text21.Text)).Text = Round(Q, CLng(Text21.Text)).Text = Round(x, CLng(Text21.Text)).Text = Round(y, CLng(Text21.Text))

'----------------------------------

'Пассивный участок траектории полета

'Эйлер

V = CDbl(Text8.Text)= CDbl(Text9.Text)= CDbl(Text10.Text)= CDbl(Text11.Text)

i = T= Mn - Rashod * T

j = TWhile y > 0= r0 * Exp(-y / 7800)= 0.5 * Cx * r * S * V ^ 2

V = V - h * (9.81 * m * Sin(Q) + Xx) / m

Q = Q - h * 9.81 * Cos(Q) / V= x + h * V * Cos(Q)= y + h * V * Sin(Q)Option1.Value = True And Round(i, 2) = Round(j, 2) Then= j + 1.Text = Text33.Text & Round(i, 2) & vbCrLf.Text = Text34.Text & Round(V, 1) & vbCrLf.Text = Text36.Text & Round(Q, 3) & vbCrLf.Text = Text40.Text & Round(x, 0) & vbCrLf.Text = Text38.Text & Round(y, 0) & vbCrLfIfAbs(Q) < 0.001 Then

str = "Значения параметров" & vbCrLf & "в наивысшей точке траектории" & vbCrLf & "методом Эйлера" & vbCrLf

str = str & "t = " & Round(i, 2) & vbCrLf= str & "V = " & Round(V, CLng(Text21.Text)) & vbCrLf= str & "Q = 0" & vbCrLf= str & "X = " & Round(x, CLng(Text21.Text)) & vbCrLf= str & "Y = " & Round(y, CLng(Text21.Text)) & vbCrLf & vbCrLfIf= i + h.Text = Round(V, CLng(Text21.Text)).Text = Round(Q, CLng(Text21.Text)).Text = Round(x, CLng(Text21.Text)).Text = Round(y, CLng(Text21.Text)).Text = Round(i, 2)

'-----------------------------

'Рунге-Кутт= CDbl(Text12.Text)= CDbl(Text13.Text)= CDbl(Text14.Text)= CDbl(Text15.Text)= T

j = T= Mn - Rashod * T

Do While y > 0= r0 * Exp(-y / 7800)= 0.5 * Cx * r * S * V ^ 2(1, 1) = h * (-9.81 * m * Sin(Q) - Xx) / m(1, 2) = -9.81 * h * Cos(Q) / V(1, 3) = h * V * Cos(Q)(1, 4) = h * V * Sin(Q)= r0 * Exp(-(y + 0.5 * K(1, 4)) / 7800)= 0.5 * Cx * r * S * (V + 0.5 * K(1, 1)) ^ 2(2, 1) = h * (-9.81 * m * Sin(Q + K(1, 2) / 2) - Xx) / m(2, 2) = -9.81 * h * Cos(Q + K(1, 2) / 2) / (V + K(1, 1) / 2)(2, 3) = h * (V + K(1, 1) / 2) * Cos(Q + K(1, 2) / 2)(2, 4) = h * (V + K(1, 1) / 2) * Sin(Q + K(1, 2) / 2)= r0 * Exp(-(y + 0.5 * K(2, 4)) / 7800)= 0.5 * Cx * r * S * (V + 0.5 * K(2, 1)) ^ 2(3, 1) = h * (-9.81 * m * Sin(Q + K(2, 2) / 2) - Xx) / m(3, 2) = -9.81 * h * Cos(Q + K(2, 2) / 2) / (V + K(2, 1) / 2)(3, 3) = h * (V + K(2, 1) / 2) * Cos(Q + K(2, 2) / 2)(3, 4) = h * (V + K(2, 1) / 2) * Sin(Q + K(2, 2) / 2)= r0 * Exp(-(y + K(3, 4)) / 7800)= 0.5 * Cx * r * S * (V + K(3, 1)) ^ 2(4, 1) = h * (-9.81 * m * Sin(Q + K(3, 2)) - Xx) / m(4, 2) = -9.81 * h * Cos(Q + K(3, 2)) / (V + K(3, 1))(4, 3) = h * (V + K(3, 1)) * Cos(Q + K(3, 2))

K(4, 4) = h * (V + K(3, 1)) * Sin(Q + K(3, 2))= V + (K(1, 1) + K(4, 1) + 2 * (K(2, 1) + K(3, 1))) / 6= Q + (K(1, 2) + K(4, 2) + 2 * (K(2, 2) + K(3, 2))) / 6= x + (K(1, 3) + K(4, 3) + 2 * (K(2, 3) + K(3, 3))) / 6

y = y + (K(1, 4) + K(4, 4) + 2 * (K(2, 4) + K(3, 4))) / 6Option2.Value = True And Round(i, 2) = Round(j, 2) Then= j + 1.Text = Text33.Text & Round(i, 2) & vbCrLf.Text = Text34.Text & Round(V, 1) & vbCrLf.Text = Text36.Text & Round(Q, 3) & vbCrLf.Text = Text40.Text & Round(x, 0) & vbCrLf.Text = Text38.Text & Round(y, 0) & vbCrLfIfAbs(Q) < 0.001 Then

str = str & "Значения параметров" & vbCrLf & "в наивысшей точке траектории" & vbCrLf & "методом Рунге-Кутта" & vbCrLf

str = str & "t = " & Round(i, 2) & vbCrLf= str & "V = " & Round(V, CLng(Text21.Text)) & vbCrLf= str & "Q = 0" & vbCrLf= str & "X = " & Round(x, CLng(Text21.Text)) & vbCrLf= str & "Y = " & Round(y, CLng(Text21.Text)).Text = strIf= i + h.Text = Round(V, CLng(Text21.Text)).Text = Round(Q, CLng(Text21.Text)).Text = Round(x, CLng(Text21.Text)).Text = Round(y, CLng(Text21.Text)).Text = Round(i, 2).ForeColor = &HC000&.Caption = "Расчет выполнен"SubSub Command2_Click()

EndSub

Список литературы

. Бахвалов Н.С. Численные методы. т.1., М.: Наука, 1975.

. Волков Е.А. Численные методы. М.: Наука, 1982

. Шалыгин А.С. и др. Автоматизация расчета траекторий ЛА: Учебн. пособие. - Л., 1990.