次の微分方程式 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
可能であれば、これをどのように修正できるのだろうか。