分类目录归档：数字之美

It is from The Book!

不可靠的直觉：Erdős与蒙迪·霍尔问题

发表于 由
2

读过《Erdős的传说》，一度对Erdős的研究风格和1475片paper心向往之，遂购得一本《数字情种》。在这书的第六章发现原来Erdős也做错了这道最初学概率论时做错的问题，顿时倍感欣慰 。。。 侵权手打版在此，请勿下载阅读：

《数字情种：埃尔德什传》第六章 得了一只羊

看完后重现了一下这个模拟，原来确实是这样：

试验次数 / 不换门成功 / 换门成功

10000 / 3350 / 6729

20000 / 6770 / 13475

30000 / 9982 / 19959

40000 /13172 / 26690

50000 / 16844 / 33117

60000 / 20019 / 40175

70000 / 23379 / 46622

80000 / 26645 / 53242

90000 / 30043 / 60100

100000 / 33415 / 66585

erdos

这个故事在Erdos的另一本传记《我的大脑敞开了：数学怪才爱多士》中也有提到。有好事者做了调查，结果只有13%的人做对这个概率问题，故此问题获得了"fool all the people all the time"的荣誉称号，这是后话。

引经据典的延伸阅读： http://en.wikipedia.org/wiki/Monty_Hall_problem

1475片paper！！！paper狂人有木有！！！

Ubuntu下安装UMFPACK的MATLAB接口

发表于 由
回复

求解大型稀疏矩阵线性方程组是一件很痛苦的事情，大型稀疏矩阵求逆是这个问题的特例。UMFPACK采用的是一种针对稀疏矩阵的LU分解方法，直接法，比较耗内存，提供了C/Fortran/MATLAB接口，如果内存够大，可以首先尝试一下，因为直接法较迭代法精度高一些。如果内存不够，一般首先进行预处理，然后选择一种迭代法进行求解，比方说基于Krylov子空间的方法。UMFPACK的Manual写得比较暧昧，下面记录一下其MATLAB接口的安装方式。虽然MATLAB这丫闭源，只叹人在江湖，身不由己 。。。不过这套SparseSuite倒是GPL的，要是真提供R接口就好了，和CSIE@UFL他家比较熟的童鞋，不妨建议一下 。。。

环境：Ubuntu 10.10 + MATLAB R2010b for UNIX

1. 从 http://www.cise.ufl.edu/research/sparse/

下好UMFPACK本身及其依赖包（均属于所谓的 SparseSuite），平行放置于一个目录。理论上这些就足够：

AMD
CAMD
CCOLAMD
CHOLMOD
COLAMD
UFconfig
UMFPACK

从 http://glaros.dtc.umn.edu/gkhome/metis/metis/download

下载metis-4.0，解压后同样放入上述目录。

2. 此时如果在MATLAB中进入UMFPACK/MATLAB目录执行 umfpack_make，会发现不能编译Mathworks他丫的MEX。提示：

Warning: You are using gcc version "4.4.4-14ubuntu5)".  The version currently supported with MEX is "4.3.4".

看来Ububtu 10.10自带的gcc版本太高，我们按照

https://help.ubuntu.com/community/MATLAB

的MEX function一节的1 - 3步先安装一个gcc 4.3.4 。。。

3. root启动MATLAB，执行

getenv('PATH')
setenv('PATH',sprintf('/home/%s/.matlab/bin:%s',getenv('USER'),getenv('PATH')));

再次 getenv('PATH') 发现环境变量已经修改成功，进入UMFPACK/MATLAB目录，执行 umfpack_make 即可。

要确认是否已经编译成功，执行 umfpack_demo 即可运行一个求解demo，还有一些对比内建浮点性能之类的信息，比较无聊。

 

另外，直接法还有一个求解器PARDISO可以选择，但是好像授权很不友好，academical也只能试用30天，简直是让人发指 。。。没有做更多的尝试。

palintiples引发的血案

发表于 由
回复

商务印书馆版《一个数学家的辩白》在pp.75给出了这样一个性质:

四位数中仅有8712和9801是其"反序数"的整数倍:
8712=4x2178, 9801=9x1089
没有比10000更小的数具有这种性质.

