Welcome to OStack Knowledge Sharing Community for programmer and developer-Open, Learning and Share
Welcome To Ask or Share your Answers For Others

Categories

0 votes
240 views
in Technique[技术] by (71.8m points)

math - How to calculate Integral estimation in Python?

How can I do the following in python3 on the provided data set listed below?

Problem
Knowing that the data.txt has 2 columns:

xValues, where ?? ≤ ?? ≤ ?? , with ?? and ?? being some constants
gOfXValues

1Compute second order estimates of ??′(??)
2Compute second order estimates of $g'(x)$ and $int_a^b g(x)dx.$
Generally, we don't know that ?? values given from a random data sample is evenly separated.

3Plot ??(??) and ??′(??) , print the integral.
4Based on the graph, what function do you think ??(??) is?
5Verify this by qualitatively comparing the exact derivative and integral of your supposed ??(??) with the numerical results obtained previously.

What I have done so far

import pandas as pd 

dataFrame = pd.read_csv('/Users/Files/data.txt', sep="s+", names=['xValue','gOfXValue'])
dataFrame.info
dataFrame.head()

dgdxArray = []
gOfXValues = []

#iterate all observations and extract the columns as the independent and dependent variable
for index, row in dataFrame.iterrows(): 
    xValue = row['xValue']
    gOfXValue = row['gOfXValue']
    
    gOfXValues.append(gOfXValue)
    
    if index > 0: 
        h = 0.05
       
        difference = gOfXValue - gOfXValues[index-1] #check the difference between Current vs Previous value
        dgdx = difference / h    #get the Derivative
        dgdxArray.append(dgdx)   #add the derivative to an array so as to plot it
        
    

dgdxArray.insert(0,0.5) #hard code values

#plot the initial values provided 
fig = plt.figure(figsize = (6,12))
ax1 = fig.add_subplot(311)
ax1.plot(dataFrame['xValue'], dataFrame['gOfXValue'])
ax1.set_title('Plot initial values x, g(x)')
ax1.set_xlabel('xValue')
ax1.set_ylabel('gOfXValue')


#Plot X value and the derivative on  y axis 
ax2 = fig.add_subplot(312)
ax2.plot(dataFrame['xValue'], dgdxArray)
ax2.set_title('Plot x, derivativeOfGOfx')
ax2.set_xlabel('xValue')
ax2.set_ylabel('gOfXValue')
fig.tight_layout() 
plt.show()

enter image description here

EDIT_1 How can I find the original function definition given that I only have access to xValue and gOfXvalue?

Edit_2 based on comments with D Stanley

#we know that in calculus g(x) = sin(x), g'(x) = cos(x), g(x)dx = - cos(x) + C

#Therefore: 
#calculate sin(x) and compare it to               g(x) provided 
#calculate f'(x)  and compare it to               g(x) provided 
#calculate integral of g(x)dx and compare it to   g(x) provided 

    constant = 2 #choose a random constant

#calculate the sin,cos and integral for existing x value  considering that the original function is sin
sinXfound = np.sin(dataFrame.xValue)
cosXfound = np.cos(dataFrame.xValue)
intXfound = - np.cos(dataFrame.xValue)  + constant

    
#create new columns in the original df with values calculate above
dataFrame['sinXfound'] = sinXfound
dataFrame['cosXfound'] = cosXfound
dataFrame['intXfound'] = intXfound

#find what is the difference between  sin,cos newly found and the original xValue provided in the request
differenceSinXfound = sinXfound - dataFrame['gOfXValue']
differenceCosXfound = cosXfound - dataFrame['gOfXValue']
differenceIntXfound = intXfound - dataFrame['gOfXValue']

#add columns to df
dataFrame['differenceSinXfound'] = differenceSinXfound
dataFrame['differenceCosXfound'] = differenceCosXfound
dataFrame['differenceIntXfound'] = differenceCosXfound

print(dataFrame)

enter image description here

Edit_3 based on Lutz answer

xValues = dataFrame.xValue
gofXValues = dataFrame.gOfXValue


