2010 五月的整月存档

正则简单笔记

发表于： 2010年05月28日分类：统计之路 2 Comments 标签: regular expression, 正则表达式, 笔记.

1. 元字符 (Metacharacter)

匹配字符串的开始(位置): ^

匹配字符串的开始(不受处理多行选项影响): \A

匹配字符串的结束(位置): $

匹配字符串的结束(不受处理多行选项影响): \z

匹配字符串的开始/结束(位置): \b

匹配字符串的结束/行尾(不受处理多行选项影响): \Z

匹配数字/字母/下划线/汉字: \w

匹配(阿拉伯)数字: \d

匹配空白符(空格/全角空格/制表符/换行符): \s

匹配横向制表符: \t

匹配竖向制表符: \v

匹配回车: \r

匹配换行符: \n

匹配换页符: \f

匹配任意字符(除换行符\n): .

前述表达式重复任意次: *

2. 元字符转义 (查找元字符本身)

\. \* \\ \^ \$ 等.

3. 字符类 (自定义的元字符): [ ]

e.g. [13579] 匹配1或3或5或7或9中的任意一个数字.

4. 限定符 (规定表达式重复次数)

重复至少0次: *

重复至少1次: +

重复0次或1次: ?

重复n次: {n}

重复n次或更多次: {n,}

重复m到n次: {m,n}

e.g. \d{2,9} 匹配2-9位阿拉伯数字.

5. 分枝: |

p.s. 分枝条件按照从左到右的顺序测试. 如果已满足了某个分枝, 则匹配结束. 故使用分枝条件时, 一定要反复斟酌各个条件的书写顺序.

6. 反义

对于元字符, 把元字符里的字母小写换为大写即表示反义. 对于字符类, [^13579] 匹配除了1或3或5或7或9以外的任意字符. 一个不错的例子:

<a[^>]+> 匹配用尖括号括起来的以a开头的字符串

7. 分组: ( )

和编程写公式时多加括号防止写错有点类似. 不过这里的现实意义强得多, 可以用来写出比较复杂的表达式.

8. 贪婪/懒惰

在正则表达式中, [最先开始的匹配] 拥有最高优先权. 其次才是默认的贪婪匹配原则. 与之对应的懒惰匹配原则:

重复0次或1次: ??

重复至少1次: +?

重复m到n次: {m,n}?

重复n次以上: {n,}?

重复任意次: *?

9. 零宽断言

(1) 零宽度正预测先行断言: (?=exp)

断言此位置之后能匹配表达式exp

e.g. \b\w+(?=ly\b) 匹配以ly结尾的字符串的前面部分.

(2) 零宽度正回顾后发断言: (?<=exp)

断言此位置之前能匹配表达式exp

e.g. (?<=\bpre)\w+\b 匹配以pre开头的字符串的后面部分.

(3) 零宽度负预测先行断言: (?!exp)

断言此位置之后不能匹配表达式exp

e.g. \b((?!exe)\w)+\b 匹配不包含连续字符串exe的字符串.

(4) 零宽度负回顾后发断言: (?<!exp)

断言此位置之前不能匹配表达式exp

e.g. (?<!010)\d{8} 匹配前面不是010的8位数字.

10. 注释

(?#注释内容)

复化Gauss-Legendre求积

发表于： 2010年05月28日分类：数字之美 0条评论标签: 勒让德, 复化, 数值积分, 高斯.

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

一般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() 如下:

^?Download compgauss.m

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点值时, 需要解一系列非线性方程组. 其解析解较难求得. 而取数值解作为公式中的参数, 不如解析解理想.

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

发表于： 2010年05月24日分类：数字之美 0条评论标签: MATLAB, Romberg's Method, 时间换空间, 空间换时间, 算法复杂度, 龙贝格算法.

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

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

……

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

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

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

^?Download romberg.m

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)}$ 的表格, 以方便完成笔头的作业. 懒到一定程度了是不是?
继续阅读’时间换空间/空间换时间？’