1

I want to double integrate the data from a text file. Could someone please help me out? Here is my data. I want the integration of second column:

0.96    2.05
1.52    2.25
2.0     2.36
2.5     2.41
2.96    2.41
3.52    2.43
4.0     2.39
5.04    2.31
6.08    2.25
7.04    2.12
8.08    2.07 
from scipy import integrate
import numpy as np

data = np.loadtxt("A_mod_intensity.txt")
xData, yData = np.hsplit(data,2)

x = xData[:,0]
y = yData[:,0]

integrate.simps(y)
Trenton McKinney
  • 56,955
  • 33
  • 144
  • 158

2 Answers2

0

First of all according to doc you need to provide x values like integrate.simps(y, x), otherwise without x it is assumed that x is 0, 1, 2, 3, 4....

Next if you need to double integrate (which is I think you meaning to integrate two times) you may use code below. In this code we are integrating two times computing indefinite integral numerically. Value of definite integral on the full range will be the last value in array. You may also try code on repl.it.

from scipy import integrate
import numpy as np

data = np.loadtxt("A_mod_intensity.txt")
xData, yData = np.hsplit(data,2)

x = xData[:,0]
y = yData[:,0]

res = y
for i in range(2):
    res = np.array([
        integrate.simps(res[:cnt], x[:cnt]) for cnt in range(1, res.size + 1)
    ])
    print(
        'integrated', i + 1, 'times, definite integral value is', res[-1],
        ', indefinite integral values are', res
    )
Arty
  • 14,883
  • 6
  • 36
  • 69
0

Taking the 2nd integral

  • In order to take the 2nd integral, there must be a function to take the 2nd integral of.
  • In this section of the answer, the 2nd integral of the function must be taken, and then the data is passed to that result.
    • The 2nd integral of the data is not what is occurring here.
  • From your previous question, there are multiple files, that are 0 degrees through 180 degrees, for magnetic field intensity in mT (milliTesla).
  • Assuming a function F(θ) = B * sin(θ), the 2nd integral is -B * sin(θ).
    • This 2nd integral can then be applied to your data
      • With B as intensity and θ extracted from the the df_dict.keys().
      • Convert degrees to radians, with np.deg2rad
  • The shape of the dataframe for this question is different than the shape of the dataframe for the previous question.
import pandas as pd
import numpy as np
from pathlib import path

# create the path to the files in the current working directory
p = Path.cwd()

# get a generator of all the files
files = list(p.glob('*.44V*'))

# load the files into a dict of pandas.DataFrames
dfd = {f'{file.suffix.split("_")[-1]}': pd.read_csv(file, sep='\\s+', header=None) for file in files}

# iterate through the dfd and add a column for for 
for k, v in dfd.items():
    
    # add column headers to the data
    dfd[k].columns = ['mag_field', 'intensity']

    # extract the degree value from the key and convert to radians
    rad = np.deg2rad(int(k[:-3]))
    
    # add a column with the 2nd integral calculation
    dfd[k]['intensity_integral'] = -v.intensity * np.sin(rad)

# display(dfd['10deg'].head())
   mag_field  intensity  integral_intensity
0   269.2575     0.0785            -0.01363
1   269.2594     0.0836            -0.01452
2   269.2613     0.0887            -0.01540
3   269.2632     0.0849            -0.01474
4   269.2651     0.0900            -0.01563

Integrating over columns with simps:

  • At this point, the issue is the shape of the data
  • yData is an (11, 1) array, but needs to be (1, 11)
  • Use .reshape on yData
  • scipy.integrate.simps
    • Integrate y(x) using samples along the given axis and the composite Simpson’s rule.
    • y: array_like. Array to be integrated
    • x: array_like, optional. If given, the points at which y is sampled.
from scipy.integrate import simps
import numpy as np

# reshape yData
y = yData.reshape((1, len(yData)))

# integrate
simps(y)

[out]:
array([23.00666667])

Use x and y

  • This data must use x and y, because the space of the points is not even.
  • The function works across an array, not a single point
y = yData.reshape((1, len(yData)))
x = xData.reshape((1, len(xData)))

simps(y, x)

[out]:
array([16.212197])

Simpson's Rule

Trenton McKinney
  • 56,955
  • 33
  • 144
  • 158