时间换空间/空间换时间？

发表于： 2010年05月24日 分类：数字之美 0条评论 标签: , , , , , .

经常看到网上有人在讨论算法时用时间换空间/空间换时间的说法. 可惜神人们使用的例子都很专业, 再加上一两句晦涩而潇洒的写法, 我等凡人基本上就很难看懂了. 今天在实现数值积分中的龙贝格 (Romberg) 算法时, 算是遇到了一个这方面有点关系的问题, 简单记录一下.

首先上一段800字的原理说明:

……

由于作者心里很清楚读者将会人肉Skip上面一段, 所以这里直接人肉省略其中的797字.

如果执意想看一下, 参考 Wikipedia 好了.

下面是一个很好懂的 MATLAB 实现:

?Download romberg.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63

function [ ] = romberg(aa, bb, kk, iter)
% Romberg Integral Method
% Coded by Nan.Xiao 2010-05-23
% Step 1. Calc TrapeStep
% Step 2. Calc Romberg Method
% Step 3. Extract Diagnal & Output
function [ ] = trapestep(a, b, k)
    format long
%     function y = f(x)
%         y=(2/sqrt(pi))*exp(-x);
%     end
%     Examples Can Also Be:
%     function y = f(x)
%         y=x*sin(x);
%     end
%     Or
    function y = f(x)
        y=x*sqrt(1+x.^2);
    end
k1=2^(k-1);
temp=zeros(k,k1);
mysum=zeros(k,1);
t=zeros(k+1,1);
for i=1:k
    for j=1:(2^(i-1))
    temp(i,j)=f(a+(((2*j-1)*(b-a))/(2^i)));
    end
end
for p=1:k
    mysum(p,1)=sum(temp(p,:));
end
t(1,1)=0.5*(f(a)+f(b));
for q=2:k+1
    t(q,1)=(0.5*t(q-1,1))+(((b-a)*mysum(q-1,1))/(2^(q-1)));
end
end
    format long
trapestep(aa, bb, kk);
diagonal=zeros(kk,1);
tt=zeros(kk,kk);
for v=1:kk
    tt(v,1)=t(v,1);
end
for r=2:kk
for s=2:r
    tt(r,s)=((4^(s-1)*tt(r,s-1))-tt(r-1,s-1))/(4^(s-1)-1);
end
end
save('tt.txt', 'tt', '-ascii', '-double'); % Save the Table
% Extract Table Diagonals
for w=1:kk
    diagonal(w,1)=tt(w,w);
end
disp('Romberg Table Diagonals:')
disp(diagonal);
% Output Final Result
cc=2;
while abs(diagonal(cc,1)-diagonal(cc-1,1))>=iter || abs(diagonal(cc+1,1)-diagonal(cc,1))>=iter
    cc=cc+1;
end
disp('The Final Result is:');
disp(diagonal(cc,1));
end

与网上其他程序不同, 我在这里没有做太多的判断, 将二分的次数k作为一个参数手动设定, 是为了观察计算的具体过程. 同时, 为了方便移植, 仅使用最基本的矩阵操作语句, 没有使用 MATLAB 中的 feval() 函数, 要改变被积函数, 直接修改函数体. 最终目的是, 我要取得计算产生的含有T_m^{(k)}的表格, 以方便完成笔头的作业. 懒到一定程度了是不是?

如果我们把前两个作为Example的被积函数注释掉, 计算第三个积分\int_0^3 x\sqrt{1+x^2}\,dx, 设定iter=1e-6, 将 k 值提高至23, 输入语句 romberg(0, 3, 23, 1e-6), 将会华丽地出现如下错误：

Out of memory. Type HELP MEMORY for your options.

这是因为, 变步长梯形法的递推公式

T_{2n}=\frac{1}{2}T_n+\frac{h}{2}\sum_{k=0}^{n-1}f(x_{k+\frac{1}{2}})

在懒人版的 trapestep() 函数中实现时, 将直接去生成一个 k 行, 2^(k-1) 列的矩阵. 而 2^(23-1)=4194304, MATLAB 觉得似乎没那么大的内存来生成一个有着四百万列的矩阵(64位版本未测).

好, 大不了不生成了, 改掉上面完全不管内存的写法, 特别是改写 25-27 行的 for 循环:

?Download romberg2.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59

function [ ] = romberg2(aa, bb, kk, iter)
% Romberg Integral Method
% Coded by Nan.Xiao 2010-05-23
% Step 1. Calc TrapeStep
% Step 2. Calc Romberg Method
% Step 3. Extract Diagnal & Output
function [ ] = trapestep2(a, b, k)
format long
      function y = f(x)
          y=(2/sqrt(pi))*exp(-x);
      end
%     Examples Can Also Be:
%     function y = f(x)
%         y=x.*sin(x);
%     end
%     Or
%     function y = f(x)
%         y=x*sqrt(1+x.^2);
%     end
temp=zeros(k,1);
t=zeros(k+1,1);
for i=1:k
    for j=1:(2^(i-1))
    temp(i,1)=temp(i,1)+f(a+(((2*j-1)*(b-a))/(2^i)));
    end
end
mysum=temp;
t(1,1)=0.5*(f(a)+f(b));
for q=2:k+1
    t(q,1)=(0.5*t(q-1,1))+(((b-a)*mysum(q-1,1))/(2^(q-1)));
end
end
    format long
trapestep2(aa, bb, kk);
diagonal=zeros(kk,1);
tt=zeros(kk,kk);
for v=1:kk
    tt(v,1)=t(v,1);
end
for r=2:kk
for s=2:r
    tt(r,s)=((4^(s-1)*tt(r,s-1))-tt(r-1,s-1))/(4^(s-1)-1);
end
end
save('tt.txt', 'tt', '-ascii', '-double'); % Save the Table
% Extract Table Diagonals
for w=1:kk
    diagonal(w,1)=tt(w,w);
end
disp('Romberg Table Diagonals:')
disp(diagonal);
% Output Final Result
cc=2;
while abs(diagonal(cc,1)-diagonal(cc-1,1))>=iter || abs(diagonal(cc+1,1)-diagonal(cc,1))>=iter
    cc=cc+1;
end
disp('The Final Result is:');
disp(diagonal(cc,1));
end

理论上只要CPU时间足够, 这样做基本上可以得出正确结果. 我们需要的计算过程保存在 tt.txt 中.

写这个小函数的时候, 还遇到了一个有意思的事情. 在使用 while 语句判断哪个是最终积分结果时, 不仅要与前一项作差, 还要与后一项作差, 否则, 一旦遇到前两项之差较小, 但与真实值相差很大的情况时将判断失误. 试一下第二个被积函数就知道了. 计算过程大致是这样的:


-0.000000000000001
0.000000000000001
-7.018385351885766
-6.266954014124826
-6.283266463460164
-6.283185210826781
-6.283185307207264
-6.283185307179584

......

如果作的是前一种简单判断, 显然将输出 0.000000000000001 作为最终积分结果.

如果你有兴趣的话, 这里附上一份 Dirty Work.

写程序是一门艺术, 而不仅仅是一门技术(又有几人能够达到大神的Literate Programming境界呢？), 在这方面, 用一句名言来形容一下自己, 仍然是 Too simple, sometimes naive.

0 Responses to “时间换空间/空间换时间？”

本文的Feed Trackback地址
  • 暂无评论

留下回复

:wink: :-| :-x :twisted: :) 8-O :( :roll: :-P :oops: :-o :mrgreen: :lol: :idea: :-D :evil: :cry: 8) :arrow: :-? :?: :!: