Проверка корректности

Вначале убедимся в корректности нашего метода (построим 3 траектории: с помощью явного задания уравнений движения, с помощью ODE без сопротивления воздуха и с помощью ODE с учетом сопротивления (но с зануленным баллистическим коэффициентом)). Кроме того необходимо сказать что данные масса, калибр, начальная скорость и баллистический коэффициент взяты для снаряда 76-мм дивизионной пушки образца 1939 года (УСВ):

m-файл:

 

function polet

clear all

global Cx S Pv m a V0 g y

V0 = 680; % начальная скорость

a = 45; % угол бросания (возвышения орудия)

g = 9.8; % ускорение свободного падения

m = 6.3; % масса снаряда

cal = 76; % калибр (в мм)

S = pi * ( cal / 1000 ) ^ 2 % площадь поперечного сечения снаряда (мидель)

Cx = 0; % аэродинамический(баллистический) коэффициент

Pv = 1.225; % плотность воздуха (на уровне моря)

NU = [0; 0]; % координаты в начальный момент времени

t = 100; % время, в течение которого ведется расчет

 

hold

 

% проверка (явное задание)

tp = 0:t;

x = V0 * cosd(a) * tp;

y = V0 * sind(a) * tp - g .* tp.^2 / 2;

plot (x, y, 'r*')

 

%без сопротивления

[T,Y]=ode45(@Simple, [0 t], NU);

plot(Y(:, 1), Y(:, 2), 'b')

 

%с сопротивлением

[T2,Y2]=ode45(@Complex, [0 t], NU);

plot(Y2(:, 1), Y2(:, 2), 'g.')

 

xlabel 'Range'

ylabel 'Height'

axis([0 50000 0 50000])

grid

 

function Simple = Simple(t, x)

global m a V0 g

Simple = [ V0 * cosd(a) ; V0 * sind(a) - g * t ];

 

function Complex = Complex(t, x)

global Cx S Pv m a V0 g

Complex = [ ( 2*m*V0*cosd(a) ) / ( Cx*S*Pv*V0*cosd(a).*t + 2*m ) ;

( 2*m*V0*sind(a) ) / ( Cx*S*Pv*V0*sind(a).*t + 2*m ) - g.*t ];

Графическая интерпретация:

Увеличим:

Как видно из графиков, снаряд данного орудия в отсутствии сопротивления воздуха способен пролететь 45 км. Однако добавление сопротивления среды катастрофически сказывается на дальности:

Получаем вполне реалистичную картину, дальность полета данного снаряда (при угле возвышения 45 градусов) составляет около 5.5 км.

В дальнейшем не будем строить проверочные графики.