Задача о глубинной бомбе

clear all
r=0.5; s=%pi*r^2; g=9.8;
h=0; ro=1000; mu=1.002; c=0.55;
v=4/3*%pi*r^3;
plo=1000;
FA=plo*g*v;
m=3000*v;
k2=0.5*c*s*ro;
k1=6*%pi*mu*r;
v0=0;
vx(1)=v0;
ax(1)=(m*g-FA-k1*vx(1)-k2*vx(1)^2)/m;
x(1)=0;
y(1)=0;
dt=0.1;
dc=0.1;
i=1;
for t=0: dt:14,
    k2=0.5*c*s*ro;
    vx(i+1)=vx(i)+ax(i)*dt; 
    x(i+1)=x(i)+vx(i)*dt+0.5*ax(i)*dt^2;    
    ax(i+1)=(m*g-FA-k1*vx(i)-k2*vx(i)^2)/m;    
i=i+1;
end
 
i=1;
x2(1)=0; 
vx2(1)=v0;
ax2(1)=(m*g-FA-k1*vx2(1)-k2*vx2(1)^2)/m; 
c=0.4;
for t=0: dt:14,
       k2=0.5*c*s*ro;
    vx2(i+1)=vx2(i)+ax2(i)*dt; 
    x2(i+1)=x2(i)+vx2(i)*dt+0.5*ax2(i)*dt^2;    
    ax2(i+1)=(m*g-FA-k1*vx2(i)-k2*vx2(i)^2)/m;    
i=i+1;
end
 
i=1;
x3(1)=0;
vx3(1)=v0;
ax3(1)=(m*g-FA-k1*vx3(1)-k2*vx3(1)^2)/m;
c=0.045;
for t=0: dt:14,
       k2=0.5*c*s*ro;
    vx3(i+1)=vx3(i)+ax3(i)*dt; 
    x3(i+1)=x3(i)+vx3(i)*dt+0.5*ax3(i)*dt^2;
    ax3(i+1)=(m*g-FA-k1*vx3(i)-k2*vx3(i)^2)/m;
i=i+1;
end
t=0:dt:14+dt,
figure(1); 
subplot(2,2,1);
plot(t,ax,t,ax2,t,ax3); 
xlabel('t');
ylabel('ax');
subplot(2,2,2); 
plot(t,x,t,x2,t,x3);
xlabel('t'); 
ylabel('x');
subplot(2,2,3);
plot(t,vx,t,vx2,t,vx3); 
xlabel('t');
ylabel('vx');

clear all
clc
 
r = 0.5; % радиус сечения, м
S = pi*r^2; % площадь сечения, м
g = 9.8; % ускорение своб падения, м/c^2
p = 1000; % плотность воды, кг/м^3
mu = 1.002; % коэф-т трения 
c = 0.55; % коэф-т формы
V = 0.9; % объем бомбы? м^3
Fa = p*g*V; % Архимедова сила
m = 3000*V; % масса бомбы
Ft = m*g; % сила тяжести
 
% ускорение (производная от скорости)
dv = @(t,v) (Ft - Fa - pi*mu*r*v - 0.5*p*c*S*v^2)/m;
Tx = [0 6]; % пределы значений t
v0 = 0; % начальное знач-е скорости
[t V] = ode45(dv,Tx,v0); % услиенно решаем методом Рунге-Кутта 4-го порядка
for i = 2:length(t)
    h(i) = trapz(t(1:i),V(1:i)); % интеграл скорости (глубина)
end
a = diff(V)./diff(t); % считаем производную скорости (ускорение)
% рисуем
plot(t, V, t, h, t(2:end),a) 
grid on,
legend('Скорость, v(t)','Глубина, h(t)','Ускорение, a(t)','Location','NorthWest')
xlabel('t')

clear all
clc
 
r = 0.5; % радиус сечения, м
S = pi*r^2; % площадь сечения, м
g = 9.8; % ускорение своб падения, м/c^2
p = 1000; % плотность воды, кг/м^3
mu = 1.002; % коэф-т трения
V = 0.9; % объем бомбы? м^3
Fa = p*g*V; % Архимедова сила
m = 3000*V; % масса бомбы
Ft = m*g; % сила тяжести
 
H = 30; % заданная глубина, м
c = [0.55, 0.4, 0.045]; % коэф-ты формы
 
Tmax = 50;
T = linspace(0,Tmax,500);
for i = 1:length(c)
    dv = @(t,v) (Ft - Fa - pi*mu*r*v - 0.5*p*c(i)*S*v^2)/m;  
    [t V] = ode23t(dv,[0 Tmax],0); % услиенно решаем методом Рунге-Кутта 4-го порядка
    h = 0;
    for j = 2:length(t)
        h(j) = trapz(t(1:j),V(1:j)); % считаем глубину
    end
    tau(i) = spline(h,t,H); % время за которое достигнем нужной глубины
    % сохраняем подпись
    str(i) = {['c=',num2str(c(i)),' \tau=',num2str(fix(tau(i)*100)/100),'c.']};
    dept(:,i) = spline(t,h,T); % сохраняем h(t)
end
 
% рисуем
yh = H*ones(1,length(tau));
plot(T, dept, tau, yh,'or'), grid on
ylim([0 H+5])
legend(str,'Location','NorthWest')
xlabel('Время, с')
ylabel('Глубина, м')

clear all
clc
 
r = 0.5; % радиус сечения, м
S = pi*r^2; % площадь сечения, м
g = 9.8; % ускорение своб падения, м/c^2
p = 1000; % плотность воды, кг/м^3
mu = 1.002; % коэф-т трения
V = 0.9; % объем бомбы? м^3
Fa = p*g*V; % Архимедова сила
m = 3000*V; % масса бомбы
Ft = m*g; % сила тяжести
 
H = 30; % заданная глубина, м
c = linspace(0.01,1.5,20); % коэф-ты формы
 
Tmax = 50;
T = linspace(0,Tmax,500);
for i = 1:length(c)
    dv = @(t,v) (Ft - Fa - pi*mu*r*v - 0.5*p*c(i)*S*v^2)/m;  
    [t V] = ode23t(dv,[0 Tmax],0); % услиенно решаем методом Рунге-Кутта 4-го порядка
    h = 0;
    for j = 2:length(t)
        h(j) = trapz(t(1:j),V(1:j)); % считаем глубину
    end
    tau(i) = spline(h,t,H); % время за которое достигнем нужной глубины
end
 
plot(c, tau,'-*r'), grid on
set(gca, 'XTick',0:0.1:1.5)
legend(['Заданная глубина H = ',num2str(H),' м'],'Location','NorthWest')
xlabel('Коэффициент формы, Сх')
ylabel('Время, с')