找回密码
 欢迎注册
查看: 17930|回复: 12

[原创] 在阶乘中间插值,能否插出零点五根号派?

[复制链接]
发表于 2019-1-19 19:35:56 | 显示全部楼层 |阅读模式
一、

假如有人这样问我:已知$f$是一个函数,$f(0)=1$,$f(1)=1$,问$f(0.5)$等于多少?

那么我就会毫不犹豫地回答:$f(0.5)=1$。

假如那个人又问我:为什么这样回答?

那么我就会这样回答:因为要使得数列【$1$、$x$、$1$】的二阶差分为$0$,那么$x$只能等于$1$呀~

二、

假如有人这样问我:已知$f$是一个函数,$f(0)=1$,$f(1)=1$,$f(2)=2$,问$f(0.5)$和$f(1.5)$分别等于多少?

那么我就会毫不犹豫地回答:$f(0.5)=0.875$,$f(1.5)=1.375$。

假如那个人又问我:为什么这样回答?

那么我就会这样回答:因为要使得数列【$1$、$x$、$1$、$y$、$2$】的三阶差分恒为$0$,那么$x$和$y$只能等于$0.875$和$1.375$呀~

不信你看:

零阶差分: $1.000$,$0.875$,$1.000$,$1.375$,$2.000$
一阶差分: $-0.125$,$0.125$,$0.375$,$0.625$
二阶差分:   $0.250$,$0.250$,$0.250$
三阶差分:      $0.000$,$0.000$

看到了吧?三阶差分恒为$0$哦~这才是一个光滑的函数嘛~

三、

假如有人这样问我:已知$f$是一个函数,$f(0)=1$,$f(1)=1$,$f(2)=2$,$f(3)=6$,问$f(0.5)$、$f(1.5)$和$f(2.5)$分别等于多少?

那么我就会毫不犹豫地回答:$f(0.5)=1$,$f(1.5)=1.25$,$f(2.5)=3.5$。

假如那个人又问我:为什么这样回答?

那么我就会这样回答:因为要使得数列【$1$、$x$、$1$、$y$、$2$、$z$、$6$】的四阶差分恒为$0$,那么$x$、$y$和$z$只能等于$1$、$1.25$和$3.5$呀~

不信你看:
  1.   1  1  1   1.25    2   3.5   6
  2.    0  0  0.25  0.75  1.5  2.5
  3.     0  0.25  0.5  0.75  1
  4.      0.25 0.25 0.25 0.25
  5.          0    0    0
复制代码

看到了吧?四阶差分恒为$0$哦~这才是一个光滑的函数嘛~

四、

假如有人这样问我:已知$f$是一个函数,$f(0)=1$,$f(1)=1$,$f(2)=2$,$f(3)=6$,$f(4)=24$,问$f(0.5)$、$f(1.5)$、$f(2.5)$和$f(3.5)$分别等于多少?

那么我就会毫不犹豫地回答:$f(0.5)=0.6484375$,$f(1.5)=1.4609375$,$f(2.5)=3.1484375$,$f(3.5)=12.2109375$。

假如那个人又问我:为什么这样回答?

那么我就会这样回答:因为要使得数列【$1$、$w$、$1$、$x$、$2$、$y$、$6$、$z$、$24$】的五阶差分恒为$0$,那么$w$、$x$、$y$和$z$只能等于$0.6484375$、$1.4609375$、$3.1484375$和$12.2109375$呀~这才是一个光滑的函数嘛~

五、

……$f(5)=120$……

……$1.8515625,\ 0.9453125,\ 3.6640625,\ 11.0078125,\ 54.9765625$……

……

依次类推,最终$f(0.5)$的极限会不会等于$\sqrt{\pi}/2$呢?
单选投票, 共有 4 人参与投票

距结束还有: 8053 天9 小时51 分钟

