使用来自数组的参数解决ODEs。

huangapple go评论72阅读模式
英文:

Solving ODEs with parameters from an array

问题

I've already solved these 2 differential equations with constant coefficients for ne, Zeff, and Epar using odeint but now I want to solve them passing values at each time step for the coefficients from an array. How can this be done?

def f(Q, t):
    # constants
    ne = 5E19
    ln = 17
    permit = 8.8542E-12
    m_e = 9.1094E-31
    c = 2.9979E8
    e = 1.6E-19
    B0 = 6
    R0 = 0.935
    Zeff = 3
    Epar = 0.2
    Er = (e**3*ne*ln)/(4*np.pi*permit**2*m_e*c**2)
    D = Epar/Er
    alpha = 1+Zeff
    Fgy = (2*permit*B0**2)/(3*ne*m_e*ln)
    Fgc = Fgy*((m_e*c)/(e*B0*R0))**2
    nur = ne*e**4*ln/(4*np.pi*permit**2*m_e**2*c**3)

    # assign each ODE to a vector element
    qp, q = Q
    gamma = (1+q**2)**(1/2)
    qpp2 = q**2-qp**2
    # define each ODE
    dqpdt = (D-gamma*(gamma+alpha)*qp/(q**3)-(Fgc+Fgy*qpp2/q**4)*gamma**4*(((gamma**2)-1)**(1/2)/gamma)**3*qp/q)*nur
    dqdt = (D*qp/q-gamma**2/q**2-(Fgc+Fgy*qpp2/q**4)*gamma**4*(((gamma**2)-1)**(1/2)/gamma)**3)*nur

    return [dqpdt, dqdt]

I've tried using args and odeint but with the simulation time equal to the time step but unsuccessfully.

英文:

I've already solved these 2 differential equations with constant coefficients for ne, Zeff and Epar using odeint but now I want to solve them passing values at each time step for the coefficients from an array. How can this be done?

def f(Q, t):
    # constants
    ne=5E19
    ln=17
    permit=8.8542E-12
    m_e=9.1094E-31
    c=2.9979E8
    e=1.6E-19
    B0=6
    R0=0.935
    Zeff=3
    Epar=0.2
    Er=(e**3*ne*ln)/(4*np.pi*permit**2*m_e*c**2)
    D = Epar/Er
    alpha = 1+Zeff
    Fgy = (2*permit*B0**2)/(3*ne*m_e*ln)
    Fgc = Fgy*((m_e*c)/(e*B0*R0))**2
    nur=ne*e**4*ln/(4*np.pi*permit**2*m_e**2*c**3)
    

    # assign each ODE to a vector element
    qp, q = Q
    gamma=(1+q**2)**(1/2)
    qpp2=q**2-qp**2
    # define each ODE
    dqpdt = (D-gamma*(gamma+alpha)*qp/(q**3)-(Fgc+Fgy*qpp2/q**4)*gamma**4*(((gamma**2)-1)**(1/2)/gamma)**3*qp/q)*nur
    dqdt = (D*qp/q-gamma**2/q**2-(Fgc+Fgy*qpp2/q**4)*gamma**4*(((gamma**2)-1)**(1/2)/gamma)**3)*nur

    return [dqpdt,dqdt]

I've tried using argsand odeint but with the simulation time equal to the time step but unsuccesfully.

答案1

得分: 1

无法完成,因为您无法控制积分中使用的时间步长。您只能控制输出的时间步长,那里的值是从内部积分步长进行插值的(odeint 在概念上是反向进行的,它试图插值下一个输出值,如果不可能,则将内部步长提前,直到可能为止。这样可以在调用底层Fortran例程时获得一个输出值,具有可变数量,可能为零,的内部步长)。

因此,如果您想要与时间有关的参数,您需要将它们作为时间相关的函数传递,以便积分器可以在积分区间的任何时间找到参数。如果参数以函数表的形式给出,您必须提供时间连续函数的最佳猜测(值连续是可选的,零保持或阶跃函数也是可能的),这个过程被称为插值。NumPy 和 SciPy 提供了多种不同插值方法的过程。

英文:

It can't be done because you have no control over the time steps used in the integration. You only control the time steps for the output, the values there get interpolated from the internal integration steps (odeint does this "conceptually" backwards, it tries to interpolate the next output value, and if that is not possible advances the internal step until it is possible. This gives one output value per call of the underlying Fortran routine, with a variable number, possibly zero, of internal steps).

Thus if you want time-dependent parameters, you need to pass them as time-dependent functions, so that the integrator can find the parameters for any time in the integration interval. If the parameters are given as a function table, you have to provide the best guess for a time-continuous function (value-continuous is optional, zero-hold or step functions are possible), this process is called interpolation. Numpy and scipy provide several procedures for different interpolation methods.

huangapple
  • 本文由 发表于 2023年6月12日 19:24:45
  • 转载请务必保留本文链接:https://go.coder-hub.com/76456191.html
匿名

发表评论

匿名网友

:?: :razz: :sad: :evil: :!: :smile: :oops: :grin: :eek: :shock: :???: :cool: :lol: :mad: :twisted: :roll: :wink: :idea: :arrow: :neutral: :cry: :mrgreen:

确定