这个性质是作为不"严肃"定理的例子给出的. 原始出处是Rouse Ball的《Mathematical Recreations and Essays》. 像这样满足"自身是其'反序数'的整数倍"的整数被称为"palintiples". 我们不妨写一小段程序继续探索更大的"palintiples"(当然要剔除一些琐碎的特殊情况):

?[Copy to clipboard]View Code RSPLUS
nixu=function(x){
options(scipen=20);
a=min(x)
b=max(x)
n=length(x)
numSlot=x
numRev<-function(num){
rem=0
while(num>0){
rem = (rem * 10) + (num %% 10)
num = as.integer(num / 10)
}
return(rem)
}
for(j in 1:n){
numSlot[j]=numRev(numSlot[j])
}
numSlot[]=x/numSlot[]
digits=as.integer(log10(b))
for(k in 0:digits){
numSlot[][numSlot[]==10^k]=0.1
}
numSlot[][seq(a, b)/10==as.integer(seq(a, b)/10)]=0.2
numSlot[]=numSlot[]==as.integer(numSlot[])
x[numSlot[]==1]
}
# x=seq(1, 10000); nixu(x)

受计算资源限制, 我只验证了10亿以内的整数. 结果分布如下:

8712 9801
87912 98901
879912 989901
8799912 9899901
87128712 87999912 98019801 98999901
871208712 879999912 980109801 989999901

下面是其他人给出的, 100亿以内的结果:

8712 9801
87912 98901
879912 989901
8799912 9899901
87128712 87999912 98019801 98999901
871208712 879999912 980109801 989999901
8712008712 8791287912 8799999912 9801009801 9890198901 9899999901

可以看到结果中始终没有出现除了87*12和98*01这两个类型之外的数字. 是不是可以猜想, 扩展到n位时, 87*12和98*01始终将保持这种性质呢?

和数学中众多的"好"问题一样, 连这个问题也是早就有人研究过了. Dan Hoey于1991-1992年间在这里给出了palintiples的一般形式. 并且说明了全部的palintiples并没有"form a regular language"的事实. 他的证明很漂亮. 这里还有9801这个数字一些有意思的其他性质.

后来居然还有一群无聊人士继续对这个问题的扩展进行了深入研究. 2007年, 有人在《Mathematics Magazine》上发了一片内容和回顾兼有而标题有点幽默的:《Digit Reversal Without Apology》, 在这里可以看到. 此文分类Number Theory, 我看着看着就睡着了. 文章最后给出了几个Open Questions, 有能力的同学可以挑战一下.

一个随便的quote都能引发血案, 哈代无愧为一代纯数学宗师.

《一个数学家的辩白》那些牛x闪闪的片段

发表于 由
1

Hardy我买到的是商务印书馆的那一版, 正文字体好像是方正博雅宋, 较一般书宋要扁一些, 行气足一些. 可能是为了缩减成本, 版面整体看起来很囧, 下次最好用大一点的纸张, 以免出现不能容忍的, pp. 57中英对照诗文里单字换行的扯淡效果.

薄薄的一本小书, 只有约摸100页, C. P. Snow所撰前言和正文几乎一样长, 正文由29节构成. 在我看来, 哈代几乎是以一种那个领域中精神领袖的高度写就了全文. 读完他的"辩白", 更是能体会到他作为一个纯粹数学研究者的强烈自尊和无可言说的寂寞. 同时, 字里行间也浸透了他在年事已高无法继续做出创造性工作之时的伤感之情, 让人想起了简媜的: "深情即是一桩悲剧 必得以死来句读." 哈代以一种非平凡的方式, 为(纯)数学和自己进行了"辩白", 同时夹杂了一些对数学, 对美丑标准的讨论. 有志于进行纯数学研究的同学们可以一读, 此书也可以作为那些做着更偏重"应用"数学的同仁们的一本参考读物. 遗憾的是里面对哈代和Ramanujan的合作情况介绍得还是太少, 这也许需要去那本《The Man Who know Infinity》中寻找了.

摘录一些片段于此:

政治家轻视政论家, 画家瞧不起艺术评论家, 生理学家、物理学家或者数学家常常具有同感; 再没有比对解释者所怀有的藐视更深刻, 或者从总体上看更合理的了. 解释、评判、欣赏都是二流头脑干的事情.

为什么对数学进行严肃的研究确实是值得的? 什么能合理地证明数学家生涯的价值? 从大体上来看, 我的答案类似于数学家常被期待做出的那种: 我认为它是值得做的, 这可以有大量证明.

