python - solve_ivp 错误:所需的步长小于数字之间的间距
问题描述
我一直在尝试在 Python 中实现不稳定冰川流的模型,使用 RK45 方法解决 scipy 中的 ODE。原始模型出版物可在此处找到。
现在,我想我明白错误发生了什么,但我找不到修复它的方法。我不知道它是来自我的实现还是来自 ODE 本身。我已经多次通过单位,检查所有时间都以秒为单位,所有距离以米为单位等等。我尝试过使用不同的 t_eval 甚至某些常量的不同值,但无法解决我的问题。
我首先创建了一个包含所有常量的类。
import numpy as np
import scipy.integrate
import matplotlib.pyplot as plt
import astropy.units as u
SECONDS_PER_YEAR = 3600*24*365.15
class Cst:
#Glenn's flow Law
A = 2.4e-25
n = 3.
#Standard physical constants
g = 10.#*(u.m)*(u.second**-2)
rho = 916#*(u.kilogram*(u.m**-3))
#Thermodynamics
cp = 2000#**(u.Joule)*(u.kilogram**-1)*(u.Kelvin**-1)
L = 3.3e5#*(u.Joule)*(u.kilogram**-1)
k = 2.1 #*(u.Watt)*(u.m**-1)*'(u.Kelvin**-1)'
DDF = 0.1/SECONDS_PER_YEAR #*(u.m)*(u.yr**-1)*'(u.Kelvin**-1)
K = 2.3e-47#*((3600*24*365.15)**9)#*((u.kilogram**-5)*(u.m**2)*(u.second**9))
C = 9.2e13#*((u.Pascal)*(u.Joule)*(u.m**-2))
#Weertman friction law
q = 1
p = 1/3
R = 15.7#*((u.m**(-1/3))*(u.second**(1/3)))
d = 10#*u.m
sin_theta = 0.05
Tm = 0+273.15 #*u.Kelvin
T_offset = -10+273.15#*u.Kelvin
w = 0.6 #u.m
Wc = 1000.#*u.m
#Velocities
u1 = 0/SECONDS_PER_YEAR #m/s
u2 = 100/SECONDS_PER_YEAR # m/s
#Dimensionless parameters
alpha = 5.
然后我声明了论文中指定的问题特定参数:
#All values are from Table 1
a0 = 1./SECONDS_PER_YEAR#* m/s (u.meter*((u.second)**-1))
l0 = 10000#*(u.meter)
E0 = 1.8e8#(Cst.g*Cst.sin_theta*a0*(l0**2))/(Cst.L*Cst.K))**(1/Cst.alpha)#*(u.Joule/u.m**2)
T0 = 10#E0/(Cst.rho*Cst.cp*Cst.d)#*u.Kelvin
w0 = 0.6#E0/(Cst.rho*Cst.L)#*u.m
N0 = 0.5#Cst.C/E0#*u.Pascal
H0 = 200 #((Cst.R*(Cst.C**Cst.q)*(a0**Cst.p)*(l0**Cst.p))/(Cst.rho*Cst.g*Cst.sin_theta*(E0**Cst.q)))**(1/(Cst.p+1))
t0 = 200 #H0/a0
u0 = 50/SECONDS_PER_YEAR#((Cst.rho*Cst.g*Cst.sin_theta*(E0**Cst.q)*a0*l0)/(Cst.R*(Cst.C**Cst.q)))**(1/(Cst.p+1))
Q0 = (Cst.g*Cst.sin_theta*a0*(l0**2))/Cst.L
S0 = ((Cst.g*Cst.sin_theta*a0*(l0**2)*Cst.Wc)/(Cst.L*Cst.K*((Cst.rho*Cst.g*Cst.sin_theta)**(1/2))))**(3/4)
lamb = ((2.*Cst.A*(Cst.rho*Cst.g*Cst.sin_theta)**Cst.n)*(H0**(Cst.n+1)))/((Cst.n+2)*u0)
chi = N0/(Cst.rho*Cst.g*H0)
gamma = 0.41
kappa = 0.7
phi = 0.2
delta = 66
mu = 0.2
定义模型:
def model(t, x):
#Initial values
H_hat = x[0]
E_hat = x[1]
#Thickness
H = H_hat*H0
#Enthalpy
E_hat_plus = max(E_hat, 0)
E_hat_minus = min(E_hat, 0)
E_plus = E_hat_plus*E0
E_minus = E_hat_minus*E0
a_hat = 1.
theta_hat = Cst.sin_theta/Cst.sin_theta
l_hat =l0/l0
T_a = 0+273.15
T = -10+273.15
# Equation 3
m_hat = (Cst.DDF*(T_a-Cst.T_offset))/a0
S_hat = 0.
T_a_hat = T_a/T0
#Equation A7
if E_plus > 0:
N = min(H/chi, 1./E_plus)
else:
N = H/chi
phi = min(1., E_plus/(H/chi))
#Equation 8
inv_p = 1./Cst.p
u = (Cst.rho*Cst.g*Cst.sin_theta/Cst.R * H * (N**(-Cst.q)))**inv_p
#Equation A7
beta = min(max(0, (u-Cst.u1)/(Cst.u2-Cst.u1)), 1)
#Equation A4
dHdt_hat = (
a_hat - m_hat
+ 1./l_hat*(
theta_hat**inv_p
* H_hat**(1.+inv_p)
* N**(-Cst.q*inv_p)
+ lamb*(theta_hat**Cst.n)
)
)
#Equation A5
dEdt_hat = 1./mu*(
theta_hat**(1+inv_p) * H_hat**(1.+inv_p) * N**(-Cst.q*inv_p)
+ gamma
+ kappa*(E_hat_minus - T_a_hat)/H_hat
- 1./l_hat * (
theta_hat * E_hat_plus**Cst.alpha
+ phi * theta_hat**(1./2) * S_hat**(4/3.)
)
+ delta * beta * m_hat
)
return [dHdt_hat, dEdt_hat]
最后调用它:
tmax = 200*SECONDS_PER_YEAR# *u.years
t = np.linspace(0, tmax, 10000)
sol = scipy.integrate.solve_ivp(model, t_span=[t[0], t[-1]], y0=[1, 1], t_eval=t, method='RK23')
print(sol)
哪个产量
message: 'Required step size is less than spacing between numbers.'
nfev: 539
njev: 0
nlu: 0
sol: None
status: -1
success: False
t: array([0.])
t_events: None
y: array([[1.],
[1.]])
y_events: None
解决方案
推荐阅读
- r - 查找每组列中的最后一个非零元素,填充不同的列
- mongodb - 无法过滤带有空格的项目 MongoDB
- gis - X/Y 坐标到经度/纬度,.net 核心应用程序包
- matlab - 如何在 MATLAB 中为 4D 数据构造分段多项式(三次样条)?
- sql - 为什么这个 Oracle 查询这么慢?
- circleci - 运行测试时 circleci python -t 标志不起作用
- scala - sbt testOnly 带有标签的排除列表不起作用
- r - 具有正向参数的 R stats::step 函数未优化 LR 模型(AIC)
- typescript - 设置语言工作者重试开始计数?
- hugo - 如何获取 Hugo v0.59 中特定部分的前三个帖子?