odeint 函数中的数值问题取决于初始值
Numerical issue in odeint function depending on the initial values
我尝试调试我的代码(求解耦合 ODE)。我发现 odeint 函数有取小数的特定限制 因为这个问题(可能是)我的图对于不同的初始条件是不同的。
问题在第 20 行:x = np.linspace(-30., -50., 5)
问题在第 25 行:state0 = [gi,Li,60.]
当我在 x (-30) 和 state0 (60) 中有初始条件时,代码正常工作:
但是当我更改初始条件 x (-25) 并且在 state0 (50) 中代码不起作用时:
我在 Fortran 中解决了相同的代码,它适用于两个限制。
谁能建议如何解决这个问题。
import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as si
#function for solving differential equations:
def func1(state, x):
g = state[0]
L = state[1]
u = state[2]
phi1 = 1.645
phi2 = 2.* 1.202
N = (-2.*g/(6.*np.pi + 5.*g))*(18./(1. - 2.*L) + 5.*np.log(1.- 2.*L) - phi1 + 6. )
B = (-(2. - N)*L) - ((g/np.pi)* (5.*np.log(1.-2.*L) - phi2 + (5.*N/4.)))
udx = -2.*(1.- ( (3.* 0.5* L)/(3.* 0.330))) #5.1
gdx = (udx) * (2. + N)*g #5.1a
Ldx = (udx) *( B) #5.1b
print gdx, Ldx, udx
return gdx, Ldx, udx
gt = 10.**(-60)
Lt = (1.202*gt)/np.pi
x = np.linspace(-30., -50., 500) #initial conditions
for i in np.arange(len(x)):
gi= np.float64(gt * np.exp(-4.*x[i]))
Li = np.float64(Lt * np.cosh(4.*x[i]))
break
state0 = [gi,Li,60.] #initial conditions
G= si.odeint(func1, state0, x)
a = np.transpose(G) # a[0]=g; a[1]=lambda (L); a[2]=u
plt.plot(a[1],a[0])
plt.title('RGIII Tragectory')
plt.xlabel('L ---->')
plt.ylabel('g ----->')
plt.show()
确实,odeint
求解器正在努力解决这个 ODE 系统。显然,这个 ODE 描述的运动在时间上非常不均匀。解决方法是告诉求解器使用更严格的错误控制:
G = si.odeint(func1, state0, x, rtol=1e-10)
这会为您的两组参数生成一个螺旋图。
一般来说,当odeint
表现怪异时,首先要做的就是请求它的infodict:
G, info = si.odeint(func1, state0, x, full_output=True)
infodict 中有很多信息,但粗略地查看 info['hu']
(使用导致问题的设置)已经显示了所使用的时间步长的荒谬值;就好像求解器决定这个解决方案不会去任何地方,并且在一个时间步长内覆盖了请求的间隔;画一条小线段。
我尝试调试我的代码(求解耦合 ODE)。我发现 odeint 函数有取小数的特定限制 因为这个问题(可能是)我的图对于不同的初始条件是不同的。
问题在第 20 行:x = np.linspace(-30., -50., 5)
问题在第 25 行:state0 = [gi,Li,60.]
当我在 x (-30) 和 state0 (60) 中有初始条件时,代码正常工作:
但是当我更改初始条件 x (-25) 并且在 state0 (50) 中代码不起作用时:
我在 Fortran 中解决了相同的代码,它适用于两个限制。 谁能建议如何解决这个问题。
import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as si
#function for solving differential equations:
def func1(state, x):
g = state[0]
L = state[1]
u = state[2]
phi1 = 1.645
phi2 = 2.* 1.202
N = (-2.*g/(6.*np.pi + 5.*g))*(18./(1. - 2.*L) + 5.*np.log(1.- 2.*L) - phi1 + 6. )
B = (-(2. - N)*L) - ((g/np.pi)* (5.*np.log(1.-2.*L) - phi2 + (5.*N/4.)))
udx = -2.*(1.- ( (3.* 0.5* L)/(3.* 0.330))) #5.1
gdx = (udx) * (2. + N)*g #5.1a
Ldx = (udx) *( B) #5.1b
print gdx, Ldx, udx
return gdx, Ldx, udx
gt = 10.**(-60)
Lt = (1.202*gt)/np.pi
x = np.linspace(-30., -50., 500) #initial conditions
for i in np.arange(len(x)):
gi= np.float64(gt * np.exp(-4.*x[i]))
Li = np.float64(Lt * np.cosh(4.*x[i]))
break
state0 = [gi,Li,60.] #initial conditions
G= si.odeint(func1, state0, x)
a = np.transpose(G) # a[0]=g; a[1]=lambda (L); a[2]=u
plt.plot(a[1],a[0])
plt.title('RGIII Tragectory')
plt.xlabel('L ---->')
plt.ylabel('g ----->')
plt.show()
确实,odeint
求解器正在努力解决这个 ODE 系统。显然,这个 ODE 描述的运动在时间上非常不均匀。解决方法是告诉求解器使用更严格的错误控制:
G = si.odeint(func1, state0, x, rtol=1e-10)
这会为您的两组参数生成一个螺旋图。
一般来说,当odeint
表现怪异时,首先要做的就是请求它的infodict:
G, info = si.odeint(func1, state0, x, full_output=True)
infodict 中有很多信息,但粗略地查看 info['hu']
(使用导致问题的设置)已经显示了所使用的时间步长的荒谬值;就好像求解器决定这个解决方案不会去任何地方,并且在一个时间步长内覆盖了请求的间隔;画一条小线段。