使用微分方程的终端速度

Question

我是 Juia lang 的新手，正在尝试使用 Julia 求解以下微分方程以求出球的最终速度。

F = - m * g - 1/2 rho * v² Cd * A

这是我写的代码：

# Termal velocity of a falling ball

using DifferentialEquations
using Plots

g  = 9.8                # Accelaration of gravity
p  = 1.2                # Density of air
m  = 0.100              # A 100 g ball
r  = 0.10               # 10 cm radius
Cd = 0.5                # Drag coeficient for a small spherical object
y0 = 1000.0             # Initial height of the body (1000 m)
v0 = 10.0               # Initial velocity of the body (10 m/s^2, going up)
A  = pi*r^2;            # Cross-section area of the body;

u0 = [v0;y0]            # Initial Conditions
tspan = (0.0,5.0)       # Time span to solve for
p = [g;p;m;Cd;A]

function Terminal_Velocity(du,u,p,t)
 du[1] = u[1]                                                      # velocity 
 du[2] = -1.0 * p[1] - 0.5 * (p[2]/p[3]) * (u[1]^2) * p[4] * p[5]  # acceleration
end



prob = ODEProblem(Terminal_Velocity,u0,tspan,p)
sol = solve(prob)

plot(sol,vars=(0,1))

我认为问题在于我将 y0 作为加速度的初始条件而不是高度。但是我还不能很好地理解语法。

我的出发点是这篇文章：https://nbviewer.jupyter.org/github/JuliaLang/ODE.jl/blob/master/examples/Terminal_Velocity.ipynb

提前感谢您的帮助。

Answer 1

我认为主要错误来自翻转登录：

du[2] = -1.0 * p[1] - 0.5 * (p[2]/p[3]) * (u[1]^2) * p[4] * p[5]

应该是：

du[2] = +1.0 * p[1] - 0.5 - sign(u[2]) * (p[2]/p[3]) * (u[1]^2) * p[4] * p[5]

但是，p和rho也很容易混淆，因为你在设置ODE参数时重新分配它。

我稍微更改了 ODE 的设置（即 u[1] 现在是位移）。这应该有效：

# Termal velocity of a falling ball
using DifferentialEquations
using Plots

g  = 9.8                # Accelaration of gravity
rho  = 1.2                # Density of air
m  = 0.100              # A 100 g ball
r  = 0.10               # 10 cm radius
Cd = 0.5                # Drag coeficient for a small spherical object
y0 = 1000.0             # Initial height of the body (1000 m)
v0 = 10.0               # Initial velocity of the body (10 m/s^2, going up)
A  = pi*r^2            # Cross-section area of the 
u0 = [y0, v0]            # Initial Conditions
tspan = (0.0,5.0)       # Time span to solve for
p = [g rho m Cd A]

function Terminal_Velocity(du,u,p,t)
 (g, rho, m, Cd, A) = p
 du[1] = u[2]                                                      # velocity 
 du[2] = -g - 0.5 * sign(u[2]) * (rho/m) * (u[2]^2) * Cd * A  # acceleration
end



prob = ODEProblem(Terminal_Velocity,u0,tspan,p)
sol = solve(prob)

p1 = plot(sol, vars=(1), label="Displacement")
p2 = plot(sol, vars=(2), label="Velocity")

plot(p1, p2)

编辑：修正符号错误。

Answer 2

您的示例中有几个错误。他们中的大多数与编程无关，而是与物理和数学有关。

您忽略了阻力项中的符号变化。此外，您在 F 方程中指定的阻力项有一个额外的错误（额外的 1/m）。

您似乎混淆了速度和加速度。 du[2] 是加速度，因为它是速度的导数 (u[2])。您正在使用 u[1] 作为速度。

du[1] = u[1] 给出了 u[1] 的指数增长，你想要的是 du[1] = u[2] 这就是说位置受速度影响。

u0 = [v0;y0]顺序颠倒，u[1]是y坐标，u[2]是速度。

我能看到的唯一编程错误是在选择要绘制的变量时使用基于 0 的索引。

修正以上几点后，你得到：

using DifferentialEquations
using Plots

g  = 9.8                # Accelaration of gravity
p  = 1.2                # Density of air
m  = 0.100              # A 100 g ball
r  = 0.10               # 10 cm radius
Cd = 0.5                # Drag coeficient for a small spherical object
y0 = 1000.0             # Initial height of the body (1000 m)
v0 = 10.0               # Initial velocity of the body (10 m/s^2, going up)
A  = pi*r^2;            # Cross-section area of the body;

u0 = [y0;v0]            # Initial Conditions
tspan = (0.0,5.0)       # Time span to solve for
p = [g;p;m;Cd;A]

function Terminal_Velocity(du,u,p,t)
 du[1] = u[2]                                                      # velocity 
 du[2] = - p[1] - sign(u[2]) * 0.5 * (p[2]/p[3]) * (u[2]^2) * p[4] * p[5]  # acceleration
end

prob = ODEProblem(Terminal_Velocity,u0,tspan,p)
sol = solve(prob)

plt1 = plot(sol; vars=1) 
plt2 = plot(sol; vars=2) 
plot(plt1, plt2)

可以更进一步，使用回调来确保符号更改不会导致数字错误。

为此，请将 solve 行替换为

cond(u, t, i) = u[2]
callback = ContinuousCallback(cond, nothing)
sol = solve(prob; callback=callback)

使用微分方程的终端速度

Terminal Velocity using Differential Equation

julia

scientific-computing

differentialequations.jl