firstDiffArray = []
def calculate_ALL_Divided_Differences():
    for index, row in dataFrame.iterrows():
        if index > 0:
            xNow = row['xValue']
            gNow = row['gOfXValue']
            difference = (gofXValues[indexNow] - gofXValues[indexNow - 1]) / (xValues[indexNow] - xValues[indexNow -1])
            firstDiffArray.append(difference)
            

firstDividedDifference = (gofXValues[1] - gofXValues[0]) / (xValues[1] - xValues[0])

x0 = xValues[0]  
gOfXZero = gofXValues[0]


#Apply Newton's divided difference interpolation formula
for index, row in dataFrame.iterrows(): 
    
    if index > 0:
        xNow = row['xValue']
        gNow = row['gOfXValue']
        x_Minus_x0 = xNow - xValues[0]
        x_Minus_x1 = xNow - xValues[1]
        #Newton's divided difference interpolation formula is
        #f(x) = y0+(x-x0) f [x0,x1]+ (x-x0) * (x-x1) * f [x0,x1,x2]
        
        divided_Difference_Interpolation = gOfXZero + (xNow - x0) * firstDividedDifference + x_Minus_x0 * x_Minus_x1

        

DataSet

0.000000000000000000e+00 1.000000000000000000e+00
3.157379551346525814e-02 1.031568549764810605e+00
6.314759102693051629e-02 1.063105631312673660e+00
9.472138654039577443e-02 1.094579807794844983e+00
1.262951820538610326e-01 1.125959705067717476e+00
1.578689775673262907e-01 1.157214042967250833e+00
1.894427730807915489e-01 1.188311666489717755e+00
2.210165685942568070e-01 1.219221576847691280e+00
2.525903641077220652e-01 1.249912962370308467e+00
2.841641596211873511e-01 1.280355229217014390e+00
3.157379551346525814e-01 1.310518031874168710e+00
3.473117506481178118e-01 1.340371303404112702e+00
3.788855461615830977e-01 1.369885285416546861e+00
4.104593416750483836e-01 1.399030557732340974e+00
4.420331371885136140e-01 1.427778067710209653e+00
4.736069327019788444e-01 1.456099159207016047e+00
5.051807282154441303e-01 1.483965601142838819e+00
5.367545237289094162e-01 1.511349615642326727e+00
5.683283192423747021e-01 1.538223905724288576e+00
5.999021147558398770e-01 1.564561682511917962e+00
6.314759102693051629e-01 1.590336691936528268e+00
6.630497057827704488e-01 1.615523240908179226e+00
6.946235012962356237e-01 1.640096222927107217e+00
7.261972968097009096e-01 1.664031143110431099e+00
7.577710923231661955e-01 1.687304142609184154e+00
7.893448878366314814e-01 1.709892022391332755e+00
8.209186833500967673e-01 1.731772266367076707e+00
8.524924788635619421e-01 1.752923063833377038e+00
8.840662743770272280e-01 1.773323331215339138e+00
9.156400698904925139e-01 1.792952733082778582e+00
9.472138654039576888e-01 1.811791702421020833e+00
9.787876609174229747e-01 1.829821460135725886e+00
1.010361456430888261e+00 1.847024033772298957e+00
1.041935251944353436e+00 1.863382275431223256e+00
1.073509047457818832e+00 1.878879878861460462e+00
1.105082842971284007e+00 1.893501395714874080e+00
1.136656638484749404e+00 1.907232250945481322e+00
1.168230433998214579e+00 1.920058757338174438e+00
1.199804229511679754e+00 1.931968129152434877e+00
1.231378025025145151e+00 1.942948494867437148e+00
1.262951820538610326e+00 1.952988909015839214e+00
1.294525616052075501e+00 1.962079363094462847e+00
1.326099411565540898e+00 1.970210795540986215e+00
1.357673207079006072e+00 1.977375100766707305e+00
1.389247002592471247e+00 1.983565137236369846e+00
1.420820798105936644e+00 1.988774734587002602e+00
1.452394593619401819e+00 1.992998699778669724e+00
1.483968389132867216e+00 1.996232822271006846e+00
1.515542184646332391e+00 1.998473878220378808e+00
1.547115980159797566e+00 1.999719633693477938e+00
1.578689775673262963e+00 1.999968846894156327e+00
1.610263571186728138e+00 1.999221269401275647e+00
1.641837366700193535e+00 1.997477646416338626e+00
1.673411162213658709e+00 1.994739716020657028e+00
1.704984957727123884e+00 1.991010207442792002e+00
1.736558753240589281e+00 1.986292838338002742e+00
1.768132548754054456e+00 1.980592311082403967e+00
1.799706344267519631e+00 1.973914308085537694e+00
1.831280139780985028e+00 1.966265486126021811e+00
1.862853935294450203e+00 1.957653469715929573e+00
1.894427730807915378e+00 1.948086843500509424e+00
1.926001526321380775e+00 1.937575143700825064e+00
1.957575321834845949e+00 1.926128848607841171e+00
1.989149117348311346e+00 1.913759368137436745e+00
2.020722912861776521e+00 1.900479032456751760e+00
2.052296708375241696e+00 1.886301079693208704e+00
2.083870503888706871e+00 1.871239642738459663e+00
2.115444299402172490e+00 1.855309735160411755e+00
2.147018094915637665e+00 1.838527236237377238e+00
2.178591890429102840e+00 1.820908875129262583e+00
2.210165685942568015e+00 1.802472214201578105e+00
2.241739481456033189e+00 1.783235631518890418e+00
2.273313276969498808e+00 1.763218302525168424e+00
2.304887072482963983e+00 1.742440180929283100e+00
2.336460867996429158e+00 1.720921978814716535e+00
2.368034663509894333e+00 1.698685145993306556e+00
2.399608459023359508e+00 1.675751848623608709e+00
2.431182254536824683e+00 1.652144947115186779e+00
2.462756050050290302e+00 1.627887973340858441e+00
2.494329845563755477e+00 1.603005107179614530e+00
2.525903641077220652e+00 1.577521152413588590e+00
2.557477436590685826e+00 1.551461512003107668e+00
2.589051232104151001e+00 1.524852162764468444e+00
2.620625027617616620e+00 1.497719629475682934e+00
2.652198823131081795e+00 1.470090958436002904e+00
2.683772618644546970e+00 1.441993690505579018e+00
2.715346414158012145e+00 1.413455833652134119e+00
2.746920209671477320e+00 1.384505835032010967e+00
2.778494005184942495e+00 1.355172552633428618e+00
2.810067800698408114e+00 1.325485226510211501e+00
2.841641596211873289e+00 1.295473449634670038e+00
2.873215391725338463e+00 1.265167138398678670e+00
2.904789187238803638e+00 1.234596502792368877e+00
2.936362982752268813e+00 1.203792016290152311e+00
2.967936778265734432e+00 1.172784385474099356e+00
2.999510573779199607e+00 1.141604519424951558e+00
3.031084369292664782e+00 1.110283498911275091e+00
3.062658164806129957e+00 1.078852545407476660e+00
3.094231960319595132e+00 1.047342989971558280e+00
3.125805755833060751e+00 1.015786242013636764e+00
3.157379551346525925e+00 9.842137579863630137e-01
3.188953346859991100e+00 9.526570100284420528e-01
3.220527142373456275e+00 9.211474545925236734e-01
3.252100937886921450e+00 8.897165010887251313e-01
3.283674733400387069e+00 8.583954805750482198e-01
3.315248528913852244e+00 8.272156145259004223e-01
3.346822324427317419e+00 7.962079837098479107e-01
3.378396119940782594e+00 7.654034972076313448e-01
3.409969915454247769e+00 7.348328616013215520e-01
3.441543710967712943e+00 7.045265503653301842e-01
3.473117506481178562e+00 6.745147734897881664e-01
3.504691301994643737e+00 6.448274473665716044e-01
3.536265097508108912e+00 6.154941649679892546e-01
3.567838893021574087e+00 5.865441663478661027e-01
3.599412688535039262e+00 5.580063094944210933e-01
3.630986484048504881e+00 5.299090415639968743e-01
3.662560279561970056e+00 5.022803705243168437e-01
3.694134075075435231e+00 4.751478372355316671e-01
3.725707870588900406e+00 4.485384879968926652e-01
3.757281666102365580e+00 4.224788475864115211e-01
3.788855461615830755e+00 3.969948928203856919e-01
3.820429257129296374e+00 3.721120266591415593e-01
3.852003052642761549e+00 3.478550528848133316e-01
3.883576848156226724e+00 3.242481513763912915e-01
3.915150643669691899e+00 3.013148540066936665e-01
3.946724439183157074e+00 2.790780211852837978e-01
3.978298234696622693e+00 2.575598190707169000e-01
4.009872030210087424e+00 2.367816974748317982e-01
4.041445825723553043e+00 2.167643684811096927e-01
4.073019621237018661e+00 1.975277857984218954e-01
4.104593416750483392e+00 

