Numerical Integration MATLAB vs Python/VBA

I have a set of Acceleration-Time data which I am numerically integrating twice to get Velocity-Time and Displacement-Time. When using the cumtrapz function on Matlab it gives me the expected results from the data given. Upon doing the same thing on Python & VBA the results, specifically the displacement-time graph, keeps on rising and rising resulting in non-coinciding values, as shown in the graph below.

What could be the problem in the Python/VBA code for it to be rising up like that, and what does MATLAB do differently which achieves the correct value.

[Torsten]

What could be the problem in the Python/VBA code for it to be rising up like that, and what does MATLAB do differently which achieves the correct value.

How could we know without seeing neither your data nor your code ? And wouldn't this be a question for a Python forum ?

Answer 1

[John D'Errico]

This is difficult to say, because we know virtually nothing.

First, does the problem arise from the use of subtly different data? Far too often, someone sees a difference between what they claim to be the identical mathematics in two different environments, but in fact, they truncated the data to only a few decimal digits when they moved the data. I know, you would NEVER have done that. Or would you?

Next, is this a question of the use of different integration methods. yes, I know people think the integral is the integral. BUT IT ISN'T! ANY numerical integration ends up being just an approximation at some level. Cumtrapz uses a cumulative trapezoidal rule. Did the python tool use a lower order integration (perhaps rectangle rule?) or a higher order integration (perhaps something like a Simpson's rule?)

We don't know any of this. And we are not given that information. The difference LOOKS like a method difference to my eyes, with a gradually growing bias. But if you truncated the accelerations to only a few digits, moving all the accelerations in essentially one direction, that too might have a similar impact.

And of course, we dont even have any data, So there is nothing we can even compare. If you want a better answer, then you need to provide your data. At the very least then we could test it out, changing the method, or truncating the numbers, to see if we could replicate that behavior. Until then, its just an educated guess.

[Gianluca]

Essentially I have time and force data on an excel sheet. There is vast amount (around 10500 cells each of data) which would be quite impossible to attach here. But I double checked my data multiple times and I used the same data on the Matlab file attached. Hope this helps clear it a bit more. If not, thanks for the help anyways.

This is the VBA code I am using to integrate the excel data.

Sub IntegrateForceTime()
    Dim ws As Worksheet
    Set ws = ThisWorkbook.Sheets("VBA")
    
    Dim lastRow As Long
    lastRow = ws.Cells(ws.Rows.Count, "A").End(xlUp).Row
    
    Dim i As Long
    Dim time1 As Double, time2 As Double
    Dim force1 As Double, force2 As Double
    Dim mass As Double
    Dim trapeziumArea As Double
    Dim velocity As Double
    Dim displacement As Double
    
    mass = ws.Cells(2, 3).Value ' Mass of Jumper
    
    ' Set Initial Velocity & Displacement to zero
    velocity = 0
    displacement = 0
    
    ' Set up headers for velocity and displacement columns
    ws.Cells(1, 4).Value = "Velocity (m/s)"
    ws.Cells(1, 5).Value = "Displacement (m)"
    
    For i = 2 To lastRow
        time1 = ws.Cells(i, 1).Value
        force1 = ws.Cells(i, 2).Value
        If i < lastRow Then
            time2 = ws.Cells(i + 1, 1).Value
            force2 = ws.Cells(i + 1, 2).Value
        Else
            time2 = time1
            force2 = force1
        End If
        
        'Integrating acceleration to get velocity using trapezium rule
        trapeziumArea = ((force1 + force2) / 2) * (time2 - time1) / mass
        velocity = velocity + trapeziumArea
        
        'Integrating velocity to get displacement using trapezium rule
        trapeziumArea = ((velocity + (velocity - trapeziumArea)) / 2) * (time2 - time1)
        displacement = displacement + trapeziumArea
        
        ' Output velocity and displacement
        ws.Cells(i, 4).Value = velocity
        ws.Cells(i, 5).Value = displacement
    Next i
End Sub

[Torsten]

Instead of using "cumtrapz", I suggest you try the following code in MATLAB.

If you get the same results as with VBA, you know that the software is not responsible for the different results in MATLAB and VBA, but your integration code.

% Set Initial Velocity & Displacement to zero
    velocity = 0;
    displacement = 0;
    
    mass = data(2,3); % Mass of Jumper
    for i = 2 : lastRow
        time1 = data(i, 1);
        force1 = data(i, 2).Value;
        if i < lastRow 
            time2 = data(i + 1, 1);
            force2 = data(i + 1, 2);
        else
            time2 = time1;
            force2 = force1;
        end
        
        %Integrating acceleration to get velocity using trapezium rule
        trapeziumArea = ((force1 + force2) / 2) * (time2 - time1) / mass;
        velocity = velocity + trapeziumArea;
        
        %Integrating velocity to get displacement using trapezium rule
        trapeziumArea = ((velocity + (velocity - trapeziumArea)) / 2) * (time2 - time1);
        displacement = displacement + trapeziumArea;
        
        % Output velocity and displacement
        data(i, 4) = velocity;
        data(i, 5) = displacement;
    end