我得到以下微分方程 y'' + a y' - y + by^3 = c cos(kx) 初始条件 y(0) = y'(0) = 0 和参数值 a = 0.05, b = k = 1 和 c = 0.5。
现在，我要做的是降低求解器的相对容差，直到数值解 y(t) 与前一次迭代的变化在 x = x_max 时小于 10^−8。
这是我的代码：
from scipy.integrate import odeint
from scipy.optimize import brent
import numpy as np
def de( Y, t ): # de as an array
a = 0.05 # parameters
b = 1.0
c = 0.5
k = 1.0
return np.array([Y[1], Y[0] - a*Y[1] - b*(Y[0])**3 - c*np.cos(k*t)])
def Ydiff(h, *args): # calculates difference between the solution Y1
de = args[0] # (with default relative error) and Y2 (with new
t = args[1] # relative error)
y0 = args[2]
Y1 = args[3]
Y2 = odeint( de, y0, t, full_output=0, rtol=h)
return np.linalg.norm(Y1[-1,0] - Y2[-1,0]) # linalg.norm isn't
# necessary here
t = np.linspace( 0, 100, 10000 )
y0 = np.array([ 0., 0. ]) # initial conditions [y(0),y'(0)]
Y1 = odeint( de, y0, t, full_output=0)
hmin = brent(Ydiff, args=(de,t,y0,Y1), brack=(0,1),full_output=False)
print Ydiff(hmin,de,t,y0,Y1)
我得到了输出：
lsoda-- rtol(i1) is r1 .lt. 0.0
In above message, I1 = 1
In above message, R1 = -0.1618034000000E+01
Illegal input detected (internal error).
Run with full_output = 1 to get quantitative information.
以不同的 R1 值重复...
lsoda-- run aborted.. apparent infinite loop
我不明白为什么我得到“检测到非法输入”和“明显的无限循环”。后者是新的，我不确定为什么会这样。
任何帮助将不胜感激。
PS我很抱歉格式。我是新手，我不知道如何更改它以使问题看起来不错。
-编辑-
我尝试了建议的“full_output = 1”并将 Ydiff 的第 5 行更改为
Y2,p = odeint( de, y0, t, full_output=1,rtol=h)
print p
print
我会得到类似的东西
{'nfe': array([0, 0, 0, ..., 0, 0, 0]), 'nje': array([0, 0, 0, ..., 0, 0, 0]), 'tolsf': array([ 0., 0., 0., ..., 0., 0., 0.]), 'nqu': array([0, 0, 0, ..., 0, 0, 0]), 'lenrw': 52, 'tcur': array([ 0., 0., 0., ..., 0., 0., 0.]), 'hu': array([ -4.10454254e+77, 0.00000000e+00, 0.00000000e+00, ...,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00]), 'imxer': -1, 'leniw': 22, 'tsw': array([ 0., 0., 0., ..., 0., 0., 0.]), 'message': 'Illegal input detected (internal error).', 'nst': array([ 0, -1217644544, 166928456, ..., 2133,
2133, 2133]), 'mused': array([0, 0, 0, ..., 0, 0, 0])}
或者
{'nfe': array([ 7, 9, 11, ..., 4547, 4547, 4547]), 'nje': array([ 0, 0, 0, ..., 12, 12, 12]), 'tolsf': array([ 0., 0., 0., ..., 0., 0., 0.]), 'nqu': array([2, 2, 2, ..., 7, 7, 7]), 'lenrw': 52, 'tcur': array([ 1.37153782e-02, 2.74214740e-02, 4.11275699e-02, ...,
1.00011148e+02, 1.00011148e+02, 1.00011148e+02]), 'hu': array([ 0.0137061 , 0.0137061 , 0.0137061 , ..., 0.07718142,
0.07718142, 0.07718142]), 'imxer': -1, 'leniw': 22, 'tsw': array([ 0. , 0. , 0. , ..., 76.44146089,
76.44146089, 76.44146089]), 'message': 'Integration successful.', 'nst': array([ 3, 4, 5, ..., 2185, 2185, 2185]), 'mused': array([1, 1, 1, ..., 1, 1, 1])}
这很奇怪。今天我没有收到无限循环错误，但我仍然收到这些消息。为了查看发生了什么，我再次将同一行编辑为：
if p['message'] != 'Integration successful.':
print '=======================\nkey \t 2-norm\n'
for k in p.keys():
if type(p[k]) == np.ndarray:
if np.linalg.norm(p[k]) == 0.:
print k,'\t ',np.linalg.norm(p[k],2)
if np.linalg.norm(p[k]) == np.inf:
print k,'\t ',np.linalg.norm(p[k],2)
if np.linalg.norm(p[k]) == np.nan:
print k,'\t ',np.linalg.norm(p[k],2)
print '\n'
我得到了输出（除了'nje'的一次出现为0）：
lsoda-- rtol(i1) is r1 .lt. 0.0
In above message, I1 = 1
In above message, R1 = -0.7482574658091E-07
Illegal input detected (internal error).
Run with full_output = 1 to get quantitative information.
=======================
key 2-norm
tolsf 0.0
tcur 0.0
hu 0.0
tsw 0.0
我想知道如果可能的话如何解决这个问题。