您所在的用户组没有投票权限
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2019-1-19 20:26:25 | 显示全部楼层
本帖最后由 .·.·. 于 2019-1-19 20:33 编辑

我用脏脏的R代码给出了一个解答
  1. > optim(sapply(1:3+0.5,gamma),function(x)sum(diff(diff(diff(diff(c(1,x[1],1,x[2],2,x[3],6)))))^2),method='BFGS',control=list(reltol=1e-60))
  2. $par
  3. [1] 1.00 1.25 3.50
复制代码
四阶:
  1. >  optim(sapply(1:4+0.5,gamma),function(x)sum(diff(diff(diff(diff(diff(c(1,x[1],1,x[2],2,x[3],6,x[4],24))))))^2),method='BFGS',control=list(reltol=1e-16))
  2. $par
  3. [1]  0.6484375  1.4609375  3.1484375 12.2109375
复制代码
五阶:
  1. >  optim(sapply(1:5+0.5,gamma),function(x)sum(diff(diff(diff(diff(diff(diff(c(1,x[1],1,x[2],2,x[3],6,x[4],24,x[5],120)))))))^2),method='BFGS',control=list(reltol=1e-16))
  2. $par
  3. [1]  1.8515625  0.9453125  3.6640625 11.0078125 54.9765625
复制代码
六阶
  1. >  optim(sapply(1:6+0.5,gamma),function(x)sum(diff(diff(diff(diff(diff(diff(diff(c(1,x[1],1,x[2],2,x[3],6,x[4],24,x[5],120,x[6],720))))))))^2),method='BFGS',control=list(reltol=1e-16))
  2. $par
  3. [1]  -3.583008   2.756836   2.370117  12.819336  49.541992 302.350586
复制代码
七阶
  1. > optim(sapply(1:7+0.5,gamma),function(x)sum(diff(diff(diff(diff(diff(diff(diff(diff(c(gamma(1),x[1],gamma(2),x[2],gamma(3),x[3],gamma(4),x[4],gamma(5),x[5],gamma(6),x[6],gamma(7),x[7],gamma(8))))))))))^2),method='BFGS',control=list(reltol=1e-16))
  2. $par
  3. [1]   26.291016   -5.390625    6.896484    8.292969   57.689453  272.476562 1965.294922
复制代码
不准备继续插值了
毕竟荣格效应摆在那里的……

点评

你这样解释,我就找对参考文献了,就知道【不能插出零点五根号派】的原因所在了。谢谢你~  发表于 2019-1-20 23:16
https://baike.baidu.com/item/%E9%BE%99%E6%A0%BC%E7%8E%B0%E8%B1%A1/5473475?fr=aladdin理论上这才是这个效应的正确命名……毕竟这个名字是当年帮室友写计算原理的时候偷学来的……记混了也是情有可原的  发表于 2019-1-20 15:25
就是多项式差值查到最后总会出问题。 一般差值,为了稳健都是指明一点值与各阶导数 ——如果指明n个点,进行拟合的时候多半会得到一个奇奇怪怪的波动剧烈的函数  发表于 2019-1-20 15:24
虽然你是对的,但是我不懂荣格效应,你能解释一下什么是荣格效应吗?  发表于 2019-1-20 06:19
看样子你对。看来伽马函数不是这么插值插出来的。  发表于 2019-1-19 20:40

评分

参与人数 1威望 +3 金币 +4 贡献 +3 经验 +3 鲜花 +3 收起 理由
KeyTo9_Fans + 3 + 4 + 3 + 3 + 3 点评里的链接有很大的参考价值

查看全部评分

毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
回复

使用道具 举报

发表于 2019-1-19 21:42:23 | 显示全部楼层
这个显然不行吧,你这种插值本质上也只是用多项式逼近而已,根据斯特灵公式,$n !$显然不能展开为收敛的多项式级数。

