龙格-库塔程序输出错误的图形

Question

我正在创建一个程序来比较两个图，使用Runge-Kutta方法来求解ODE。理想情况下，生成的图将直接叠加在一起，但龙格-库塔函数输出的图不正确。我是错误地填充了数组，还是调用新值有问题？

这是我的程序：

#Import correct libraries and extensions
import numpy as np
from scipy.integrate import odeint
from matplotlib import pyplot as plt
import warnings
def fxn():
warnings.warn("deprecated", DeprecationWarning)
with warnings.catch_warnings():
warnings.simplefilter("ignore")
fxn()
#Define conditions and store values
h=0.02 #Given values for ODE
x0=1
y0=5
xpoints=[x0] #array for storing X values
ypoints=[y0] #Array for storing Y values
#State equations and define functions
def dy_dx(y,x):
return x / (np.exp(x) - 1) #Provided Equation
def RungeKuttaFehlberg(x,y):
return x / (np.exp(x) - 1)
#Calculates k1-k4, x and y solutions
def RKFAlg(x0,y0,h):
k1 = RungeKuttaFehlberg(x0,y0)
k2 = RungeKuttaFehlberg(x0+(h/2),y0+((h/2)*k1))
k3 = RungeKuttaFehlberg(x0+(h/2),y0+((h/2)*k2))
k4 = RungeKuttaFehlberg(x0+h,y0+(h*k3))
y1 = y0+(h/6)*(k1+(2*k2)+(2*k3)+k4)
x1 = x0+h
x1 = round(x1,2)
print("Point (x(n+1),y(n+1)) =",(x1,y1))
return((x1,y1)) #Returns as ordered pair
#Define range for number of calculations
for i in range(2000):
print(f"Y{i+1}".format(i)) #Solution value format
x0,y0 = RKFAlg(x0,y0,h) #Calls RKF Function
xpoints.append(x0) #Saves values into array
ypoints.append(y0)
y0 = 1
ODEy1 = odeint(dy_dx,y0,xpoints)
#Runge-Kutta Graph
plt.plot(xpoints,ypoints,'b:',linewidth = 1) #command to plot lines using various colors and widths
plt.suptitle("RKF Graph")
plt.xlabel("x Points")
plt.ylabel("y Points")
plt.show()
#ODE graph
plt.plot(xpoints,ODEy1,'g-',linewidth=1)
plt.suptitle("ODE Graph")
plt.xlabel("x Points")
plt.ylabel("y Points")
plt.show()
#Function for plotting RKF and ODE graph
plt.plot(xpoints,ODEy1,'g-',linewidth=2,label="ODE")
plt.plot(xpoints,ypoints,'b:',linewidth=3,label="Runge-Kutta")
plt.suptitle("ODE and RKF Comparison")
plt.legend(bbox_to_anchor=(.8,1),loc=0,borderaxespad=0)
plt.xlabel("X")
plt.ylabel("Y")
plt.show()
#Function for plotting the difference graph
diff = [] #array to store difference
for i in range(len(xpoints)):
diff.append(ypoints[i]-ODEy1[i])
plt.plot(xpoints,diff)
plt.suptitle("Difference")
plt.xlabel("x Points")
plt.ylabel("RKF and ODE diff.")
plt.show()

Answer 1

在python中，缩进定义了程序的结构、循环中的命令块以及循环后的命令块。

您的目标是通过步骤h在时间间隔内迭代，为x和y构建数组。在这个循环之后，您希望调用odeint一次，以获得x数组上的引用值。

因此将此部分更改为

#Define range for number of calculations
x,y = x0,y0  # for systems use y0.copy()
for i in range(2000):
print(f"Y{i+1}".format(i)) #Solution value format
x,y = RKFAlg(x,y,h) #Calls RKF Function
xpoints.append(x) #Saves values into array
ypoints.append(y)
ODEy1 = odeint(dy_dx,ypoints[0],xpoints)

为了保持一致性，最好保持集成商的接口紧密，

def RK4Alg(func,x0,y0,h):
k1 = func(x0,y0)
...

然后，您可以将相同的rhs函数dy_dx传递给这两个方法，并可以删除错误命名的RungeKuttaFehlberg。

对于合理的参考解决方案，应将误差容限设置为明显低于测试方法预期误差的值。

opts = {"atol":1e-8, "rtol":1e-11}
ODEy1 = odeint(dy_dx,y0,xpoints,**opts)