博文
Romberg算法求数值积分的原理与实现(2006-10-18 22:52:00)
摘要：数值积分实际上都是基于插值的，我们算的都是离散的一些个点与点对应的函数值，于是要求连续的积分是不可能的。
记等距分点a<x0<x1<...<b，距离h=(b-a)/n，n为分段数，设它的插值函数存在且为L(x)，我们用SL(x)dx的积分代替Sf(x)dx的积分值，这就是数值积分的最基本原理。
Newton-Cores公式，给出上面的积分为：
In=(b-a)∑[ C(n,k)f(xk) ]
其中C(n,k)是Cores系数，设x=a+th，则有
C(n,k)=(-1)^(n-k)/(nk!(n-k)! S∏(t-j)dt j=0 to n,j!=k
当n=1时，算出C(1,0)=C(1,1)=1/2，即熟悉的梯形公式
I1=(b-a)(f(a)+f(b))/2
当n=2时，算出C(2,0)=1/6,C(2,1)=4/6,C(2,2)=1/6，这时就是Simpson公式，记得原来数学分析书上有提到过，它是按二次曲线插值的，有比较好的效率
I2=(b-a)(f(a)+4f((a+b)/2)+f(b))/6
这样一直算下去，会有更好的公式，且当n为偶数时，具有n+1次代数精度，但是当n>=8时，系数出现负项，将会不稳定，因此实际中只使用低次的Newton_Cores公式。
解决的办法有几种，一种是进行复化的求积公式，跟低次分段插值一样，将积分区间分成若干小段，每小段做低次积分，再求和；另一种办法就是用所谓的变步长积分。
其实所谓变步长法，就是我们熟知的递推法，在这里将看到递推法的妙处。我们拿到一个积分，和被积区间(a,b)，如何确定h呢？递推，先用较大的h试，如果不满意就将它细分，直到满意为止。
如果一开始只有一个分点，再将它们二分得四个，一直二分下去，那么原梯形公式将得到递推梯形公式，它是从n个分段二分到2n个的递推关系：
T2n=Tn/2+h/2∑f(xk+1/2) [1]
右边求和部分是新插入的结点，于是就很方便程序控制，误差控制就考虑|T2n-Tn|的值，事后估计误差。但是直接用梯形法，收敛速度很慢，进一步优化发现，梯形法的余项展开为：
T(h)=I+a1h^2+a2h^4+a3h^6+...
如果h比较小，那T(h)=I+O(h^2)，为了提......
牛顿插值算法与实现(2006-09-15 00:09:00)
摘要：牛顿真是牛，拉格朗日插值法只能算是数学意义上的插值，从插值基函数的巧妙选取，已经构造性的证明了插值法的存在性和惟一性，但是从实现的角度看并不很好，而牛顿很好的解决了这个问题。
牛顿插值是基于下面这些的公式：
f[x0,x1,...xk]=(f[x1,...xk]-f[x0,...xk-1])/(xk-x0)
f[x]=f(x)
f(x)=f[x0]+f[x0,x1](x-x0)+f[x0,x1,x2](x-x0)(x-x1)+...f[x0,...xn](x-x0)...(x-xn-1)+Rn(x)
前两个是均差的递推关系式，而后一个就是牛顿插值公式，其中N(x)=f(x)-Rn(x)，即目标多项式，Rn(x)是n阶插值余项，我们就是用N(x)去近似f(x)。
可以构造这样一个均方差表：
xk f(xk) 一阶均差 二阶均差 ...
x0 f(x0)
x1 f(x1) f[x0,x1]
x2 f(x2) f[x1,x2] f[x0,x1,x2]
...
如果有n个点插值，表会有(n*n)/2+n个表项，如果直接编程会有O(n*n)的空间复杂度，编程时做个简单的改进，不难发现在这个表中只有部分数据有用，对角线（斜行）它们是目标值，用来表示多项式的，左边的两纵行（实际上只需要x一行）以及最底下的一行，表示当前插值的状态。经过改进后只需要O(n)的空间复杂度。
两个过程：
1，新增加一个点时的更新。只须更新最底下一行数据，其递推关系由均差公式给出，最后算出高一队的均差值，需时O(n)
2，插入点完成后如何计算多项式在另外给定点的值N(x)。
由牛顿插值公式，最终的表达式为：
N(x)=f[x0]+f[x0,x1](x-x0)+f[x0,x1,x2](x-x0)(x-x1)+...f[x0,...xn](x-x0)...(x-xn-1)
如果直接将它展开，再算实在麻烦，实际上大可不必这样做，还记得多项式求值的秦九......
误差的积累(2006-09-02 00:22:00)
摘要：这学期开了一门我很喜欢的课程《数值分析》，我感觉它讨论的都是如何控制误差，还有如何提高算法效率。
教科书上的例子，计算定积分I(n)=e^-1Sx^ne^xdx[0,1] 积分区间[0,1]并估计误差。
首先是找到递推公式I(n)=1-nI(n-1)，然后求出I(0)=1-1/e，用泰勒公式以部分和近似代替取值为0.6321，且知道误差限为0.25*10^-4，然后以递推公式算I(1),I(2),...I(n)，并求它们的误差限，结果是求到I(8)的时候出现了负值，从原始积分看不可能为负值，于是就考虑误差从哪里‘积累’起来了？
原来是递推公式有问题，这样的东西以前从来没考虑过，从操作次数和效率上看这个递推式是不错的算法，但从误差角度看很烂，可以算出e(n)=n!*e(0)=n!*0.25*10^-4，阶乘啊，多么恐怖，即使是e(0)有非常大的精确度，在阶乘面前也是弱不禁风，算不了几下，指数效应就出来了，所以整个算法彻底失败。
误差的累积是算法上的问题，所以要控制误差就需要改进算法，做的改进是，将它反过来，化成：
I(n-1)=(1-I(n))/n
一眼看去，反过来不是很好的办法，如果我要求I(n)就要先求I(n+1)，如此没有个头了，其实不然，从原积分看知道I(n)单调减，且有下界，所以它的I-n图是一个无穷趋于0的图像，所以取一个很大的n，比如n=1000，设I(1000)=0，那误差也是很小的，然后反过来递推，由于这里误差没有累积，所以算1000前面的所有值，误差也可以在控制之中，有相当的精准度，而整个算法也是线性的，相当巧妙。......