2
\$\begingroup\$

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 ? or Does it looks like a good code in terms of implementation of the trapezoidal rule ?

\$\endgroup\$
0
\$\begingroup\$

A few suggestions:

  • Constant variable names: Any values that do not change in the course of the program should be in UPPERCASE, indicating that they are constants.
  • Unneeded imports: Running your code through pylint, you had two unneeded imports. numpy was not used at all, and ylabel was not used either. Also, this could have just been my machine but the line from pylab ... did not compile, I had to include from matplotlib.pylab import ... in order to run the code.
  • Since this code looks like a one time use, I didn't bother wrapping it in a main guard, but clicking on the link can tell you the benefits of wrapping code that isn't in a function, in a main guard. An excerpt from an answer:

Having a main guard clause in a module allows you to both run code in the module directly and also use procedures and classes in the module from other modules. Without the main guard clause, the code to start your script would get run when the module is imported.

  • You can also convert this into a function, but I will leave that to you :)

Refactored code:

""" Import Statements """
from matplotlib.pylab import loadtxt, plot, show, xlabel, 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

X_INITIAL = TIME[0]  #x_initial
X_FINAL = TIME[-1] #x_final
DIFFERENCE = (X_FINAL - X_INITIAL) // N #difference between each step

S = 0
for k in range(1, N):
    S += FUNCTION[X_INITIAL + k * DIFFERENCE]

#the integral value of the velocity
TOTAL_DISTANCE = DIFFERENCE * (1/2 * (FUNCTION[X_INITIAL] + FUNCTION[X_FINAL]) + S)

#the initial value
DISTANCE = [0]
for k in range(N):
    d = 1/2 * DIFFERENCE * (FUNCTION[X_INITIAL + k * DIFFERENCE] + FUNCTION[X_INITIAL + (k + 1) * DIFFERENCE])
    DISTANCE.append(DISTANCE[-1] + d)

#Create plot graph
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()
\$\endgroup\$
  • \$\begingroup\$ You sure that constant names have to bee uppercaee ? \$\endgroup\$ – Reign Jul 18 at 11:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.