Trapezoidal rule for set of data

Question

Here is the question from the book of Mark Newman-Computational Physics Exc 5.1

a) Read in the data and, using the trapezoidal rule, calculate from them the approximate distance traveled by the particle in the x direction as a function of time. See Section 2.4.3 on page 57 if you want a reminder of how to read data from a file.

b) Extend your program to make a graph that shows, on the same plot, both the original velocity curve and the distance traveled as a function of time.

the text file of the data.

0   0
1   0.069478
2   0.137694
3   0.204332
4   0.269083
5   0.331656
6   0.391771
7   0.449167
8   0.503598
9   0.554835
10  0.602670
11  0.646912
12  0.687392
13  0.723961
14  0.756491
15  0.784876
16  0.809032
17  0.828897
18  0.844428
19  0.855608
20  0.862439
21  0.864945
22  0.863172
23  0.857184
24  0.847067
25  0.832926
26  0.814882
27  0.793077
28  0.767666
29  0.738824
30  0.706736
31  0.671603
32  0.633638
33  0.593065
34  0.550118
35  0.505039
36  0.458077
37  0.409488
38  0.359533
39  0.308474
40  0.256576
41  0.204107
42  0.151330
43  0.098509
44  0.045905
45  -0.006228
46  -0.057640
47  -0.108088
48  -0.157338
49  -0.205163
50  -0.251347
51  -0.295685
52  -0.337984
53  -0.378064
54  -0.415757
55  -0.450909
56  -0.483382
57  -0.513052
58  -0.539809
59  -0.563563
60  -0.584234
61  -0.601764
62  -0.616107
63  -0.627235
64  -0.635136
65  -0.639814
66  -0.641289
67  -0.639596
68  -0.634786
69  -0.626922
70  -0.616085
71  -0.602366
72  -0.585872
73  -0.566720
74  -0.545039
75  -0.520970
76  -0.494661
77  -0.466272
78  -0.435970
79  -0.403929
80  -0.370330
81  -0.335357
82  -0.299201
83  -0.262054
84  -0.224114
85  -0.185575
86  -0.146636
87  -0.107492
88  -0.068339
89  -0.029370
90  0.009227
91  0.047268
92  0.084574
93  0.120970
94  0.156290
95  0.190375
96  0.223073
97  0.254244
98  0.283753
99  0.311479
100 0.337308

my code

from numpy import array
from pylab import loadtxt, plot, show, xlabel, ylabel, title,legend

values = loadtxt("velocities.txt", float)

time = (values[:, 0]).astype(int)  # time values in float, using astypr to convert them into integer
velocity = values[:, 1]  # velocity values

function = dict(zip(time, velocity))

N = len(time) - 1  # step size

a = time[0]  #x_initial
b = time[-1] #x_final
h = (b - a) // N   #difference between each step

S = 0
for k in range(1, N):
    S += function[a + k * h]

total_distance = h * (1/2 * (function[a] + function[b]) + S)  #the integral value of the velocity

distance = [0]  #the initial value 
for k in range(N):
    d = 1/2 * h * (function[a + k*h] + function[a + (k+1) * h])
    distance.append(distance[-1] + d)

plot(time, distance, "g--", label = "position")
plot(time, velocity, "b-", label= "velocity")
legend(loc='upper left')
xlabel("Time(s)")
title("Velocity vs Time and Distance vs Time")
show()

Is the code correct? Does it looks like good code in terms of implementation of the trapezoidal rule?

Answer 1

Just a few things I noticed.

Utilize built in functions

This

S = 0
for k in range(1, N):
    S += function[a + k * h]

can be this

S = sum(function[a + k * h] for k in range(1, N))

Python3's sum takes an iterable, and returns the sum of all the values in that iterable. So, you can pass in the for loop and it will return the sum for you. Looks neater, too.

Operator Spacing

There should be a space before and after every operator and equal sign. Here's an example from your code:

d = 1/2 * h * (function[a + k*h] + function[a + (k+1) * h])

This should be

d = 1/2 * h * (function[a + k * h] + function[a + (k + 1) * h])

Spacing just helps you and other programmers read your code.

This, however, doesn't apply when passing default parameters. Have a look:

plot(time, distance, "g--", label = "position") # WRONG
plot(time, distance, "g--", label="position") # RIGHT

Variable Naming

A lot of single character variable names can confuse you, and make a program hard to read. I see you have comments such as #x_initial and #x_final. Why not use those as variable names? They're a lot more concise than a and b. h can also be something like diff.

Answer 2

Δ​t

Since you have the timestamps of the velocity measurements, there is no need to calculate the step size. For all you know, the steps are not even all equal. You can use np.diff to calculate the differences between all the data points:

time_differences = (np.diff(time))

or in pure python:

time_differences = [b - a for a, b in zip(time, time[1:])]

or using the pairwise itertools recipe:

time_differences = [b - a for a, b in pairwise(time)]

distance

How much the particle travelled in a certain period is then as simple as velocity[1:] * dt. To calculate the position from this, you take the cumulative sum:

distance= (velocity[1:] * time_differences).cumsum()

or in pure python using itertools.accumulate:

from itertools import accumulate
distance = list(
    accumulate(
        v * time_difference
        for v, time_difference in zip(velocity[1:], time_differences)
    )
)

formatting

I use black (with a maximum line length of 79) to take care of formatting like spacing between operators.

pylab

Pylab takes care of a lot of things for you, but I prefer importing numpy, matplotlib etc. individually. Especially from pylab import ... can become problematic if you have a variable which you want to call title or plot. Which happens a lot for me.