没有数学家会允许自己忘记这一点: 数学比任何其他的艺术或科学都更加是年轻人的游戏.

我所知道的每个真正有天分的年轻数学家都对数学忠心不二, 而且原因不是由于他们缺少抱负, 相反倒是由于他们充满抱负. 他们都认识到, 如果有成名成家的道路的话, 它就着落在数学上面.

抱负是高尚的激情, 它可以合理地以许多形式显现出来. 在阿提拉和拿破仑的抱负中有一些高尚的东西: 但最高尚的志向是在死后留下某些具有永恒价值的东西.

可能会有很多值得高度尊敬的动机促使人们开展研究工作, 但其中有三个比其他的更为重要. 第一个(没有它则其他的都谈不上)是智力上的好奇心, 知道真理的欲望.

在所有的成就中, 数学成就的生命力乃是最持久的. ... 阿基米德(Archimedes)在埃斯库罗斯(Aeschylus)被忘却后, 仍将长留在人们的记忆里, 因为语言会灭亡, 而数学思想则不会.

一位数学家就像一位画家或诗人, 是模式(pattern)的创造者. 如果他的模式比画家或诗人的模式的生命更加长久的话, 那是因为他的模式使用思想(idea)所造就的.

... 这些思想, 就像色彩或者字词一样, 必须以和谐的方式统一起来. 优美性是第一道检验标准: 这个世界没有为丑陋数学准备长久的地盘.

最好的数学既是严肃的, 又是优美的 ... 一个数学定理的"严肃性", 不在于它在实践上的效果——那往往是可以忽略不计的, 而在于它所关联着的数学思想是"重要的". ... 一个严肃的数学定理, 联系着重要的数学思想, 很可能导致数学, 乃至其他科学的重要进展.

任何一流的数学定理都有某种程度的普遍性, 但若普遍性过强, 则定理无可避免地趋向平淡乏味. ... "被恰当的特殊性所限制了的巨大的普遍性, 才是富有成果的概念."

任何名副其实的职业数学家生涯的价值, 是不可能在其工作的"有用性"的基础上得到证明的.
(评: 从这个角度, 炸药奖如果真的设立数学奖, 那还真是对数学的侮辱.)

我相信数学实在是外在于我们的, 我们的作用是发现或观察它, 我们证明了的, 我们卖弄般地描述为我们"创造"的定理, 只不过是我们的观察笔记.

数学家和物理学家在立场上的差别, 也许比通常所想的要小. 在我看来, 其中最重要的差别似乎是: 数学家与实在有更直接的接触. ... 实在论观点对于数学实在比对于物理实在更有可能是正确的, 因为数学对象与它们看上去的样子是那样的一致. ... 现代物理学或许最为符合某种唯心论哲学的框架 ... 对我来说, 纯粹数学像是岩石, 所有的唯心论都在它上面碰得头破血流.
(评: 某些物理真的相当唯心)

正是数学中常识和枯燥的东西, 对实际生活才具有重要性.

如果我想要的是合情合理地舒适和快乐的生活, 那我的选择就是正确的. 事务律师和股票经纪人和出版家常常会过上舒适和快乐的生活, 但很难看出世界因他们的存在而变得更加丰富多彩了.

如果在伦敦的一个纪念柱上刻一幅雕像, 我们是宁愿让这个柱子高得使雕像看不见、还是愿意让柱子低得足以使雕像的容貌可以被辨认出来呢? 我将挑选第一种可能, 而可以想见, 斯诺博士会挑选第二种.

复化Gauss-Legendre求积

发表于 由
回复

高斯型求积公式是数值积分中一个比较成熟的想法. 速度快, 精度高, 公式优雅.

一般Gauss-Legendre求积公式

对于一般的Gauss型求积公式
\int_a^b f(x) \rho (x) dx \approx \sum_{k=0}^{n} A_k f(x_k)
取权函数\rho (x)=1, 积分区间设定为\left[-1, 1 \right], 则得到一般的Gauss-Legendre求积公式
\int_{-1}^{1} f(x) dx \approx \sum_{k=0}^{n} A_k f(x_k).
取Legendre多项式的零点作为Gauss点, 通过计算得到Gauss点个数n=2和n=3时的求积公式
\int_{-1}^{1} f(x) dx \approx f(-\frac{1}{\sqrt{3}})+f(\frac{1}{\sqrt{3}})
以及
\int_{-1}^{1} f(x) dx \approx \frac{5}{9} f(-\frac{\sqrt{15}}{5})+\frac{8}{9}f(0)+\frac{5}{9}f(\frac{\sqrt{15}}{5}).