与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…
Welcome To Ask or Share your Answers For Others

1 Answer

0 votes
by (71.8m points)

Your derivative is wrong, at this level it should be

(gofx[i]-gofx[i-1]) / (x[i]-x[i-1])

But this is only a first order approximation of the derivative, the task asks for a second error order. That is, for the derivative at x[i], you have to take the interpolation polynomial through the points x[i-1], x[i], x[i+1] and their values

g[x[i]] + g[x[i],x[i+1]] * (x-x[i]) + g[x[i],x[i-1],x[i+1]] * (x-x[i])*(x-x[i-1])

and compute the derivative of it at x=x[i]. Or alternatively, from the Taylor expansion you know that

    (gofx[i]-gofx[i-1]) / (x[i]-x[i-1]) = g'(x[i]) - 0.5*g''(x[i])*(x[i]-x[i-1])+...
    (gofx[i+1]-gofx[i]) / (x[i+1]-x[i]) = g'(x[i]) + 0.5*g''(x[i])*(x[i+1]-x[i])+...

Combining both you can eliminate the term with g''(x[i]).

so if

dx = x[1:]-x[:-1]
dg = g[1:]-g[:-1]

are the simple differences, then the first order derivative with second error order is

dg_dx = dg/dx
diff_g = ( dx[:-1]*(dg_dx[1:]) + dx[1:]*(dg_dx[:-1]) ) / (dx[1:]+dx[:-1])

