解决二次阻力耦合微分方程。

huangapple go评论88阅读模式
英文:

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(&#39;x (m)&#39;); ylabel(&#39;y (m)&#39;)

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

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

huangapple
  • 本文由 发表于 2020年10月18日 15:21:41
  • 转载请务必保留本文链接:https://go.coder-hub.com/64410747.html
匿名

发表评论

匿名网友

:?: :razz: :sad: :evil: :!: :smile: :oops: :grin: :eek: :shock: :???: :cool: :lol: :mad: :twisted: :roll: :wink: :idea: :arrow: :neutral: :cry: :mrgreen:

确定