我对编码完全陌生,我想用数值方法解决这5个微分方程。我拿了一个python template并把它应用到我的案例中。以下是我所写内容的简化版本:
import numpy as np
from math import *
from matplotlib import rc, font_manager
import matplotlib.pyplot as plt
from scipy.integrate import odeint
#Constants and parameters
alpha=1/137.
k=1.e-9
T=40.
V= 6.e-6
r = 6.9673e12
u = 1.51856e7
#defining dy/dt's
def f(y, t):
A = y[0]
B = y[1]
C = y[2]
D = y[3]
E = y[4]
# the model equations
f0 = 1.519e21*(-2*k/T*(k - (alpha/pi)*(B+V))*A)
f1 = (3*B**2 + 3*C**2 + 6*B*C + 2*pi**2*B*T + pi**2*T**2)**-1*(-f0*alpha/(3*pi**3) - 2*r*(B**3 + 3*B*C**2 + pi**2*T**2*B) - u*(D**3 - E**3))
f2 = u*(D**3 - E**3)/(3*C**2)
f3 = -u*(D**3 - E**3)/(3*D**2)
f4 = u*(D**3 - E**3)/(3*E**2) + r*(B**3 + 3*B*C**2 + pi**2*T**2*B)/(3*E**2)
return [f0, f1, f2, f3, f4]
# initial conditions
A0 = 2.e13
B0 = 0.
C0 = 50.
D0 = 50.
E0 = C0/2.
y0 = [A0, B0, C0, D0, E0] # initial condition vector
t = np.linspace(1e-15, 1e-10, 1000000) # time grid
# solve the DEs
soln = odeint(f, y0, t, mxstep = 5000)
A = soln[:, 0]
B = soln[:, 1]
C = soln[:, 2]
D = soln[:, 3]
E = soln[:, 4]
y2 = [A[-1], B[-1], C[-1], D[-1], E[-1]]
t2 = np.linspace(1.e-10, 1.e-5, 1000000)
soln2 = odeint(f, y2, t2, mxstep = 5000)
A2 = soln2[:, 0]
B2 = soln2[:, 1]
C2 = soln2[:, 2]
D2 = soln2[:, 3]
E2 = soln2[:, 4]
y3 = [A2[-1], B2[-1], C2[-1], D2[-1], E2[-1]]
t3 = np.linspace(1.e-5, 1e1, 1000000)
soln3 = odeint(f, y3, t3)
A3 = soln3[:, 0]
B3 = soln3[:, 1]
C3 = soln3[:, 2]
D3 = soln3[:, 3]
E3 = soln3[:, 4]
#Plot
rc('text', usetex=True)
plt.subplot(2, 1, 1)
plt.semilogx(t, B, 'k')
plt.semilogx(t2, B2, 'k')
plt.semilogx(t3, B3, 'k')
plt.subplot(2, 1, 2)
plt.loglog(t, A, 'k')
plt.loglog(t2, A2, 'k')
plt.loglog(t3, A3, 'k')
plt.show()
我得到以下错误:
^{pr2}$对于某些参数集,在使用odeint中的mxstep
(也尝试过hmin
和hmax
但没有发现任何区别),虽然错误仍然存在,但我的图看起来很好并且没有受到影响,但大多数时候都是这样。
有时,我得到的错误要求我使用odeint选项full_output=1
运行,这样我得到:
A2 = soln2[:, 0]
TypeError: tuple indices must be integers, not tuple
我不明白这意味着什么。在
我想知道问题在哪里,如何解决。 odeint适合我的工作吗?在
在
lsoda
警告中,t
表示当前时间值,h
表示当前积分步长。由于舍入误差(即r1 + r2 == r1
),步长大小与当前时间之和变得非常接近于当前时间。这类问题通常发生在您正在集成的ode表现不佳时。在在我的机器上,警告消息似乎只在计算
soln2
时出现。在这里,我绘制了警告发生的区域中每个参数的值(注意,为了清晰起见,我切换到了线性轴)。我用红线标记了lsoda
错误第一次出现的时间点(在r1 = 0.2135341098625E-06
)上:警告信息的出现与E中的“扭结”不谋而合,这不是巧合!在
我怀疑这是因为扭结代表了E的梯度中的一个奇点,这迫使解算器采取越来越小的步骤,直到步长达到浮点精度的极限。事实上,你可以在D中看到另一个拐点,它与B中的“摆动”相吻合,可能是由相同的现象引起的。在
关于如何在集成ode时处理奇点的一些一般性建议,请看一下section 5.1.2 here(或任何关于ode的像样的教科书)。在
使用
full_output=True
时出现错误,因为在本例中,odeint
返回一个包含解决方案的元组和一个包含附加输出信息的dict
。试着像这样打开元组:相关问题 更多 >
编程相关推荐