博文

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),为了提......

阅读全文(14377) | 评论:7 | 复制链接

牛顿插值算法与实现(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)
如果直接将它展开,再算实在麻烦,实际上大可不必这样做,还记得多项式求值的秦九......

阅读全文(17684) | 评论:3 | 复制链接

误差的积累(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前面的所有值,误差也可以在控制之中,有相当的精准度,而整个算法也是线性的,相当巧妙。......

阅读全文(5207) | 评论:1 | 复制链接