Терминальная скорость с использованием дифференциального уравнения

Question

Я новичок в Джуе Лэнг и пытаюсь решить следующие дифференциальные уравнения, чтобы найти предельную скорость шара, используя Юлию.

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

В вашем примере несколько ошибок. Большинство из них связаны не с программированием, а с физикой и математикой.

Вы игнорируете изменение знака в термине перетаскивания. Кроме того, термин перетаскивания, указанный в уравнении F, имеет дополнительную ошибку (дополнительно 1 / м).

Кажется, вы путаете скорость и ускорение. 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)

Answer 2

Я думаю, что основная ошибка возникла из-за перевернутого знака:

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)

Редактировать: Исправлена ​​ошибка знака.