对于某些函数,欧拉能否比龙格-库塔更好?
Can Euler be better than Runge-Kutta for some functions?
我正在尝试解决 Steven Strogatz 的非线性动力学和混沌中的练习。习题2.8.3、2.8.4、2.8.5中,初值问题预计分别实现欧拉法、改进欧拉法、龙格-库塔(四阶)法dx/dt = - X; x(0) = 1 求 x(1)。
经分析,答案是1/e。而且我发现了每种方法中获得的错误。令我惊讶的是,我在 Euler 中得到的错误比在改进的 Euler 和 Runge-Kutta 中要少!
我的代码如下所示。对不起,我的破旧。
from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt
to = 0
xo = 1
tf = 1
deltaT = np.zeros([5])
errorE = np.zeros([5])
errorIE = np.zeros([5])
errorRK = np.zeros([5])
for j in range(0,5):
n = pow(10,j)
deltat = (tf - to)/(n)
print ("delta t is",deltat)
deltaT[j] = deltat
t = np.linspace(to,tf,n)
xE = np.zeros([n])
xIE = np.zeros([n])
xRK = np.zeros([n])
xE[0] = xo
xIE[0] = xo
xRK[0] = xo
for i in range (1,n):
#Regular Euler
xE[i] = deltat*(-xE[i-1]) + xE[i-1]
#Improved Euler
IEintermediate = deltat*(-xIE[i-1]) + xIE[i-1]
xIE[i] = xIE[i-1] - deltat*(xIE[i-1] + IEintermediate)/2
#Runge-Kutta fourth order
k1 = -deltat*xRK[i-1]
k2 = -deltat*(xRK[i-1] + k1/2)
k3 = -deltat*(xRK[i-1] + k2/2)
k4 = -deltat*(xRK[i-1] + k3)
xRK[i] = xRK[i-1] + (k1 + 2*k2 + 2*k3 + k4)/6
print (deltat,xE[i],xIE[i],xRK[i])
errorE[j] = np.exp(-1) - xE[n-1]
errorIE[j] = np.exp(-1) - xIE[n-1]
errorRK[j] = np.exp(-1) - xRK[n-1]
错误:
对于 delT = 1.0
- 欧拉误差为-0.6321205588285577
- I.Euler 错误为 -0.6321205588285577
- RK 误差为 -0.6321205588285577
delT = 0.1
- 欧拉误差-0.019541047828557645
- I.Euler错误-0.039348166379443716
- RK 错误-0.03869055002863331
delT = 0.01
- 欧拉 -0.0018501964782845493
- I.Euler-0.003703427083890265
- RK -0.0036972498815148747
delT = 0.001
- 欧拉 -0.0001840470877806366
- I.Euler-0.00036812480143849635
- RK-0.00036806344222467535
delT = 0.0001
- 欧拉-1.839504510836587e-05
- I.Euler-3.67903967520844e-05
- RK -3.678978357835039e-05
这是合法的吗?如果不是,为什么会这样?
您只进行 n-1
个步长 h=1/n
的积分步骤,因此您计算
exp(-(n-1)/n)=1/e*exp(1/n)
有近似值
1/e + 1/e*1/n
报告的错误值正是 -h/e
,这是一阶的,因此被 1 阶欧拉方法明显扭曲。更确切地说,欧拉值是
(1-1/n)^(n-1) = exp((n-1)*(-1/n-1/(2n^2)+O(1/n^3))
= 1/e*exp(1/(2n)+..)
= 1/e + h/(2e) + ...
如果您调整代码使额外的步骤达到时间 1
,您将得到正确的错误图片。
delta t is 1.0
Euler 0.0 0.367879441171
imp. Euler 0.5 -0.132120558829
Runge-Kutta 4 0.375 -0.00712055882856
delta t is 0.1
Euler 0.3486784401 0.0192010010714
imp. Euler 0.368540984834 -0.00066154366211
Runge-Kutta 4 0.367879774412 -3.33241056083e-07
delta t is 0.01
Euler 0.366032341273 0.00184709989821
imp. Euler 0.367885618716 -6.17754474969e-06
Runge-Kutta 4 0.367879441202 -3.09130498977e-11
delta t is 0.001
Euler 0.367695424771 0.000184016400479
imp. Euler 0.367879502531 -6.13592486265e-08
Runge-Kutta 4 0.367879441171 -4.05231403988e-15
delta t is 0.0001
Euler 0.367861046433 1.83947385133e-05
imp. Euler 0.367879441785 -6.13176398545e-10
Runge-Kutta 4 0.367879441171 -2.6645352591e-15
我正在尝试解决 Steven Strogatz 的非线性动力学和混沌中的练习。习题2.8.3、2.8.4、2.8.5中,初值问题预计分别实现欧拉法、改进欧拉法、龙格-库塔(四阶)法dx/dt = - X; x(0) = 1 求 x(1)。
经分析,答案是1/e。而且我发现了每种方法中获得的错误。令我惊讶的是,我在 Euler 中得到的错误比在改进的 Euler 和 Runge-Kutta 中要少!
我的代码如下所示。对不起,我的破旧。
from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt
to = 0
xo = 1
tf = 1
deltaT = np.zeros([5])
errorE = np.zeros([5])
errorIE = np.zeros([5])
errorRK = np.zeros([5])
for j in range(0,5):
n = pow(10,j)
deltat = (tf - to)/(n)
print ("delta t is",deltat)
deltaT[j] = deltat
t = np.linspace(to,tf,n)
xE = np.zeros([n])
xIE = np.zeros([n])
xRK = np.zeros([n])
xE[0] = xo
xIE[0] = xo
xRK[0] = xo
for i in range (1,n):
#Regular Euler
xE[i] = deltat*(-xE[i-1]) + xE[i-1]
#Improved Euler
IEintermediate = deltat*(-xIE[i-1]) + xIE[i-1]
xIE[i] = xIE[i-1] - deltat*(xIE[i-1] + IEintermediate)/2
#Runge-Kutta fourth order
k1 = -deltat*xRK[i-1]
k2 = -deltat*(xRK[i-1] + k1/2)
k3 = -deltat*(xRK[i-1] + k2/2)
k4 = -deltat*(xRK[i-1] + k3)
xRK[i] = xRK[i-1] + (k1 + 2*k2 + 2*k3 + k4)/6
print (deltat,xE[i],xIE[i],xRK[i])
errorE[j] = np.exp(-1) - xE[n-1]
errorIE[j] = np.exp(-1) - xIE[n-1]
errorRK[j] = np.exp(-1) - xRK[n-1]
错误:
对于 delT = 1.0
- 欧拉误差为-0.6321205588285577
- I.Euler 错误为 -0.6321205588285577
- RK 误差为 -0.6321205588285577
delT = 0.1
- 欧拉误差-0.019541047828557645
- I.Euler错误-0.039348166379443716
- RK 错误-0.03869055002863331
delT = 0.01
- 欧拉 -0.0018501964782845493
- I.Euler-0.003703427083890265
- RK -0.0036972498815148747
delT = 0.001
- 欧拉 -0.0001840470877806366
- I.Euler-0.00036812480143849635
- RK-0.00036806344222467535
delT = 0.0001
- 欧拉-1.839504510836587e-05
- I.Euler-3.67903967520844e-05
- RK -3.678978357835039e-05
这是合法的吗?如果不是,为什么会这样?
您只进行 n-1
个步长 h=1/n
的积分步骤,因此您计算
exp(-(n-1)/n)=1/e*exp(1/n)
有近似值
1/e + 1/e*1/n
报告的错误值正是 -h/e
,这是一阶的,因此被 1 阶欧拉方法明显扭曲。更确切地说,欧拉值是
(1-1/n)^(n-1) = exp((n-1)*(-1/n-1/(2n^2)+O(1/n^3))
= 1/e*exp(1/(2n)+..)
= 1/e + h/(2e) + ...
如果您调整代码使额外的步骤达到时间 1
,您将得到正确的错误图片。
delta t is 1.0
Euler 0.0 0.367879441171
imp. Euler 0.5 -0.132120558829
Runge-Kutta 4 0.375 -0.00712055882856
delta t is 0.1
Euler 0.3486784401 0.0192010010714
imp. Euler 0.368540984834 -0.00066154366211
Runge-Kutta 4 0.367879774412 -3.33241056083e-07
delta t is 0.01
Euler 0.366032341273 0.00184709989821
imp. Euler 0.367885618716 -6.17754474969e-06
Runge-Kutta 4 0.367879441202 -3.09130498977e-11
delta t is 0.001
Euler 0.367695424771 0.000184016400479
imp. Euler 0.367879502531 -6.13592486265e-08
Runge-Kutta 4 0.367879441171 -4.05231403988e-15
delta t is 0.0001
Euler 0.367861046433 1.83947385133e-05
imp. Euler 0.367879441785 -6.13176398545e-10
Runge-Kutta 4 0.367879441171 -2.6645352591e-15