CODEWITHSUNDEEP
pythonenumerate
Yalom
Yalom

Reputation: 113

Why do we use enumerate in this program?

i`m having a class called numerical methods, where we learn how to write programs for certain problems in physics. We had to write 4 programs which could solve ODEs (implicit/explicit euler, velocity-verlet, implicit midpoint rule), now we have to calculate the error by using |y_N - y(T)|. We already have a template which we need to fill out. This is the code which we have to complete. 

def ex2_d():
T = 0.2
y0 = np.array([0.3, 0.0])

all_methods = [explicit_euler, implicit_euler, implicit_mid_point, velocity_verlet]
all_rhs = 3*[pendulum_rhs] + [pendulum_verlet_rhs]
resolutions = 2**np.arange(4, 11)

_, y_exact = ode45(pendulum_rhs, (0.0, T), y0, reltol=1e-12)

for method, rhs in zip(all_methods, all_rhs):
    error = np.empty(resolutions.size)
    for k, N in enumerate(resolutions):
        # TODO: Berechen Sie die Lösung und den Fehler
        error[k] = np.absolute(methode())

    rate = convergence_rate(error, resolutions)
    print(method)
    print("rate: " + str(rate) + "\n")

The only thing I need to fill out is the TODO part. But I don`t understand, the for loop, which is looping over k and N in enumerate(resolution), and why is the resolution array declared as it is anyways?

Thank you in advance for your help!

Upvotes: 0

Views: 107

Answers (2)

Jesse Bakker
Jesse Bakker

Reputation: 2623

In numerically solving an ODE, you want to have doubling resolutions (halving step sizes), to find the convergence rate, using the standard method:

(u_h - u_(h/2))/(u_(h/2) - u_(h/4)) = 2^p + O(h)

with u_h the numerical solution at a step h, u_(h/2) the solution with a step h/2 (e.g. double resolution) and u_(h/4) the solution with a step h/4 (e.g. again double resolution). The order of the error is p, which gives a convergence rate of h^p

This is why the resolutions are declared as 2**np.arange(4,11), which gives[ 16, 32, 64, 128, 256, 512, 1024]`. (You can use other grid sizes, which will change the formula accordingly. For more information, see this.

To store the errors in a list, you need the corresponding indices of the resolutions, which is why enumerate is used:

enumerate(resolutions) -> [(0,16), (1,32), (2,64), (3,128), (4,256), (5,512), (6,1024)]

which is unpacked by the for loop:

iteration    k    N
1            0    16
2            1    32

etc.

Upvotes: 1

ImportanceOfBeingErnest
ImportanceOfBeingErnest

Reputation: 339470

The aim of this excercise is to compare different methods for solving the differential equation given by pendulum_rhs.

The quantity by which the comparison takes place is the convergence rate. In order to determine this rate you need to solve the DE with variing resolution (of the underlying grid) and compute the error for every resolution.

The resolutions to use are given: resolutions =[16, 32, 64, ...].

So for a given method method, you iterate over the resolutions:

for k in range(len(resolutions)):
    N = resolutions[k]
    # calculate the result using N
    result = method(..., N, ...)
    #store it in an array called 
    error[k] = np.abs(y_exact - result)

Upvotes: 0

Related Questions