解决_ivp错误:所需步长小于数字之间的间距

2024-03-29 10:29:28 发布

您现在位置:Python中文网/ 问答频道 /正文

我一直在尝试用Python实现一个不稳定冰川流的模型,用RK45方法在scipy中求解ODE。 可以在here找到原始模型发布

现在,我想我明白了错误是怎么回事,但我找不到修复它的方法。 我不知道它是来自我的实现还是来自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   


Tags: hatplussina0yearsecondsrhoper