7

私はコーディングがまったく初めてで、これらの 5 つの微分方程式を数値的に解きたいと思っています。Python テンプレートを使用して、自分のケースに適用しました。これが私が書いたものの簡略化されたバージョンです：

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()

次のエラーが表示されます。

lsoda--  warning..internal t (=r1) and h (=r2) are         
   such that in the machine, t + h = t on the next step 
   (h = step size). solver will continue anyway         
  In above,  R1 =  0.2135341098625E-06   R2 =  0.1236845248713E-22

いくつかのパラメーターのセットについて、mxstepin odeint で遊んでみると (これも試しhminてみhmaxましたが違いに気づきませんでした)、エラーは解決しませんが、私のグラフは見栄えがよく、影響を受けていませんが、ほとんどの場合は影響を受けています。odeint オプションfull_output=1を指定して実行するように要求されるエラーが表示されることがあります。

A2 = soln2[:, 0]
TypeError: tuple indices must be integers, not tuple

検索しても意味がわかりません。

どこに問題があり、どのように解決するのかを理解したいと思います。odeint は私がやろうとしていることにも適していますか?

4

1 に答える 1

12

そのlsoda警告でtは、現在の時間値をh参照し、現在の統合ステップ サイズを参照します。ステップ サイズがゼロに近づきすぎたため、丸め誤差 (つまりr1 + r2 == r1) により、現在の時間にステップ サイズを加えたものが現在の時間と等しいと評価されます。この種の問題は通常、統合する ODE の動作が不適切な場合に発生します。

私のマシンでは、警告メッセージは計算時にのみ発生するようsoln2です。ここでは、警告が発生している領域の各パラメーターの値をプロットしました (わかりやすくするために、線形軸に切り替えたことに注意してください)。lsodaエラーが最初に発生した時間ステップ ( r1 = 0.2135341098625E-06) を赤い線でマークしました。

ここに画像の説明を入力

警告メッセージの出現が E の「ねじれ」と一致するのは偶然ではありません。

キンクが E の勾配の特異点を表し、ステップ サイズが浮動小数点精度の限界に達するまで、ソルバーがますます小さなステップを取らざるを得なくなっているのではないかと私は考えています。実際、おそらく同じ現象によって引き起こされた、B の「ぐらつき」と一致する D の別の変曲点を見ることができます。

ODE を統合するときに特異点を処理する方法に関する一般的なアドバイスについては、セクション 5.1.2 を参照してください(または ODE に関する適切な教科書)。

full_output=Trueこの場合、ソリューションを含むタプルと追加の出力情報を含むodeintを返すため、エラーが発生していました。dict次のようにタプルを解凍してみてください。

soln, info = odeint(..., full_output=True)
于 2015-01-23T00:15:04.817 に答える