python-3.x - 如何在 Python 中计算积分估计？

Question

如何在 python3 中对下面列出的提供的数据集执行以下操作？

问题
知道data.txt有 2 列：

xValues，其中 ≤ ≤ ，并且是一些常数
gOfXValues

1计算 '() 的二阶估计
2计算 $g'(x)$ 和 $\int_a^bg(x)dx.$ 的二阶估计
通常，我们不知道从随机数据样本中给出的值是均匀分开。

3绘制 () 和 '() ，打印积分。
4根据图，你认为()是什么函数？
5通过定性比较假设 () 的精确导数和积分与之前获得的数值结果来验证这一点。

到目前为止我做了什么

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()

EDIT_1 鉴于我只能访问 xValue 和 gOfXvalue，我如何才能找到原始函数定义？

Edit_2 基于 D Stanley 的评论

#我们知道在微积分中 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)

Edit_3 基于 Lutz 的回答

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

        

数据集

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 1.790911248707376391e-01
4.136167212263949011e+00 1.614727637626225398e-01
4.167741007777413742e+00 1.446902648395883562e-01
4.199314803290879361e+00 1.287603572615404479e-01
4.230888598804344980e+00 1.136989203067911847e-01
4.262462394317809711e+00 9.952096754324846195e-02
4.294036189831275330e+00 8.624063186256325508e-02
4.325609985344740060e+00 7.387115139215894022e-02
4.357183780858205679e+00 6.242485629917493561e-02
4.388757576371671298e+00 5.191315649949035382e-02
4.420331371885136029e+00 4.234653028407053821e-02
4.451905167398601648e+00 3.373451387397807810e-02
4.483478962912066379e+00 2.608569191446230562e-02
4.515052758425531998e+00 1.940768891759592218e-02
4.546626553938997617e+00 1.370716166199725805e-02
4.578200349452462348e+00 8.989792557207887391e-03
4.609774144965927967e+00 5.260283979342972316e-03
4.641347940479392697e+00 2.522353583661263166e-03
4.672921735992858316e+00 7.787305987243531291e-04
4.704495531506323047e+00 3.115310584367314561e-05
4.736069327019788666e+00 2.803663065220618478e-04
4.767643122533254285e+00 1.526121779621192331e-03
4.799216918046719016e+00 3.767177728993265085e-03
4.830790713560184635e+00 7.001300221330386542e-03
4.862364509073649366e+00 1.122526541299739833e-02
4.893938304587114985e+00 1.643486276363004261e-02
4.925512100100580604e+00 2.262489923329280561e-02
4.957085895614045334e+00 2.978920445901367398e-02
4.988659691127510953e+00 3.792063690553715283e-02
5.020233486640975684e+00 4.701109098416056398e-02
5.051807282154441303e+00 5.705150513256296296e-02
5.083381077667906922e+00 6.803187084756523451e-02
5.114954873181371653e+00 7.994124266182545124e-02
5.146528668694837272e+00 9.276774905451867781e-02
5.178102464208302003e+00 1.064986042851256975e-01
5.209676259721767622e+00 1.211201211385396492e-01
5.241250055235233241e+00 1.366177245687767439e-01
5.272823850748697971e+00 1.529759662277010435e-01
5.304397646262163590e+00 1.701785398642742253e-01
5.335971441775628321e+00 1.882082975789789447e-01
5.367545237289093940e+00 2.070472669172214175e-01
5.399119032802559559e+00 2.266766687846610839e-01
5.430692828316024290e+00 2.470769361666227404e-01
5.462266623829489909e+00 2.682277336329232931e-01
5.493840419342954640e+00 2.901079776086668005e-01
5.525414214856420259e+00 3.126958573908157346e-01
5.556988010369884989e+00 3.359688568895683458e-01
5.588561805883350608e+00 3.599037770728926722e-01
5.620135601396816227e+00 3.844767590918209965e-01
5.651709396910280958e+00 4.096633080634712876e-01
5.683283192423746577e+00 4.354383174880819274e-01
5.714856987937211308e+00 4.617760942757109799e-01
5.746430783450676927e+00 4.886503843576732731e-01
5.778004578964142546e+00 5.160343988571614027e-01
5.809578374477607277e+00 5.439008407929837308e-01
5.841152169991072896e+00 5.722219322897904581e-01
5.872725965504537626e+00 6.009694422676585823e-01
5.904299761018003245e+00 6.301147145834531393e-01
5.935873556531468864e+00 6.596286965958874093e-01
5.967447352044933595e+00 6.894819681258308464e-01
5.999021147558399214e+00 7.196447707829856100e-01
6.030594943071863945e+00 7.500870376296912001e-01
6.062168738585329564e+00 7.807784231523084983e-01
6.093742534098795183e+00 8.116883335102823560e-01
6.125316329612259914e+00 8.427859570327489447e-01
6.156890125125725532e+00 8.740402949322825243e-01
6.188463920639190263e+00 9.054201922051545726e-01
6.220037716152655882e+00 9.368943686873262289e-01
6.251611511666121501e+00 9.684314502351897280e-01
6.283185307179586232e+00 9.999999999999997780e-01

标签： python-3.xmathdifferential-equations

Answer 1

您的导数是错误的，在这个级别应该是

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

但这只是导数的一阶近似，任务要求第二个错误阶。也就是说，对于 处的导数x[i]，您必须通过点x[i-1], x[i], x[i+1]及其值进行插值多项式

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])

并计算它在 处的导数x=x[i]。或者，从泰勒展开式你知道

    (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])+...

将两者结合起来，您可以消除带有 的术语g''(x[i])。

因此，如果

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])

这样写是为了使作为凸组合的性质变得明显。

---

对于积分，累积梯形求积应该足够了。

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

如果要将反导数用作函数（表），请使用累积和。

您可能想直接将数据提取到 numpy 数组中，pandas 中应该有这样做的函数。

---

总的来说，我得到了简短的脚本

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()

导致情节集合

首先表明数据确实是函数g(x)=1+sin(x)，导数正确地看起来像余弦，积分是x+1-cos(x)。