java - 求解二次阻力的耦合微分方程

Question

目标

我一直在尝试求解并绘制以下属于二次阻力的耦合微分方程：

方程中的变量定义为：

c = 0.0004
m = 0.0027
starting position of projectile= (0,0.3)
velocity on horizontal component = 2.05
velocity on vertical component = 0.55
g = 9.81

问题

我似乎无法正确求解方程并且在编程中出现了一些错误

我试过的

我曾尝试在 MatLab 上使用在线代码，并在 Matematica 上尝试过，但没有一个能够编写这个方程。我也尝试过在 Python 上查看 SciPy，但它似乎不起作用。

任何人都可以为我正确编码如何正确编码吗？

标签： javapythonmatlabdragprojectile

Answer 1

您可以使用许多 MATLAB 内置 ODE 求解器。ode45通常是一个很好的起点。

您有两个位置和两个速度（总共 4 个状态），因此您需要将 4 个 ODE 传递给求解器ode45（每个状态一个导数）。如果x(1)是 x 位置，x(2)是 y 位置，x(3)是 x 速度，x(4)是 y 速度，那么x(1)isx(3)的导数， is 的导数和x(2)的x(4)导数是两个阻力方程给出的导数.x(3)x(4)

最后，MATLAB 实现可能如下所示：

c = 0.0004;
m = 0.0027;
p0 = [0; 0.3]; % starting positions 
v0 = [2.05; 0.55]; % starting velocities
g = -9.81; % gravitational acceleration

tspan = [0 5];
x0 = [p0; v0]; % initial states
[t_sol, x_sol] = ode45(@(t,x) drag_ode_fun(t,x,c,m,g), tspan, x0);

function dxdt = drag_ode_fun(t,x,c,m,g)
   dxdt = zeros(4,1);
   dxdt(1) = x(3);
   dxdt(2) = x(4);
   dxdt(3) = -c/m*x(3)*sqrt(x(3)^2+x(4)^2);
   dxdt(4) = g-c/m*x(4)*sqrt(x(3)^2+x(4)^2);
end

您可以按如下方式绘制结果：

figure; 
subplot(3,1,1); grid on;
plot(x_sol(:,1), x_sol(:,2))
xlabel('x (m)'); ylabel('y (m)')

subplot(3,1,2); grid on;
plot(t_sol, x_sol(:,3))
xlabel('time'); ylabel('v_x (m/s)')

subplot(3,1,3); grid on;
plot(t_sol, x_sol(:,4))
xlabel('time')
ylabel('v_y (m/s)')