如何在循环中创建多个x，y点的图形？

Question

所以我在堆栈上看到了RK4的这段代码，我发现它非常有用。然而，我不能想出一种方法来为x的每个增量(H)的每个y值绘制图。

def f(x,y):
    return 2*x**2-4*x+y

def RK4(x0,y0):
    while x0 < b:
        k1 = h*f(x0,y0)
        k2 = h*f(x0+0.5*h,y0+0.5*k1)
        k3 = h*f(x0+0.5*h,y0+0.5*k2)
        k4 = h*f(x0+h,y0+k3)
        y0+=(k1+2*k2+2*k3+k4)/6
        x0+=h
    return y0

b=3
h=0.001
print(RK4(1,0.7182818))

Answer 1

您可以将列表中的每个点作为元组追加，然后在元组列表上执行线条绘制操作。你可以在下面的注释代码中找到它。

import matplotlib.pyplot as plt


def f(x, y):
    return 2 * x ** 2 - 4 * x + y


def RK4(x0, y0):
    pts = []  # empty list
    while x0 < b:
        k1 = h * f(x0, y0)
        k2 = h * f(x0 + 0.5 * h, y0 + 0.5 * k1)
        k3 = h * f(x0 + 0.5 * h, y0 + 0.5 * k2)
        k4 = h * f(x0 + h, y0 + k3)
        y0 += (k1 + 2 * k2 + 2 * k3 + k4) / 6
        x0 += h
        pts.append((x0, y0))  # appending the tuple
    plt.plot(*zip(*pts))  # plotting the list of tuple
    plt.show()
    return y0


b = 3
h = 0.001
print(RK4(1, 0.7182818))

您可以看到如下图

​

Answer 2

从设计的角度来看，最好是将RK4代码和绘图代码分开，数值求解器不应该关心以后如何使用它的结果。

下一个决定是关于时间数组的构造，它可以传递给RK4方法，也可以在内部构造并返回，这两种方法都有优点。如果速度是一个问题，数组应该以它们的最终形式显式构造(参见example on math.SE)，为了方便起见，人们也可以递增地构造它们。因此，代码可以更改为

def RK4(f,x0,y0,xb,dx):
    x, y = [x0],[y0]
    while x0 < xb:
        k1 = dx*f(x0,y0)
        k2 = dx*f(x0+0.5*dx,y0+0.5*k1)
        k3 = dx*f(x0+0.5*dx,y0+0.5*k2)
        k4 = dx*f(x0+dx,y0+k3)
        y0 += (k1+2*k2+2*k3+k4)/6
        x0 += dx
        x.append(x0); y.append(y0) # for vector y use y0.copy()
    return x,y

然后调用为

x,y = RK4(f=f,x0=1.0,y0=0.7182818,xb=3.0,dx=1e-3)
plt.plot(x,y)
#title, axis labels
plt.grid(); plt.show()