用mathematica自动推导数值积分的余项

Question

比如要计算$\int_a^b f(x)dx$，如果输入intrmd[ {a,b}, {1,1}/2, {a,b} ]，表示取点a和b处的函数值，加权系数均为1/2，去逼近区间[a,b]上的积分，即$\int_a^b f(x)dx\approx \frac{f(a)+f(b)}{2}$。这是梯形公式，其输出应该为$-\frac{(b-a)^3}{12}f''(\eta)$。

而如果输入intrmd[ {a, (a+b)/2, b}, {1,4,1}/6, {a,b} ]，则表示$\int_a^b f(x)dx\approx \frac{f(a)+4f((a+b)/2)+f(b)}{6}$，这是辛普森公式，其输出应为$-\frac{(b-a)^5}{90}f^{(4)}(\eta)$。

如果输入intrmd[ {$\frac{a+b}{2}-\frac{b-a}{2\sqrt{3}}$, $\frac{a+b}{2}+\frac{b-a}{2\sqrt{3}}$}, {1,1}/2, {a,b} ]，这是两点的高斯求积公式，其输出应为$\frac{(b-a)^4}{135}f^{(4)}(\eta)$。

系数的符号可以不管，输出的函数是几阶导数似乎可以令$f(x)$分别取$1,x,x^2,\dots$来估计（代数精度的定义），但是系数不知该如何计算。可否利用mathematica来代替手工繁琐的推导呢？比如输入intrmd[ Range[0,1,1/30], Total@NestList[RotateRight[#, 3] &, PadRight[{1, 3, 3, 1}, 30 + 1],  10 - 1], {0,1} ]，这是把$[0,1]$分为10等分，每区间内采用3阶牛顿法的数值积分，其结果应该是什么呢？

Answer 1

首先纠正一下你的加权系数给的不对，梯形公式是$\int_a^b f(x)dx\approx \frac{b-a}{2}(f(a)+f(b))$，辛普森公式是$\int_a^b f(x)dx\approx \frac{b-a}{6}(f(a)+4f(\frac{a+b}{2})+f(b))$，其余同理

然后是代码

intrmd[pts_,weights_,{a_,b_}]:=Module[{m},
(*先计算代数精确度m*)
m=NestWhile[#+1&,1,FullSimplify[Integrate[x^#,{x,a,b}]==weights.pts^#]&]-1;
(*根据插值多项式误差估计积分得到余项*)
(Derivative[m+1][f][\[Eta]]/(m+1)! Integrate[Evaluate[Times@@(x-pts)*x^(m+1-Length[pts])],{x,a,b}]//FullSimplify
]

大概的解释就是我之前说的那些这里就先不写了，晚上下课回来有时间的话再补充下吧

这是一些运行的例子

不过这个写法目前只能处理代数精确度m大于给定点数量的情况，对于你给的那个31个点精确才到3阶的情况还处理不了，我晚上回来再看看吧

(*顺便你最后给的那个东西权重也有问题，应该除一个80才对*)

Comments

多个子区间的情况貌似就是加一下……我也是从这里现学现卖的：http://www.uio.no/studier/emner/matnat/math/MAT-INF1100/h11/kompendiet/chap12.pdf

---

嗯确实是加一下就好，但问题是怎么根据给定的pts和weights判断出这是一个分段辛普森积分的累加，这个我想了半天还是没想出来

---

能不能根据系数的周期性来解决呢？