了解scipy.integrate.odein中的时间步

Question

我试图解决一个偏微分方程使用odeint和线的方法。我的代码肯定是错的-我正在努力找出哪里出错了

我使用odeint(odefunc,y0,tspan)调用ode解算器，其中tspan = np.linspace(0.0, 0.5, 5)&y0 = 1.0*np.ones(3)

我试着在odefunc内打印t，但对输出感到困惑。尽管我正在求解t=0.5，最后一个要打印的t值是0.015081203121127767。输出的数量与tspan匹配，但是我看不出当de函数中的最后一个时间是0.015时，它怎么可能解到t=0.5。我错过了什么

我的数据元素依赖于时间，因此很难找出哪里出了问题，因为我似乎看不到任何事情都失败的时间

ETA-这是失败的，但运行这个没有一些不相关的东西，我得到警告ODEintWarning: Excess work done on this call (perhaps wrong Dfun type). Run with full_output = 1 to get quantitative information.，我假设这是问题的一部分-但它似乎没有停止代码

兆瓦

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import math
import sys

plt.interactive(False) 

sigma = 2320
rho = 1000
gravity = 9.81  # [m/s^2]
g = gravity*3600*3600  # [m/hour^2]
S = 0.01
settlingVelocity = 0.02  # [m/s]
ws = settlingVelocity*3600  # [m/hour]
n = 0.04 # [SI]

J = 400  # [Ws/m]
k = 0.02
Cstar = 0.2 * sigma  # [kg/m^3]
W = 2  # [m]

D0 = 1.2 
Lw = 20   
L = 100  

tend = 0.5  # in hours
tspan = np.linspace(0.0, tend, 5)


def d(t):  # metres
    if t < 50:  # hours
        return 0.5
    else:
        return 0.05

def Q(t):
    return 3600 * (math.sqrt(S)/n)*((W*d(t))**(5/3))/((2*d(t) + W)**(2/3))

def h(t):
    return d(t)/2

def beta(t):
    return (sigma - rho) * g * h(t)/sigma

def Omega(t):
    return rho * g * S * Q(t)  # [W/m]

def PsiTime(t):
    return rho * g * Q(t) * (D0 - d(t))/(Lw)

N = 10
X = np.linspace(0, L, N)
delX = L/ (N-1)

def odefunc(y, t):

    def zetaEh(t):
        return k * (PsiTime(t) + Omega(t)) / (J + beta(t))

    def zetaEW(t): 
        return (2*d(t)/(W + 2*d(t))) * k * Omega(t)/(J + beta(t))

    def zetaR(t):  
        return (W/(W + 2*d(t))) * k*Omega(t)/(beta(t))

    def zetaEF(t,i):
        return (W/(W + 2*d(t))) * k * Omega(t) / (J + beta(t))

    C = y[:N]
    M = y[N:]

    print("time: ", t)

    dCdt = np.zeros(X.shape)
    dMdt = np.zeros(X.shape)

    dCdt[0] = (  # forward difference for dCdx
                -Q(t) / (W*d(t)) * (C[1] - C[0]) / delX
                  + (zetaEh(t) / (W * d(t))) * ((Cstar - C[0]) / Cstar)
                - (ws * C[0] * (beta(t))) / (d(t) * (J + beta(t)))
        )
    dMdt[0] = 0

     # gully channel
    for i in range (1, N-1):  # central difference
        if M[i] + W *C[i] * ws - zetaR(t) * (Cstar - C[i]) / Cstar < 0:
            reMass = M[i] + W * C[i] * ws

            dCdt[i] = (
                    -Q(t) / (W*d(t)) * (C[i+1] - C[i - 1]) / (2*delX)
                        + 1 / (W * d(t)) * ((zetaEW(t) + zetaEF(t,i)) * (Cstar - C[i]) / Cstar
                                            + reMass * (1 - (beta(t))/ (J + beta(t))))
                        - C[i] * ws/d(t)
             )
            dMdt[i] = -M[i]
        else:
            dCdt[i] = (
                    -Q(t) / (W*d(t)) * (C[i+1] - C[i - 1]) / (2*delX)
                      + 1 / (W * d(t)) * (zetaEW(t) + zetaR(t)) * (Cstar - C[i]) / Cstar
                      - C[i] * ws / d(t)
                       )
            dMdt[i] = W * C[i] * ws - zetaR(t) * (Cstar - C[i]) / Cstar

    # Final node - backward difference
    if M[N-1] + W * C[N-1] * ws - zetaR(t) * (Cstar - C[N-1]) / Cstar < 0: 
        reMass = M[N-1] + W * C[N-1] * ws
        dCdt[N-1] = (
                -Q(t) / (W * d(t)) * (C[N-1] - C[N-2]) / delX
                + 1 / (W * d(t)) * ((zetaEW(t) + zetaEF(t, i)) * (Cstar - C[N-1]) / Cstar
                                    + reMass * (1 - (beta(t)) / (J + beta(t))))
                - C[i] * ws / d(t)
        )
        dMdt[N-1] = -M[N-1]

    else:
        dCdt[N-1] = (
                -Q(t) / (W * d(t)) * (C[N-2] - C[N - 1]) / delX
                + 1 / (W * d(t)) * (zetaEW(t) + zetaR(t)) * (Cstar - C[N-1]) / Cstar
                - C[N-1] * ws / d(t)
        )
        dMdt[N-1] = W * C[N-1] * ws - zetaR(t) * (Cstar - C[N-1]) / Cstar
    dydt = np.ravel([dCdt, dMdt])
    return dydt


init_C = 0.0 * np.ones(X.shape)
init_M = 0.0 * np.ones(X.shape)
init= np.ravel([init_C, init_M])

sol = odeint(odefunc, init, tspan)

conc = sol[:, :N]