不打印 if 语句范围内的值(Python 中的 Runge-Kutta 四阶)
Values Within range of if statement not printing (Runge-Kutta Fourth Order in Python)
我正在编写代码来执行具有自适应步长的四阶龙格-库塔数值逼近。
def f(x, t): #integrating function
return -x
def exact(t): #exact solution
return np.exp(-t)
def Rk4(x0, t0, dt): #Runge-Kutta Fourth Order
t = np.arange(0, 1+dt, dt)
n = len(t)
x = np.array([x0]*n)
x[0],t[0] = x0,t0
for i in range(n-1):
h = t[i+1]-t[i]
k1 = h*f(x[i], t[i])
k2 = h*f(x[i] + c21*k1, t[i] + c20*h)
k3 = h*f(x[i] + c31*k1 + c32*k2, t[i] + c30*h)
k4 = h*f(x[i] + c41*k1 + c42*k2 + c43*k3, t[i] + c40*h)
k5 = h*f(x[i] + c51*k1 + c52*k2 + c53*k3 + c54*k4, t[i] + c50*h)
k6 = h*f(x[i] + c61*k1 + c62*k2 + c63*k3 + c64*k4 + c65*k5, t[i] + c60*h)
x[i+1] = x[i] + a1*k1 + a3*k3 + a4*k4 + a5*k5
x5 = x[i] + b1*k1 + b3*k3 + b4*k4 + b5*k5 + b6*k6
E = abs(x[n-1]-exact(t[n-1])) #error
print(E)
if E < 10^-5: #error tolerance
return x[n-1] #approximation
print("For dt = 10e-2, x(1) =",Rk4(1.0,0.0,10e-2))
print("For dt = 10e-3, x(1) =",Rk4(1.0,0.0,10e-3))
然而,当我 运行 代码时,它打印:
5.76914409023e-08
For dt = 10e-2, x(1) = None
4.8151482801e-12
For dt = 10e-3, x(1) = None
因此,在这两种情况下,错误 E 都小于 10^-5,但它不会打印 x(1)。
^
是Python中的bitwise exclusive or,不是求幂;那是 **
。所以 10^-5
(10 xor -5) 是 -15,而你的错误 none 小于 -15。
您可能需要科学记数法常量 1e-5
或者,如果您想将其写为实际的幂运算,10 ** -5
.
对于那些关心性能的人:使用**
的表达式实际上在Python加载脚本时被编译为一个常量,因为参数都是常量,所以它不会有任何性能这样写的影响。 :-)
我正在编写代码来执行具有自适应步长的四阶龙格-库塔数值逼近。
def f(x, t): #integrating function
return -x
def exact(t): #exact solution
return np.exp(-t)
def Rk4(x0, t0, dt): #Runge-Kutta Fourth Order
t = np.arange(0, 1+dt, dt)
n = len(t)
x = np.array([x0]*n)
x[0],t[0] = x0,t0
for i in range(n-1):
h = t[i+1]-t[i]
k1 = h*f(x[i], t[i])
k2 = h*f(x[i] + c21*k1, t[i] + c20*h)
k3 = h*f(x[i] + c31*k1 + c32*k2, t[i] + c30*h)
k4 = h*f(x[i] + c41*k1 + c42*k2 + c43*k3, t[i] + c40*h)
k5 = h*f(x[i] + c51*k1 + c52*k2 + c53*k3 + c54*k4, t[i] + c50*h)
k6 = h*f(x[i] + c61*k1 + c62*k2 + c63*k3 + c64*k4 + c65*k5, t[i] + c60*h)
x[i+1] = x[i] + a1*k1 + a3*k3 + a4*k4 + a5*k5
x5 = x[i] + b1*k1 + b3*k3 + b4*k4 + b5*k5 + b6*k6
E = abs(x[n-1]-exact(t[n-1])) #error
print(E)
if E < 10^-5: #error tolerance
return x[n-1] #approximation
print("For dt = 10e-2, x(1) =",Rk4(1.0,0.0,10e-2))
print("For dt = 10e-3, x(1) =",Rk4(1.0,0.0,10e-3))
然而,当我 运行 代码时,它打印:
5.76914409023e-08
For dt = 10e-2, x(1) = None
4.8151482801e-12
For dt = 10e-3, x(1) = None
因此,在这两种情况下,错误 E 都小于 10^-5,但它不会打印 x(1)。
^
是Python中的bitwise exclusive or,不是求幂;那是 **
。所以 10^-5
(10 xor -5) 是 -15,而你的错误 none 小于 -15。
您可能需要科学记数法常量 1e-5
或者,如果您想将其写为实际的幂运算,10 ** -5
.
对于那些关心性能的人:使用**
的表达式实际上在Python加载脚本时被编译为一个常量,因为参数都是常量,所以它不会有任何性能这样写的影响。 :-)