我是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
谢谢您的帮助。
您的示例中有几个错误。它们中的大多数与编程无关,但与物理和数学相关。
您无视拖动项中的符号更改。另外,您在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)
我认为主要错误来自翻转的登录:
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)
编辑:固定符号错误。