CODEWITHSUNDEEP
pythonfunctionmatplotlibseabornrunge-kutta
user16852609
user16852609

Reputation: 11

how do i create a graph of multiple x,y points inside a loop?

so I saw this code for RK4 on stack and I found it very useful. However, I cannot figure out a way to plot for each y value at each increment(h) of x.

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))

Upvotes: 1

Views: 357

Answers (2)

Lutz Lehmann
Lutz Lehmann

Reputation: 26040

From a design perspective, it would be preferred if the RK4 code and the plotting code were separated, the numerical solver should not be concerned with how its results are used afterwards.

Then the next decision would be about the construction of the time array, it could be passed to the RK4 method, or be constructed inside and returned, both have advantages. If speed is a concern, the arrays should be constructed explicitly in their final form (see example on math.SE), for expediency one can also construct them incrementally. Thus the code could be changed as

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

and then call as

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()

Upvotes: 0

Ahmad Anis
Ahmad Anis

Reputation: 2714

You can append each point in a list as a tuple, and then perform the line plot operation on the list of tuples. You can find it in the commented code below.

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))

You can see the plot as followsenter image description here

Upvotes: 2

Related Questions