变量替换方法

上述分析中, 积分区间固定为\left[-1, 1\right], 实际应用时做变量替换
x=\frac{a+b}{2}+\frac{b-a}{2}t.
将被积区间\left[a, b \right]化为\left[-1, 1 \right].

复化Gauss-Legendre求积公式

将被积区间m等分, 记h=\frac{b-a}{m}, x_k=a+kh, k=0, 1, 2 \ldots, m. 作变换
x=\frac{x_k+x_{k-1}}{2}+\frac{h}{2} t

在每个小区间上应用Gauss-Legendre公式, 累加即得复化Gauss-Legendre求积公式
\int_a^b f(x) dx=\frac{h}{2}\sum\sum A_k f(x_k).

不妨设
G(t)=f(\frac{x_k+x_{k+1}}{2}+\frac{h}{2} t)

则有:

Gauss点个数n=2时,
\int_a^b f(x) dx \approx\frac{h}{2}\sum G(-\frac{1}{\sqrt{3}})+G(\frac{1}{\sqrt{3}}),

Gauss点个数n=3时,
\int_a^b f(x) dx \approx \frac{h}{2}\sum \frac{5}{9}G(-\frac{\sqrt{15}}{5})+\frac{8}{9}G(0)+\frac{5}{9}G(\frac{\sqrt{15}}{5}).

总结复化Gauss-Legendre求积过程如下:

1. 分割区间, 记录区间端点值；
2. 通过查表或求解非线性方程组, 在所有小区间上, 将Gauss系数和Gauss点的值代入变量替换后的公式；
3. 将所有区间的结果累加, 即得到整个区间上的积分近似值.

针对Gauss点个数n=2和n=3的复化Gauss-Legendre求积公式编写的一个简单的MATLAB函数 compgauss() 如下:

?[Copy to clipboard]View Code SCILAB
function [ ] = compgauss(a, b, n)
% Composite Gauss Integration
% Equation Type: n=2, n=3
% Coded by Nan.Xiao 2010-05-25
% Step.1 Divide Interval
% Step.2 Calculate
% Step.3 Sum Results
format long
f = @(x) exp(x).*sin(x);
h=(b-a)/n;
xk=zeros(n+1,1);
xk(1,1)=a;
xk(n+1,1)=b;
fk1=zeros(n,1);
fk2=zeros(n,1);
for i=1:n-1
    xk(i+1,1)=a+h*i;
end
for j=1:n
    fk1(j)=f((xk(j)+xk(j+1))/2+(h/2)*(-1/sqrt(3)))+...
        f((xk(j)+xk(j+1))/2+(h/2)*(1/sqrt(3)));
end
for r=1:n
    fk2(r)=(5/9)*f((xk(r)+xk(r+1))/2+(h/2)*(-sqrt(15)/5))+...
        (8/9)*f((xk(r)+xk(r+1))/2+(h/2)*(0))+...
        (5/9)*f((xk(r)+xk(r+1))/2+(h/2)*(sqrt(15)/5));
end
mysum1=h*sum(fk1)/2;
mysum2=h*sum(fk2)/2;
disp('Result of 2 Nodes:')
disp(mysum1);
disp('Result of 3 Nodes:')
disp(mysum2);
end

总结

1. Gauss求积公式较一般的机械求积公式的进步之处, 在于其针对插值型求积公式进行了改进. Gauss公式将插值节点设为未知, 成功地将代数精度由n次提高到2n+1次. 计算过程简单, 速度快, 达到要求精度所需步骤较少.

2. (复化) Gauss 求积的一个问题在于, 增加Gauss点个数, 继续求解Gauss系数和Gauss点值时, 需要解一系列非线性方程组. 其解析解较难求得. 而取数值解作为公式中的参数, 不如解析解理想.

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

发表于 由
回复

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

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

......

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

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

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

?[Copy to clipboard]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)}的表格, 以方便完成笔头的作业. 懒到一定程度了是不是?
继续阅读 →