Numerical Integration MATLAB vs Python/VBA

11 views (last 30 days)
Gianluca
Gianluca on 9 Oct 2024
Edited: Torsten on 11 Oct 2024
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.
  1 Comment
Torsten
Torsten on 9 Oct 2024
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 ?

Sign in to comment.

Answers (1)

John D'Errico
John D'Errico on 9 Oct 2024
Edited: John D'Errico on 9 Oct 2024
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.
  2 Comments
Gianluca
Gianluca on 11 Oct 2024
Edited: Gianluca on 11 Oct 2024
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
Torsten on 11 Oct 2024
Edited: Torsten on 11 Oct 2024
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

Sign in to comment.

Products


Release

R2022b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!