直觉上应该是要去对$\ln n!$进行插值而不是$n!$,但根据斯特灵公式,其实对$\ln n!$进行插值也不可能完全准确,但至少比直接对$n!$插值要好。

点评

对超越函数用多项式差值通常都是极其不稳定的  发表于 2019-1-19 21:48
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
回复

使用道具 举报

发表于 2019-1-19 21:45:27 | 显示全部楼层
$\Gamma$函数作为解析函数,其接近分式函数的表达式的极限形式为
\(\Gamma(z)=\frac{e^{-\gamma z}}z\prod_{n=1}^{\infty}(1+\frac{z}n)^{-1}e^{\frac{z}n}\)
其实其倒数的形式比较好,是个整函数
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
回复

使用道具 举报

发表于 2019-2-1 19:28:07 | 显示全部楼层
楼主真坏,刚好在六阶的时候 数据就飘了。却不说出来。

以下的代码给出 插值出来的多项式表达式,以及$f(1/2)$的精确值

  1. Table[{n, f = Expand[InterpolatingPolynomial[Range[0, n]!, x]],  t = f /. x -> 3/2, N[t, 10]}, {n, 1, 10}] // Column
复制代码


(通过双击鼠标可以点选特定的行)
{1,1,1,1.000000000}
{2,2-(3 x)/2+x^2/2,7/8,0.8750000000}
{3,(13 x)/6-(3 x^2)/2+x^3/3,1,1.000000000}
{4,9-(199 x)/12+(93 x^2)/8-(41 x^3)/12+(3 x^4)/8,83/128,0.6484375000}
{5,-35+(5033 x)/60-(567 x^2)/8+(111 x^3)/4-(41 x^4)/8+(11 x^5)/30,237/128,1.851562500}
{6,230-(16961 x)/30+(37933 x^2)/72-(11653 x^3)/48+(8537 x^4)/144-(589 x^5)/80+(53 x^6)/144,-(3669/1024),-3.583007813}
{7,-1624+(445388 x)/105-(1549387 x^2)/360+(539341 x^3)/240-(95287 x^4)/144+(8887 x^5)/80-(7151 x^6)/720+(103 x^7)/280,13461/512,26.29101563}
{8,13209-(6060127 x)/168+(56378641 x^2)/1440-(10802551 x^3)/480+(43757951 x^4)/5760-(12461 x^5)/8+(549883 x^6)/2880-(14421 x^7)/1120+(2119 x^8)/5760,-(5501853/32768),-167.9032288}
{9,-120287+(860791079 x)/2520-(3953981713 x^2)/10080+(22110414181 x^3)/90720-(526937449 x^4)/5760+(46913461 x^5)/2160-(9462317 x^6)/2880+(9289093 x^7)/30240-(217549 x^8)/13440+(16687 x^9)/45360,42222967/32768,1288.542694}
{10,1214674-(642325433 x)/180+(27086967281 x^2)/6300-(517096898113 x^3)/181440+(23497209427 x^4)/20160-(2143606915 x^5)/6912+(3153940991 x^6)/57600-(110070719 x^7)/17280+(6308927 x^8)/13440-(14417579 x^9)/725760+(16481 x^10)/44800,-(2907506455/262144),-11091.25692}
{11,-13469896+(282585480647 x)/6930-(430248854827 x^2)/8400+(6508673943193 x^3)/181440-(119110424663 x^4)/7560+(667294565603 x^5)/145152-(52735532429 x^6)/57600+(2522393789 x^7)/20160-(155221343 x^8)/13440+(499542371 x^9)/725760-(28924153 x^10)/1209600+(1468457 x^11)/3991680,6980687065/65536,106516.8314}
{12,162744945-(222664046271 x)/440+(7121098228717 x^2)/10800-(41844788861927 x^3)/86400+(2460517309163759 x^4)/10886400-(6891019525375 x^5)/96768+(680854337685337 x^6)/43545600-(279135072763 x^7)/115200+(547768029073 x^8)/2073600-(1370534447 x^9)/69120+(42483796219 x^10)/43545600-(753807931 x^11)/26611200+(16019531 x^12)/43545600,-(4732718849353/4194304),-1.128368103*10^6}
{13,-2128047987+(271430090069503 x)/40040-(43362116234713 x^2)/4752+(6642058057364639 x^3)/950400-(37418918675442277 x^4)/10886400+(2796783305681173 x^5)/2419200-(2397936803100967 x^6)/8709120+(114317105939191 x^7)/2419200-(12168423536459 x^8)/2073600+(180610549367 x^9)/345600-(283350260293 x^10)/8709120+(35771612707 x^11)/26611200-(15859335683 x^12)/479001600+(63633137 x^13)/172972800,54831333572045/4194304,1.307280864*10^7}
{14,29943053062-(319417242134749 x)/3276+(393342556439310113 x^2)/2910600-(76567449093702827 x^3)/712800+(6624721229009437643 x^4)/119750400-(47632136694776917 x^5)/2419200+(72946970408344967 x^6)/14515200-(3895040376129719 x^7)/4147200+(79047421468324999 x^8)/609638400-(54758877769571 x^9)/4147200+(85265811978673 x^10)/87091200-(18415340851 x^11)/358400+(577397295443 x^12)/319334400-(158828309957 x^13)/4151347200+(2467007773 x^14)/6706022400,-(5517754573749165/33554432),-1.644419007*10^8}
{15,-451123462672+(67512844627464781 x)/45045-(40356589305665798357 x^2)/18918900+(379015165990587837577 x^3)/216216000-(112620449550302245477 x^4)/119750400+(252732070009773097847 x^5)/718502400-(4151338577382222403 x^6)/43545600+(3575430639143017631 x^7)/186624000-(195745054926250289 x^8)/67737600+(85756141531559659 x^9)/261273600-(486334839668323 x^10)/17418240+(25163220248729183 x^11)/14370048000-(75238450631111 x^12)/958003200+(89011050168379 x^13)/37362124800-(1272153674941 x^14)/29059430400+(34361893981 x^15)/93405312000,37446858098739135/16777216,2.232006675*10^9}
{16,7245940789073-(1359549666317546843 x)/55440+(1665123679437862741457 x^2)/46569600-(183323636693452497903299 x^3)/6054048000+(1469163482186365221977431 x^4)/87178291200-(2366399994198347097329 x^5)/359251200+(46826289630262014389 x^6)/24883200-(1052313210527999927711 x^7)/2612736000+(213027478302403556459 x^8)/3251404800-(14913032249771605457 x^9)/1828915200+(378418224298648607 x^10)/487710720-(402278607448313291 x^11)/7185024000+(658986069102853 x^12)/218972160-(335376314410787 x^13)/2874009600+(27567989527297 x^14)/8941363200-(129887959248197 x^15)/2615348736000+(15549624751 x^16)/42268262400,-(69828647039070145245/2147483648),-3.251649767*10^10}
{17,-123604151490591+(579317032462666755641 x)/1361360-(54864766413333611745301 x^2)/86486400+(10028908779034404119686151 x^3)/18162144000-(138465428409646411299784093 x^4)/435891456000+(21113241681262192097743231 x^5)/163459296000-(74089942236114689755967 x^6)/1916006400+(251391850450227291156683 x^7)/28740096000-(504082131984754695881 x^8)/331776000+(467408212333386341753 x^9)/2286144000-(52115368516292115317 x^10)/2438553600+(28993695623116736327 x^11)/16765056000-(4108910966971393637 x^12)/38320128000+(933226294327961507 x^13)/186810624000-(2808811464990533 x^14)/16605388800+(10272694376985091 x^15)/2615348736000-(129994862918359 x^16)/2324754432000+(8178130767479 x^17)/22230464256000,1086809948534406710475/2147483648,5.060853197*10^11}
{18,2231697509543362-(3188233375801061607493 x)/408408+(10786527177922122290283661 x^2)/908107200-(38512384961362443174387571 x^3)/3632428800+(373640469395388731308684669 x^4)/59439744000-(6927762267466540163452521643 x^5)/2615348736000+(19465979017820380018793598371 x^6)/23538138624000-(64542319673178096146351 x^7)/328458240+(87291576859438774799326103 x^8)/2414168064000-(1523802716374561537587677 x^9)/292626432000+(517849814802894914290613 x^10)/877879296000-(8458462196309648648149 x^11)/160944537600+(26525031638613812771381 x^12)/7242504192000-(589606011929998366351 x^13)/2988969984000+(45896942681626161163 x^14)/5706215424000-(501940184278366147 x^15)/2092278988800+(38718157452340007 x^16)/7846046208000-(14829677125028593 x^17)/237124952064000+(138547156531409 x^18)/376610217984000,-(143981815027423756093635/17179869184),-8.380844667*10^12}
{19,-42519034050101744+(1464240507393101738410753 x)/9699690-(258131132034885503570056901 x^2)/1102701600+(857487960819503527422310837 x^3)/4009824000-(6087167689932311988906043769 x^4)/46702656000+(148389119509356408514528592087 x^5)/2615348736000-(433071493397110157668417626619 x^6)/23538138624000+(23898407610549175338935090071 x^7)/5230697472000-(6936501517986786620701709 x^8)/7838208000+(62213773761277197458472989 x^9)/459841536000-(2061236955680464451428681 x^10)/125411328000+(2559195166511695096958033 x^11)/1609445376000-(885361744120563503318269 x^12)/7242504192000+(14104327056699432268483 x^13)/1902071808000-(3140659110828822892501 x^14)/8966909952000+(672802487229400417 x^15)/53374464000-(1496203730800907521 x^16)/4483454976000+(1451999857330005437 x^17)/237124952064000-(40468364903219737 x^18)/582033973248000+(92079694567171 x^19)/250298560512000,157990833796706803456665/1073741824,1.471404301*10^14}
{20,852495597142800377-(14008964102294687701871803 x)/4564560+(149562539897658522707217800777 x^2)/30875644800-(12085888115355076001965759687 x^3)/2673216000+(44883417392551573656826201682771 x^4)/15878903040000-(2210180058211139429647574187487 x^5)/1743565824000+(40062119563571628583407444018133 x^6)/94152554496000-(3451142137051290051306730529821 x^7)/31384184832000+(954483496055296473187927296127 x^8)/42796615680000-(3307042066776598919704495807 x^9)/919683072000+(77684698944490702270862129 x^10)/167215104000-(103624443815140539285858737 x^11)/2145927168000+(2339950625535963598209771179 x^12)/579400335360000-(3089830238692244867604377 x^13)/11412430848000+(10867554099021629180978717 x^14)/753220435968000-(102921749742912627671 x^15)/170798284800+(384272707331917992979 x^16)/19926466560000-(12727916760447118163 x^17)/27897053184000+(17494324264828917439 x^18)/2328135892992000-(6235268597310551443 x^19)/81096733605888000+(4282366656425369 x^20)/11640679464960000,-(750417552657363912173206755/274877906944),-2.730003153*10^15}

点评

我是故意的,^_^  发表于 2019-2-8 21:20
你用工具用得66的,但我是手算,能算到5阶已经花了九牛二虎之力了,你让我算6阶,我实在是有心无力,不是故意不算的,希望你理解。  发表于 2019-2-8 09:12
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

小黑屋|手机版|数学研发网 ( 苏ICP备07505100号 )

GMT+8, 2024-5-18 11:00 , Processed in 0.055615 second(s), 20 queries .

Powered by Discuz! X3.5

© 2001-2024 Discuz! Team.

快速回复 返回顶部 返回列表