This is written so that the nature as convex combination becomes obvious.


For the integral, the cumulative trapezoidal quadrature should be sufficient.

sum( 0.5*(g[:-1]+g[1:])*(x[1:]-x[:-1]) )

Use the cumulative sum if you want the anti-derivative as function (table).

You might want to extract the data into numpy arrays directly, there should be functions in pandas that do that.


In total I get the short script

x,g = np.loadtxt('so65602720.data').T

%matplotlib inline
plt.figure(figsize=(10,3))
plt.subplot(131)
plt.plot(x,g,x,np.sin(x)+1); plt.legend(["table g values", "$1+sin(x)$"]); plt.grid(); 

dx = x[1:]-x[:-1]
dg = g[1:]-g[:-1]
dg_dx = dg/dx
diff_g = ( dx[:-1]*(dg_dx[1:]) + dx[1:]*(dg_dx[:-1]) ) / (dx[1:]+dx[:-1])

plt.subplot(132)
plt.plot(x,g,x[1:-1],diff_g); plt.legend(["g", "g'"]); plt.grid(); 

int_g = np.cumsum(0.5*(g[1:]+g[:-1])*(x[1:]-x[:-1]))

plt.subplot(133)
plt.plot(x[1:],int_g,x,x); plt.legend(["integral of g","diagonal"]); plt.grid(); 

plt.tight_layout(); plt.show()

resulting in the plot collection

plots of function, derivative and integral

showing first that indeed the data is of the function g(x)=1+sin(x), that the derivative correctly looks like the cosine and the integral is x+1-cos(x).


与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…
Welcome to OStack Knowledge Sharing Community for programmer and developer-Open, Learning and Share
Click Here to Ask a Question

...