Trapezoidal rule for a set of data

[Arman777]

Homework Statement

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


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

Relevant Equations

Tthe question taken from Mark Newman-Computational Physics Exc 5.1

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 ?

 

[scottdave]

It appears to correctly implement the trapezoidal rule to Integrate velocities to estimate position.

I ran the code, and the graph appears correct, but it is difficult to tell what the velocity is doing because of the different scales.
I created another list, called d10, which stores the distance values, divided by 10. These values are just for plotting. Then change the legend.
I'm not sure if it should say position x10 or /10. I think on a tachometer on my car the numbers are 1.0 2.0 etc
and I think it says RPM x1000. Anyway, here is what that graph looks like. With the position scaled down, it is easier to visualize what the velocity is doing.

Also, something you could do is run the formula for the trapezoidal rule in Excel. It only took me a few
minutes to copy your data into a spreadsheet and then create cells to calculate the position.
Then take the numerical output of your python code to see if it matches.

 

[scottdave]

So the max position is 2.5, which represents 25 (meters?)
Probably a better legend for position would be to say (dekameters, or x10 meters), that is if the velocity is meters/sec.

 

[scottdave]

Also you need to take care to note the difference between distance traveled and distance from start (which what I called position)

Think of this as a 25 meter swimming pool, and the swimmer starts from the block, swims down toward the end (positive velocity), then swims back toward the start (negative velocity). Your graph doesn't come back to zero, but when the swimmer makes it all the way back to the start, she would be at position zero again, yet the distance traveled would be 50 meters. What your program calculated was position, rather than distance traveled.

I just came across this post. I am now noticing it is 2.5 months old. Better late than never, I guess.

 

[Arman777]

Thanks for your reply. I ll try to look at it again. Its possible that I have made some mistakes..

 

[scottdave]

I don't know if this function is supposed to return to zero in this context, that was just an example that I came up with of the swimmer returning to the starting point. If you do want to capture "Distance Traveled" then you'd need a different formula. One way would be to integrate in sections. Have a section when the velocity is positive, then have another section when the velocity is negative, then take the absolute value of any "negative distances" and add them together.

The graph of velocity appears to be the shape of a sine wave. You could be able to fit a curve to your data and arrive at an analytical solution for the integral. This would give you an idea of the error associated with estimation.

 

[Arman777]

Have a section when the velocity is positive, then have another section when the velocity is negative, then take the absolute value of any "negative distances" and add them together.

Yes that was on my mind

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

It should be like this.

Thansks for the tips.