英文:
Solving coupled differential equations of quadratic drag
问题
目标
我一直在尝试解决并绘制属于二次阻力的以下耦合微分方程:
(图片已省略)
方程中的变量定义如下:
c = 0.0004
m = 0.0027
抛射体的初始位置=(0,0.3)
水平分量上的速度=2.05
垂直分量上的速度=0.55
g = 9.81
问题
我似乎无法正确解决这个方程,并且在编程中出现了一些错误。
我的尝试
我已经尝试过在MatLab上使用在线代码,还尝试过在Mathematica上,但它们都无法编写这个方程。我还尝试过查阅Python中的SciPy,但似乎不起作用。
请问有人可以引导我正确编写这个方程的代码吗?
英文:
Goal
I have been attempting to solve and plot the following coupled differential equation belonging to quadratic drag:
The variables from the equations are defined as:
<!-- begin snippet: js hide: false console: true babel: false -->
<!-- language: lang-html -->
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
<!-- end snippet -->
The problem
I cannot seem to solve the equation correctly and have some errors in programming
What I have tried
I have tried using the code from online on MatLab and I have tried on Matematica but none of them is able to program this equation. I have also tried looking at SciPy on Python but it does not seem to work.
Could anybody please set me in the right direction on how to code this properly?
答案1
得分: 1
你可以使用多种MATLAB内置的常微分方程(ODE)求解器。ode45
通常是一个不错的起点。
你有两个位置和两个速度(总共4个状态),因此你需要将4个ODE传递给求解器ode45
(每个状态的导数一个)。
如果x(1)
是x位置,x(2)
是y位置,x(3)
是x速度,x(4)
是y速度,那么x(1)
的导数是x(3)
,x(2)
的导数是x(4)
,x(3)
和x(4)
的导数是由你的两个阻力方程给出的。
最终,MATLAB实现可能如下所示:
c = 0.0004;
m = 0.0027;
p0 = [0; 0.3]; % 初始位置
v0 = [2.05; 0.55]; % 初始速度
g = -9.81; % 重力加速度
tspan = [0 5];
x0 = [p0; v0]; % 初始状态
[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('时间'); ylabel('v_x (m/s)')
subplot(3,1,3); grid on;
plot(t_sol, x_sol(:,4))
xlabel('时间')
ylabel('v_y (m/s)')
英文:
You can use a number of MATLAB built-in ODE solvers. ode45
is usually a good place to start.
You have two positions and two velocities (4 states total), so you need to pass 4 ODEs to the solver ode45
(one derivative for each state).
If x(1)
is the x-position, x(2)
is the y-position, x(3)
is the x-velocity, and x(4)
is the y-velocity, then the derivative of x(1)
is x(3)
, the derivative of x(2)
is x(4)
and the derivatives of x(3)
and x(4)
are the ones given by your two drag equations.
In the end, the MATLAB implementation might look like this:
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
And you can plot the results as follows:
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)')
通过集体智慧和协作来改善编程学习和解决问题的方式。致力于成为全球开发者共同参与的知识库,让每个人都能够通过互相帮助和分享经验来进步。
评论