Category Archives: Mathematics

拉格朗日四平方和定理

拉格朗日四平方和定理

每个正整数均可表示为四个整数的平方和。

Every positive integer is the sum of four squares.

例如:

  • 1=1^{2}+0^{2}+0^{2}+0^{2}
  • 2 = 1^{2}+1^{2}+0^{2}+0^{2}
  • 7 = 2^{2}+1^{2}+1^{2}+1^{2}

证明:可以直接验证如下恒等式

(x_{1}^{2}+x_{2}^{2}+x_{3}^{2}+x_{4}^{2})\cdot(y_{1}^{2}+y_{2}^{2}+y_{3}^{2}+y_{4}^{2}) = z_{1}^{2}+z_{2}^{2}+z_{3}^{2}+z_{4}^{2},其中

\begin{cases} z_{1}=x_{1}y_{1}+x_{2}y_{2}+x_{3}y_{3}+x_{4}y_{4} \\ z_{2}=x_{1}y_{2}-x_{2}y_{1}-x_{3}y_{4}+x_{4}y_{3} \\ z_{3}=x_{1}y_{3}-x_{3}y_{1}+x_{2}y_{4}-x_{4}y_{2} \\ z_{4}=x_{1}y_{4}-x_{4}y_{1}-x_{2}y_{3}+x_{3}y_{2}\end{cases}

由于 1 与 2 都明显满足这个定理,那么只需要考虑大于 2 的正整数。而这些正整数都可以分解成素数的乘积,因此,只需要证明该定理对所有的素数成立,则使用以上恒等式就可以得到最终的结论。假设 p 是一个奇素数。

由于 \{a^{2}:a\in\{0,1,\cdots,(p-1)/2\}\} 里面有 (p+1)/2 个不同的同余类,\{-b^{2}-1: b\in \{0,1,\cdots,(p-1)/2\}\} 里面也有 (p+1)/2 个不同的同余类,但是素数 p 的同余类只有 p 个,因此存在正整数 a,b\in \{0,1,\cdots, (p-1)/2\} 满足 a^{2}\equiv -b^{2}-1 (\mod p)。也就是说 a^{2}+b^{2}+1^{2}+0^{2}\equiv 0(\mod p)。令 n\in\mathbb{Z} 满足 np=a^{2}+b^{2}+1,则有 p\leq np\leq 2(p-1)^{2}/4+1<p^{2}。于是,1\leq n<p

因此存在一个 1\leq n<p 使得 np = a^{2}+b^{2}+1^{2}+0^{2} 是四个整数的平方和。于是必定存在一个最小的正整数 m 使得 1\leq m\leq n<p 使得 mp 为四个整数的平方和,不妨设为 mp=x_{1}^{2}+x_{2}^{2}+x_{3}^{2}+x_{4}^{2}

Claim. m=1

proof of the claim. 反证法,假设 1<m\leq n<p 成立。令 y_{i}=x_{i}(\mod m) 对于 i\in\{1,2,3,4\} 成立,并且 -m/2<y_{i}\leq m/2。因此,y_{1}^{2}+y_{2}^{2}+y_{3}^{2}+y_{4}^{2}\equiv(x_{1}^{2}+x_{2}^{2}+x_{3}^{2}+x_{4}^{2})\equiv mp \equiv 0(\mod m)。令 mr = y_{1}^{2}+y_{2}^{2}+y_{3}^{2}+y_{4}^{2}。因此,mr\leq 4(m/2)^{2}=m^{2}

如果 r =m,通过以上不等式得知 r=m 等价于 y_{i}=m/2 对于 i\in\{1,2,3,4\} 都成立。此时,mp = x_{1}^{2}+x_{2}^{2}+x_{3}^{2}+x_{4}^{2}\equiv 4(m/2)^{2} \equiv 0 (\mod m^{2})。因此,pm 的倍数,这与 p 是素数,m>1 矛盾。所以,r<m 成立。i.e. 1\leq r<m\leq n<p 成立。

进一步地,(mp)\cdot(mr) = (x_{1}^{2}+x_{2}^{2}+x_{3}^{2}+x_{4}^{2})\cdot(y_{1}^{2}+y_{2}^{2}+y_{3}^{2}+y_{4}^{2}) = z_{1}^{2}+z_{2}^{2}+z_{3}^{2}+z_{4}^{2},这里的 z_{i} 正如恒定式里面所定义的。由于 y_{i}\equiv x_{i}(\mod m),并且 \sum_{i=1}^{4}x_{i}^{2}\equiv 0(\mod m)。因此,z_{i}\equiv 0(\mod m) 对于 i\in\{1,2,3,4\} 都成立。所以,z_{i}=w_{i}mw_{i}\in\mathbb{Z} 对于 i\in\{1,2,3,4\} 都成立。通过 (mp)\cdot(mr) = \sum_{i=1}^{4}z_{i}^{2} 可以得到 pr=\sum_{i=1}^{4}w_{i}^{2} 成立。但是,1\leq r<m 这与 m 的最小性假设矛盾了。

因此,m=1,Claim 证明完毕。

于是,对于所有的奇素数,都可以表示为四个整数的平方之和。根据之前的分析,可以得到对于所有的正整数,都可以表示为四个整数的平方之和。Lagrange 定理证明完毕。

参考文献

  1. GTM 164, Additive Number Theory, Melvyn B.Nathanson, 1996.

传染病的数学模型

近期,国内的疫情闹得沸沸扬扬,很多省市自治区都出现了流感的患者。回想起之前在学校的时候曾经研究过微分方程和动力系统,于是整理一下相关的数学模型,分享给各位读者。笔者并不是研究这个领域的专家,并且这篇文章只是从微分方程角度出发,分析方程的性质,不一定适用于真实环境,而且真实环境比这个也复杂得多。

关于传染病的数学模型,在许多年前数学界早已做过研究,根据传染病的传播速度不同,空间范围各异,传播途径多样,动力学机理等各种因素,对传染病模型按照传染病的类型划分为 SI,SIR,SIRS,SEIR 模型。如果是按照连续时间来划分,那么这些模型基本上可以划分为常微分方程(Ordinary Differential Equation),偏微分方程(Partial Differential Equation)等多种方程模型;如果是基于离散的时间来划分,那么就是所谓的差分方程(Difference Equation)。

在本文中,将会主要介绍常微分方程中的一些传染病数学模型。在介绍方程之前,首先要介绍一些常用的符号:在时间戳 t 上,可以定义以下几种人群:

  • 易感者(susceptible):用符号 S(t) 来表示;
  • 感染者(infective):用符号 I(t) 来表示;
  • 康复者(Recoverd):用符号 R(t) 来表示;

其次,在时间戳 t 上,总人口是 N(t)=S(t)+I(t)+R(t)。如果暂时不考虑人口增加和死亡的情况,那么 N(t)\equiv N 是一个恒定的常数值。

除此之外,

  • r 表示在单位时间内感染者接触到的易感者人数;
  • 传染率:\beta 表示感染者接触到易感者之后,易感者得病的概率;
  • 康复率:\gamma 表示感染者康复的概率,有可能变成易感者(可再感染),也有可能变成康复者(不再感染)。

在进行下面的分析之前,先讲一个常微分方程的解。

Claim. 假设 x=x(t) 是关于 t 的一个方程,且满足 \frac{dx}{dt} + a_{1}x + a_{2}x^{2}=0x(0)=x_{0},那么它的解是:x(t) = \frac{e^{-a_{1}t}}{\frac{1}{x_{0}}-\frac{a_{2}}{a_{1}}(e^{-a_{1}t}-1)}.

Proof. 证明如下:

通过 \frac{dx}{dt}+a_{1}x+a_{2}x^{2}=0 可以得到 -\frac{d}{dt}\bigg(\frac{1}{x}\bigg) + a_{1}\bigg(\frac{1}{x}\bigg)+a_{2}=0;令 y = 1/x,得到 \frac{dy}{dt}-a_{1}y=a_{2}。所以,\frac{d}{dt}(e^{-a_{1}t}y) = a_{2}e^{-a_{1}t},两边积分可以得到 e^{-a_{1}t}y-y_{0}=\bigg(-\frac{a_{2}}{a_{1}}\bigg)(e^{-a_{1}t}-1),其中 y_{0}=1/x_{0}。求解之后得到:x(t) = e^{-a_{1}t}/\bigg(\frac{1}{x_{0}}-\frac{a_{2}}{a_{1}}(e^{-a_{1}t}-1)\bigg)

SI 模型(Susceptible-Infective Model)

在 SI 模型里面,只考虑了易感者和感染者,并且感染者不能够恢复,此类病症有 HIV 等;

SI_model_1
SI Model

其微分方程就是:

\begin{cases}\frac{dS}{dt} = -\frac{r\beta I}{N} S \\ \frac{dI}{dt}=\frac{r\beta I}{N}S \end{cases}

初始条件就是 S(0)=S_{0}I(0) = I_{0},并且 S(t)+I(t)=N 对于所有的 t\geq 0 都成立。

于是,把 S = N - I 代入第二个微分方程可以得到:\frac{dI}{dt} - r\beta I + \frac{r\beta}{N}I^{2}=0。因此根据前面所提到的常微分方程的解可以得到:

I(t) = \frac{NI_{0}}{I_{0}+(N-I_{0})e^{-r\beta t}}.

这个就是所谓的逻辑回归函数,而在机器学习领域,最简单的逻辑回归函数就是 \sigma(x) = 1/(1+e^{-x}) 这个定义。而 I(t) 只是做了一些坐标轴的平移和压缩而已。由于 \lim_{t\rightarrow +\infty}e^{-t}=0,所以,\lim_{t\rightarrow +\infty}I(t) = N,从而 \lim_{t\rightarrow +\infty}S(t) = 0

通过数值模拟可以进一步知道:

SI_model_graph_1
SI model 的数值模拟(一)

简单来看,在 SI 模型的假设下,全部人群到最后都会被感染。

SIS 模型(Susceptible-Infectious-Susceptible Model)

除了 HIV 这种比较严重的病之外,还有很多小病是可以恢复并且反复感染的,例如日常的感冒,发烧等。在这种情况下,感染者就有一定的几率重新转化成易感者。如下图所示:

SIS_model_1
SIS model

其微分方程就是:

\begin{cases} \frac{dS}{dt} = -r \beta S\frac{I}{N} + \gamma I \\ \frac{dI}{dt}=r\beta S \frac{I}{N} - \gamma I \end{cases},其初始条件就是 S(0)=S_{0}I(0)=I_{0}.

使用同样的方法,把 S=N-I 代入第二个微分方程可以得到:\frac{dI}{dt} - (r\beta - \gamma)I + \frac{r\beta}{N}I^{2}=0. 通过之前的 Claim 可以得到解为:

I(t) = \frac{N(r\beta-\gamma)}{r\beta}/\bigg(\bigg(\frac{N(r\beta-\gamma)}{I_{0}r\beta}-1\bigg)e^{-(r\beta-\gamma)t}+1\bigg).

从而可以得到 \lim_{t\rightarrow +\infty} I(t) = N(r\beta - \gamma)/(r\beta)\lim_{t\rightarrow +\infty} S(t) = (N\gamma)/(r\beta). 这个方程同样也是逻辑回归方程,只是它的渐近线与之前的 SI 模型有所不同。

SIS_model_graph_1
SIS model 的数值模拟(二)

SIR 模型(Susceptible-Infectious-Recovered Model)

有的时候,感染者在康复了之后,就有了抗体,于是后续就不再会获得此类病症,这种时候,考虑 SIS 模型就不合适了,需要考虑 SIR 模型。此类病症有麻疹,腮腺炎,风疹等。

SIR_model_1
SIR model

其微分方程是:\begin{cases} \frac{dS}{dt}=-r\beta S \frac{I}{N} \\ \frac{dI}{dt}=r\beta S\frac{I}{N} - \gamma I \\ \frac{dR}{dt}=\gamma I\end{cases}。其初始条件是 S(0)=S_{0}, I(0)=I_{0}, R(0)=R_{0},并且 S(t), I(t), R(t)\geq 0S(t) +I(t)+R(t)=N 对于所有的 t\geq 0 都成立。

对于这类方程,就不能够得到其解析解了,只能够从它的动力系统开始进行分析,得到解的信息。根据第一个微分方程可以得到:\frac{dS}{dt} = -r\beta S\frac{I}{N}<0,于是 S(t) 是一个严格递减函数。同时,0\leq S(t)\leq N 对于所有的 t\geq 0 都成立,于是存在 S_{\infty}\in[0,\infty] 使得 \lim_{t\rightarrow \infty}S(t)=S_{\infty}.

通过第一个微分方程和第二个微分方程可以得到:\frac{d(S+I)}{dt} = -\gamma I,因此对它两边积分得到 \int_{0}^{T} \frac{d(S+I)}{dt} = -\gamma \int_{0}^{T}I(t)dt. 左侧等于 S(T) + I(T) - S(0) - I(0),上界是 4N,因此令 T\rightarrow \infty 可以得到 \int_{0}^{\infty}I(t) dt\leq 4N/\gamma. 而 I(t)\geq 0 且是连续可微函数,因此 \lim_{t\rightarrow \infty}I(t) = 0。这意味着所有的感染人群都将康复。

由于 S(t) 是严格单调递减函数,因此从第二个微分方程可以得到:当 S(t) = N\gamma/(r\beta) 时,感染人数 I(t) 达到最大值。

SIR_model_graph_1
SIR model 的数值模拟(一)
SIR_model_graph_2
SIR model 的数值模拟(二)

其余模型

在以上的 SI,SIS,SIR 模型中,还可以把死亡因素考虑进去。除此之外,还有 SIRS 模型,SEIR 模型等,在这里就不再做赘述。有兴趣的读者可以参阅相关的参考书籍。

参考文献

  1. Introduction to SEIR Models, Nakul Chitnis, Workshop on Mathematical Models of Climate Variability, Environmental Change and Infectious Diseases, Trieste, Italy, 2017

 

用 Python 来研究数学 — SymPy 符号工具包介绍

SymPy 的简单介绍

SymPy 是一个符号计算的 Python 库,完全由 Python 写成,为许多数值分析,符号计算提供了重要的工具。SymPy 的第一个版本于 2007 年开源,并且经历了十几个版本的迭代,在 2019 年已经基于修正的 BSD 许可证开源了 1.4 版本。SymPy 的开源地址和官方网站分别是:

  1. GitHub 链接:https://github.com/sympy/sympy
  2. SymPy 官方网站:https://www.sympy.org/en/index.html
sympy_logo
SymPy 的 logo

SymPy 的 1.4 版本文档中,可以看出,SymPy 可以支持很多初等数学,高等数学,甚至研究生数学的符号计算。在初等数学和高等数学中,SymPy 可以支持的内容包括但不限于:

  1. 基础计算(Basic Operations);
  2. 公式简化(Simplification);
  3. 微积分(Calculus);
  4. 解方程(Solver);
  5. 矩阵(Matrices);
  6. 几何(geometry);
  7. 级数(Series);

在更多的数学领域中,SymPy 可以支持的内容包括但不限于:

  1. 范畴论(Category Theory);
  2. 微分几何(Differential Geometry);
  3. 常微分方程(ODE);
  4. 偏微分方程(PDE);
  5. 傅立叶变换(Fourier Transform);
  6. 集合论(Set Theory);
  7. 逻辑计算(Logic Theory)。
sympy_tutorial
SymPy 的教学目录

SymPy 的工具库介绍

SymPy 的基础计算

在数学中,基础的计算包括实数和复数的加减乘除,那么就需要在程序中描述出实数与复数。著名的欧拉公式

e^{i\pi}+1 = 0

正好用到了数学中最常见的五个实数。在 SymPy 里面,e, i, \pi, \infty 是用以下符号来表示的:其中 sympy.exp() 表示以 e 为底的函数。

sympy.exp(1), sympy.I, sympy.pi, sympy.oo

而想要计算欧拉公式的话,只需要输入下面的公式即可:

>>> sympy.exp(sympy.I * sympy.pi) + 1
0

如果需要看 e, \pi 的小数值,可以使用 evalf() 函数,其中 evalf() 函数里面的值表示有效数字的位数。例如下面就是精确到 10 位有效数字。当然,也可以不输入。

>>> sympy.E.evalf(10)
2.718281828
>>> sympy.E.evalf()
2.71828182845905
>>> sympy.pi.evalf(10)
3.141592654
>>> sympy.pi.evalf()
3.14159265358979

除此之外,如果需要查看某个实数的有效数字,也是类似操作的:

>>> expr = sympy.sqrt(8)
>>> expr.evalf()
2.82842712474619

而对于实数的加减乘除,则可以如下操作:

>>> x, y= sympy.symbols("x y")
>>> x + y
x + y
>>> x - y
x - y
>>> x * y
x*y
>>> x / y
x/y

而对于复数的加减乘除,则是类似的操作,令两个复数分别是 z_{1} = x_{1} + i y_{1}z_{2} = x_{2} + i y_{2}

>>> x1, y1, x2, y2 = sympy.symbols("x1 y1 x2 y2")
>>> z1 = x1 + y1 * sympy.I
x1 + I*y1
>>>  z2 = x2 + y2 * sympy.I
x2 + I*y2
>>> z1 + z2
x1 + x2 + I*y1 + I*y2
>>> z1 - z2
x1 - x2 + I*y1 - I*y2
>>> z1 * z2
(x1 + I*y1)*(x2 + I*y2)
>>> z1 / z2
(x1 + I*y1)/(x2 + I*y2)

对于多项式而言,有的时候我们希望将其展开,有的时候则需要将其合并,最终将其简化成最简单的形式。

>>> sympy.expand((x+1)**2)
x**2 + 2*x + 1
>>> sympy.expand((x+1)**5)
x**5 + 5*x**4 + 10*x**3 + 10*x**2 + 5*x + 1
>>> sympy.factor(x**3+1)
(x + 1)*(x**2 - x + 1)
>>> sympy.factor(x**2+3*x+2)
(x + 1)*(x + 2)
>>> sympy.simplify(x**2 + x + 1 - x)
x**2 + 1
>>> sympy.simplify(sympy.sin(x)**2 + sympy.cos(x)**2)
1

在多变量的场景下,SymPy 也可以对其中的某个变量合并同类项,同时还可以计算某个变量的某个次数所对应的系数是多少,例如:

>>> expr = x*y + x - 3 + 2*x**2 - x**2 + x**3 + y**2 + x**2*y**2
>>> sympy.collect(expr,x)
x**3 + x**2*(y**2 + 1) + x*(y + 1) + y**2 - 3
>>> sympy.collect(expr,y)
x**3 + x**2 + x*y + x + y**2*(x**2 + 1) - 3
>>> expr.coeff(x, 2)
y**2 + 1
>>> expr.coeff(y, 1)
x

有理函数形如 f(x) = p(x)/q(x),其中 p(x)q(x) 都是多项式。一般情况下,我们希望对有理函数进行简化,合并或者分解的数学计算。

在需要合并的情形下,如果想把有理函数处理成标准格式 p(x)/q(x) 并且去除公因子,那么可以使用 cancel 函数。另一个类似的就是 together 函数,但是不同之处在于 cancel 会消除公因子,together 不会消除公因子。例如:

expr = \frac{x^{2}+3x+2}{x^{2}+x}

>>> expr = (x**2 + 3*x + 2)/(x**2 + x)
>>> sympy.cancel(expr)
(x + 2)/x
>>> sympy.together(expr)
(x**2 + 3*x + 2)/(x*(x + 1))

除了合并和消除公因子之外,有的时候还希望对分子和分母进行因式分解,例如:

expr = (x**2 + 3*x + 2)/(x**2 + x)
>>> sympy.factor(expr)
(x + 2)/x
>>> expr = (x**3 + 3*x**2 + 2*x)/(x**5+x)
>>> sympy.factor(expr)
(x + 1)*(x + 2)/(x**4 + 1)
>>> expr = x**2 + (2*x+1)/(x**3+1)
>>> sympy.factor(expr)
(x**5 + x**2 + 2*x + 1)/((x + 1)*(x**2 - x + 1))

合并的反面就是部分分式展开(Partial Fraction Decomposition),它是把有理函数分解成多个次数较低的有理函数和的形式。这里需要用 apart 函数:

>>> expr = (x**4 + 3*x**2 + 2*x)/(x**2+x)
>>> sympy.apart(expr)
x**2 - x + 4 - 2/(x + 1)
>>> expr = (x**5 + 1)/(x**3+1)
>>> sympy.apart(expr)
x**2 - (x - 1)/(x**2 - x + 1)

在 SymPy 里面,同样支持各种各样的三角函数,例如:三角函数的简化函数 trigsimp,三角函数的展开 expand_trig,

>>> expr = sympy.sin(x)**2 + sympy.cos(x)**2
>>> sympy.trigsimp(expr)
1
>>> sympy.expand_trig(sympy.sin(x+y))
sin(x)*cos(y) + sin(y)*cos(x)
>>> sympy.expand_trig(sympy.cos(x+y))
-sin(x)*sin(y) + cos(x)*cos(y)
>>> sympy.trigsimp(sympy.sin(x)*sympy.cos(y) + sympy.sin(y)*sympy.cos(x))
sin(x + y)
>>> sympy.trigsimp(-sympy.sin(x)*sympy.sin(y) + sympy.cos(x)*sympy.cos(y))
cos(x + y)

同样的,在乘幂上面,同样有简化函数 powsimp,效果与之前提到的 simplify 一样。除此之外,还可以根据底数来做合并,即分别使用 expand_power_exp 函数与 expand_power_base 函数。

>>> sympy.powsimp(x**z*y**z*x**z)
x**(2*z)*y**z
>>> sympy.simplify(x**z*y**z*x**z)
x**(2*z)*y**z
>>> sympy.expand_power_exp(x**(y + z))
x**y*x**z
>>> sympy.expand_power_base(x**(y + z))
x**(y + z)

作为指数的反函数对数,sympy 也是有着类似的展开合并函数,expand_log,logcombine 承担了这样的角色。

\ln(xy) = \ln(x) + \ln(y)

\ln(x/y) = \ln(x) - \ln(y)

>>> sympy.expand_log(sympy.log(x*y), force=True)
log(x) + log(y)
>>> sympy.expand_log(sympy.log(x/y), force=True)
log(x) - log(y)

 

SymPy 的微积分工具

下面,我们会从一个最基本的函数 f(x) = 1/x 出发,来介绍 SymPy 的各种函数使用方法。如果想进行变量替换,例如把 x 变成 y,那么可以使用 substitution 方法。除此之外,有的时候也希望能够得到函数 f 在某个点的取值,例如 f(1) = 1/1 = 1,那么可以把参数换成 1 即可得到函数的取值。例如,

>>> import sympy
>>> x = sympy.Symbol("x")
>>> f = 1 / x
1/x
>>> y = sympy.Symbol("y")
>>> f = f.subs(x,y)
1/y
>>> f = f.subs(y,1)
1

在微积分里面,最常见的概念就是极限,SymPy 里面的极限函数就是 limit。使用方法如下:

>>> f = 1/x
>>> sympy.limit(f,x,0)
oo
>>> sympy.limit(f,x,2)
1/2
>>> sympy.limit(f,x,sympy.oo)
0
>>> g = x * sympy.log(x)
>>> sympy.limit(g,x,0)
0

对于函数 f(x) = 1/x 而言,它的导数计算函数是 diff,n 阶导数也可以用这个函数算。

>>> f = 1/x
>>> sympy.diff(f,x)
-1/x**2
>>> sympy.diff(f,x,2)
2/x**3
>>> sympy.diff(f,x,3)
-6/x**4
>>> sympy.diff(f,x,4)
24/x**5

提到 n 阶导数,就必须要提一下 Taylor Series 了。对于常见函数的 Taylor Series,SymPy 也是有非常简便的方法,那就是 series 函数。其参数包括 expr, x, x0, n, dir,分别对应着表达式,函数的自变量,Taylor Series 的中心点,n 表示阶数,dir 表示方向,包括”+-“,”-“,”+”,分别表示 x\rightarrow x0, x\rightarrow x0^{-}, x\rightarrow x0^{+}

sympy.series.series.series(exprx=Nonex0=0n=6dir='+') >>> g = sympy.cos(x) >>> sympy.series(g, x) 1 - x**2/2 + x**4/24 + O(x**6) >>> sympy.series(g, x, x0=1, n=10) cos(1) - (x - 1)*sin(1) - (x - 1)**2*cos(1)/2 + (x - 1)**3*sin(1)/6 + (x - 1)**4*cos(1)/24 - (x - 1)**5*sin(1)/120 - (x - 1)**6*cos(1)/720 + (x - 1)**7*sin(1)/5040 + (x - 1)**8*cos(1)/40320 - (x - 1)**9*sin(1)/362880 + O((x - 1)**10, (x, 1))

积分的计算函数是 integrate,包括定积分与不定积分:

\int\frac{1}{x}dx = \ln(x)+C

\int_{1}^{2}\frac{1}{x}dx = \ln(2)

>>> f = 1/x
>>> sympy.integrate(f,x)
log(x)
>>> sympy.integrate(f, (x,1,2))
log(2)

对于广义积分而言,就需要用到 \infty 这个概念了,但是在 SymPy 里面的写法还是一样的。

\int_{-\infty}^{0}e^{-x^{2}}dx=\frac{\sqrt{\pi}}{2}

\int_{0}^{+\infty}e^{-x}dx = 1

\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}e^{-x^{2}-y^{2}}dxdy = \pi

>>> g = sympy.exp(-x**2)
>>> sympy.integrate(g, (x,-sympy.oo,0))
sqrt(pi)/2
>>> g = sympy.exp(-x)
>>> sympy.integrate(g, (x, 0, sympy.oo))
1
>>> h = sympy.exp(-x**2 - y**2)
>>> sympy.integrate(h, (x,-sympy.oo, sympy.oo), (y, -sympy.oo, sympy.oo))
pi

 

SymPy 的方程工具

在初等数学中,通常都存在一元一次方程,一元二次方程等,并且在不同的域上有着不同的解。SymPy 里面的相应函数就是 solveset,根据定义域的不同,可以获得完全不同的解。

\{x\in\mathbb{R}: x^{3}-1=0\}

\{x\in\mathbb{C}:x^{3}-1=0\}

\{x\in\mathbb{R}:e^{x}-x=0\}

\{x\in\mathbb{R}:e^{x}-1=0\}

\{x\in\mathbb{C}:e^{x}-1=0\}

>>> sympy.solveset(sympy.Eq(x**3,1), x, domain=sympy.S.Reals)
{1}
>>> sympy.solveset(sympy.Eq(x**3,1), x, domain=sympy.S.Complexes)
{1, -1/2 - sqrt(3)*I/2, -1/2 + sqrt(3)*I/2}
>>> sympy.solveset(sympy.Eq(x**3 - 1,0), x, domain=sympy.S.Reals)
{1}
>>> sympy.solveset(sympy.Eq(x**3 - 1,0), x, domain=sympy.S.Complexes)
{1, -1/2 - sqrt(3)*I/2, -1/2 + sqrt(3)*I/2}
>>> sympy.solveset(sympy.exp(x),x)
EmptySet()
>>> sympy.solveset(sympy.exp(x)-1,x,domain=sympy.S.Reals)
{0}
>>> sympy.solveset(sympy.exp(x)-1,x,domain=sympy.S.Complexes)
ImageSet(Lambda(_n, 2*_n*I*pi), Integers)

在这里,Lambda 表示的是数学公式,第一个是自变量,第二个是函数,最后是自变量的定义域。

在线性代数中,最常见的还是多元一次方程组,那么解法是一样的:

\begin{cases}x+y-10=0 \\ x-y-2=0\end{cases}

>>> sympy.solve([x+y-10, x-y-2], [x,y])
{x: 6, y: 4}

对于三角函数,也是类似的写法:

\begin{cases} \sin(x-y)=0 \\ \cos(x+y)=0 \end{cases}

>>> sympy.solve([sympy.sin(x-y), sympy.cos(x+y)], [x,y])
[(-pi/4, 3*pi/4), (pi/4, pi/4), (3*pi/4, 3*pi/4), (5*pi/4, pi/4)]

 

SymPy 的矩阵工具

在矩阵论中,最常见的就是单位矩阵了,而单位矩阵只与一个参数有关,那就是矩阵的大小。下面就是 3*3,3*2,2*3 大小的矩阵。

>>> sympy.eye(3)
Matrix([
[1, 0, 0],
[0, 1, 0],
[0, 0, 1]])
>>> sympy.eye(3,2)
Matrix([
[1, 0],
[0, 1],
[0, 0]])
>>> sympy.eye(2,3)
Matrix([
[1, 0, 0],
[0, 1, 0]])

另外还有全零和全一矩阵,意思就是矩阵中的所有值全部是零和一。

>>> sympy.ones(2,3)
Matrix([
[1, 1, 1],
[1, 1, 1]])
>>> sympy.zeros(3,2)
Matrix([
[0, 0],
[0, 0],
[0, 0]])

而对角矩阵也可以使用 diag 轻松获得:

>>> sympy.diag(1,1,2)
Matrix([
[1, 0, 0],
[0, 1, 0],
[0, 0, 2]])

而矩阵的加法,减法,乘法,逆运算,转置,行列式,SymPy 都是可以支持的:

A = \left(\begin{array}{cc} 1 & 1 \\ 0 & 2 \end{array}\right)

B = \left(\begin{array}{cc} 1 & 0 \\ 1 & 1 \end{array}\right)

>>> A = sympy.Matrix([[1,1],[0,2]])
>>> B = sympy.Matrix([[1,0],[1,1]])
>>> A
Matrix([
[1, 1],
[0, 2]])
>>> B
Matrix([
[1, 0],
[1, 1]])
>>> A + B
Matrix([
[2, 1],
[1, 3]])
>>> A - B
Matrix([
[ 0, 1],
[-1, 1]])
>>> A * B
Matrix([
[2, 1],
[2, 2]])
>>> A * B**-1
Matrix([
[ 0, 1],
[-2, 2]])
>>> B**-1
Matrix([
[ 1, 0],
[-1, 1]])
>>> A.T
Matrix([
[1, 0],
[1, 2]])
>>> A.det()
2

在某些情况下,需要对矩阵进行加上一行或者加上一列的操作,在这里就可以使用 row_insert 或者 col_insert 函数:第一个参数表示插入的位置,第二个参数就是相应的行向量或者列向量。而在删除上就很简单了,直接使用 row_del 或者 col_del 即可。

>>> A
Matrix([
[1, 1],
[0, 2]])
>>> A.row_insert(1, sympy.Matrix([[10,10]]))
Matrix([
[ 1, 1],
[10, 10],
[ 0, 2]])
>>> A.col_insert(0, sympy.Matrix([3,3]))
Matrix([
[3, 1, 1],
[3, 0, 2]])
>>> A.row_del(0)
>>> A
Matrix([[0, 2]])
>>> A.col_del(1)
>>> A
Matrix([[0]])

在对角化方面,同样可以使用 eigenvals(),eigenvecs(), diagonalize() 函数:

>>> A
Matrix([
[1, 1],
[0, 2]])
>>> A.eigenvals()
{2: 1, 1: 1}
>>> A.eigenvects()
[(1, 1, [Matrix([
[1],
[0]])]), (2, 1, [Matrix([
[1],
[1]])])]
>>> A.diagonalize()
(Matrix([
[1, 1],
[0, 1]]), Matrix([
[1, 0],
[0, 2]]))

在 eigenvals() 返回的结果中,第一个表示特征值,第二个表示该特征值的重数。在特征向量 eigenvecs() 中,第一个表示特征值,第二个表示特征值的重数,第三个表示特征向量。在对角化 diagonalize() 中,第一个矩阵表示 P,第二个矩阵表示 DA = P*D*P^{-1}

在矩阵中,最常见的还是多元一次方程的解。如果要求 Ax =b 的解,可以有以下方案:

>>> A = sympy.Matrix([[1,1],[0,2]])
>>> A
Matrix([
[1, 1],
[0, 2]])
>>> b = sympy.Matrix([3,5])
>>> b
Matrix([
[3],
[5]])
>>> A**-1*b
Matrix([
[1/2],
[5/2]])
>>> sympy.linsolve((A,b))
{(1/2, 5/2)}
>>> sympy.linsolve((A,b),[x,y])
{(1/2, 5/2)}

 

SymPy 的集合论工具

集合论可以说是数学的基础,在任何数学的方向上都能够看到集合论的身影。在 SymPy 里面,有一个类叫做 sympy.sets.sets.set。在集合论里面,常见的就是边界,补集,包含,并集,交集等常见的操作。但是感觉 SymPy 中的集合论操作主要集中在实数域或者复数域上。

对于闭区间 I=[0,1] 和开区间 J = (0,1) 而言,在 SymPy 中使用以下方法来表示:

I = sympy.Interval(0,1)
J = sympy.Interval.open(0,1)
K = sympy.Interval(0.5,2)

其开始和结束的点可以分别使用 start 和 end 来表示:

>>> I.start
0
>>> I.end
1

其长度用 measure 来表示:

>>> I.measure
1

其闭包用 closure 来表示:

>>> I.closure
Interval(0, 1)

其内点用 interior 来表示:

>>> I.interior
Interval.open(0, 1)

判断其边界条件可以使用 left_open 或者 right_open 来做:

>>> I.left_open
False
>>> I.right_open
False

对于两个集合之间的操作,可以参考以下方法:

I = sympy.Interval(0,1)
K = sympy.Interval(0.5,2)
>>> I.intersect(K)
Interval(0.500000000000000, 1)
>>> I.union(K)
Interval(0, 2)
>>> I-K
Interval.Ropen(0, 0.500000000000000)
>>> K-I
Interval.Lopen(1, 2)
>>> I.symmetric_difference(K)
Union(Interval.Ropen(0, 0.500000000000000), Interval.Lopen(1, 2))

实数集 \mathbb{R} 在 SymPy 中用 sympy.S.Reals 来表示,自然数使用 sympy.S.Naturals,非负整数用 sympy.S.Naturals0,整数用 sympy.S.Integers 来表示。补集的计算可以用减号,也可以使用 complement 函数。

>>> sympy.S.Reals
Reals
>>> sympy.S.Reals-I
Union(Interval.open(-oo, 0), Interval.open(1, oo))
>>> I.complement(sympy.S.Reals)
Union(Interval.open(-oo, 0), Interval.open(1, oo))
>>> sympy.S.Reals.complement(I)
EmptySet()
>>> I.complement(K)
Interval.Lopen(1, 2)
>>> I.complement(sympy.S.Reals)
Union(Interval.open(-oo, 0), Interval.open(1, oo))

 

SymPy 的逻辑工具

在逻辑运算中,我们可以使用 A, B, C 来代表元素。&, |, ~, >> 分别表示 AND,OR,NOT,imply。而逻辑运算同样可以使用 sympy.simplify_logic 简化。

A, B, C = sympy.symbols("A B C")
>>> sympy.simplify_logic(A | (A & B))
A
>>> sympy.simplify_logic((A>>B) & (B>>A))
(A & B) | (~A & ~B)
>>> A>>B
Implies(A, B)

 

SymPy 的级数工具

SymPy 的级数工具有一部分放在具体数学(Concrete Mathematics)章节了。有的时候,我们希望计算某个级数是发散的,还是收敛的,就可以使用 is_convergence 函数。考虑最常见的级数:

\sum_{n=1}^{\infty}\frac{1}{n} = +\infty

\sum_{n=1}^{\infty}\frac{1}{n^{2}} = \frac{\pi^{2}}{6}

>>> n = sympy.Symbol("n", integer=True)
>>> sympy.Sum(1/n, (n,1,sympy.oo)).is_convergent()
False
>>> sympy.Sum(1/n**2, (n,1,sympy.oo)).is_convergent()
True

如果想计算出收敛级数的值,加上 doit() 函数即可;如果想计算有效数字,加上 evalf() 函数即可。

>>> sympy.Sum(1/n**2, (n,1,sympy.oo)).evalf()
1.64493406684823
>>> sympy.Sum(1/n**2, (n,1,sympy.oo)).doit()
pi**2/6
>>> sympy.Sum(1/n**3, (n,1,sympy.oo)).evalf()
1.20205690315959
>>> sympy.Sum(1/n**3, (n,1,sympy.oo)).doit()
zeta(3)

除了加法之外,SymPy 也支持连乘,其符号是 sympy.Product,考虑

\prod_{n=1}^{+\infty}\frac{n}{n+1}

\prod_{n=1}^{+\infty}\cos\left(\frac{\pi}{n}\right)

>>> sympy.Product(n/(n+1), (n,1,sympy.oo)).is_convergent()
False
>>> sympy.Product(sympy.cos(sympy.pi/n), (n, 1, sympy.oo)).is_convergent()
True

 

SymPy 的 ODE 工具

在常微分方程(Ordinary Differential Equation)中,最常见的就是解方程,而解方程主要是靠 dsolve 函数。例如想求解以下的常微分方程:

df/dx + f(x) = 0,

d^{2}f/dx^{2} + f(x) = 0

d^{3}f/dx^{3} + f(x) = 0

可以使用 dsolve 函数:

>>> f = sympy.Function('f')
>>> sympy.dsolve(sympy.Derivative(f(x),x) + f(x), f(x))
Eq(f(x), C1*exp(-x))
>>> sympy.dsolve(sympy.Derivative(f(x),x,2) + f(x), f(x))
Eq(f(x), C1*sin(x) + C2*cos(x))
>>> sympy.dsolve(sympy.Derivative(f(x),x,3) + f(x), f(x))
Eq(f(x), C3*exp(-x) + (C1*sin(sqrt(3)*x/2) + C2*cos(sqrt(3)*x/2))*sqrt(exp(x)))

而常微分方程对于不同的方程类型也有着不同的解法,可以使用 classify_ode 来判断常微分方程的类型:

>>> sympy.classify_ode(sympy.Derivative(f(x),x) + f(x), f(x))
('separable', '1st_exact', '1st_linear', 'almost_linear', '1st_power_series', 'lie_group', 'nth_linear_constant_coeff_homogeneous', 'separable_Integral', '1st_exact_Integral', '1st_linear_Integral', 'almost_linear_Integral')
>>> sympy.classify_ode(sympy.Derivative(f(x),x,2) + f(x), f(x))
('nth_linear_constant_coeff_homogeneous', '2nd_power_series_ordinary')
>>> sympy.classify_ode(sympy.Derivative(f(x),x,3) + f(x), f(x))
('nth_linear_constant_coeff_homogeneous',)

 

SymPy 的 PDE 工具

在偏微分方程(Partitial Differential Equation)中,同样可以直接求解和判断偏微分方程的类型,分别使用函数 pdsolve() 和 classify_pde()。假设 f = f(x,y) 是一个二元函数,分别满足以下偏微分方程:

\partial f/\partial x + \partial f/\partial y =0

\partial f/\partial x + \partial f/\partial y + f = 0

\partial f/\partial x + \partial f/\partial y + f + 10 = 0

>>> f = sympy.Function("f")(x,y)
>>> sympy.pdsolve(sympy.Derivative(f,x)+sympy.Derivative(f,y),f)
Eq(f(x, y), F(x - y))
>>> sympy.pdsolve(f.diff(x)+f.diff(y)+f,f)
Eq(f(x, y), F(x - y)*exp(-x/2 - y/2))
>>> sympy.pdsolve(f.diff(x)+f.diff(y)+f+10,f)
Eq(f(x, y), F(x - y)*exp(-x/2 - y/2) - 10)

查看类型就用 classify_pde() 函数:

>>> sympy.classify_pde(f.diff(x)+f.diff(y)+f)
('1st_linear_constant_coeff_homogeneous',)
>>> sympy.classify_pde(f.diff(x)+f.diff(y)+f+10,f)
('1st_linear_constant_coeff', '1st_linear_constant_coeff_Integral')
>>> sympy.classify_pde(f.diff(x)+f.diff(y)+f+10,f)
('1st_linear_constant_coeff', '1st_linear_constant_coeff_Integral')

不过目前的 PDE 解法貌似只支持一阶偏导数,二阶或者以上的偏导数就不支持了。

 

SymPy 的数论工具

在数论中,素数就是一个最基本的概念之一。而素数的批量计算,比较快的方法就是筛法(sieve method)。在 sympy 中,同样有 sympy.sieve 这个工具,用于计算素数。如果想输出前100个素数,那么

>>> sympy.sieve._reset()
>>> sympy.sieve.extend_to_no(100)
>>> sympy.sieve._list
array('l', [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 631])

如果想输出一个区间内的所有素数,可以使用 primerange(a,b) 函数:

>>> [i for i in sympy.sieve.primerange(10,100)]
[11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]

search() 函数是为了计算某个数附近是第几个素数:

>>> sympy.sieve.search(10)
(4, 5)
>>> sympy.sieve.search(11)
(5, 5)

如果只想获得第 n 个素数,则使用函数 sympy.ntheory.generate.prime(n) 即可。如果是希望计算 x 后面的下一个素数,使用 sympy.ntheory.generate.nextprime(x) 即可。判断 x 是否是素数,可以使用 sympy.ntheory.generate.isprime(x)。

>>> sympy.ntheory.generate.prime(10)
29
>>> sympy.ntheory.generate.nextprime(10)
11
>>> sympy.ntheory.generate.nextprime(11)
13
>>> sympy.ntheory.generate.isprime(11)
True
>>> sympy.ntheory.generate.isprime(12)
False

除此之外,SymPy 的数论方法还有很多,需要读者根据 SymPy 的官方文档自行探索。

 

SymPy 的范畴论工具

SymPy 还支持范畴论(Category Theory)的一些计算方法,在这里简要地列举一下。

>>> A = sympy.categories.Object("A")
>>> B = sympy.categories.Object("B")
>>> f = sympy.categories.NamedMorphism(A,B,"f")
>>> f.domain
Object("A")
>>> f.codomain
Object("B")

由于范畴论是数学的“黑话”,因此其余方法留给范畴论的科研人员自行翻阅。

总结:

整体来看,SymPy 是一个非常卓越的 Python 开源符号计算库。在符号计算领域,不仅支持常见的微积分,线性代数,几何运算,还支持集合论,微分方程,数论等诸多数学方向。后续笔者将会持续跟进并研究这一卓越的开源工具库。

 

参考文献:

  1. Meurer A, Smith C P, Paprocki M, et al. SymPy: symbolic computing in Python[J]. PeerJ Computer Science, 2017, 3: e103.
  2. GitHub:https://github.com/sympy/sympy
  3. SymPy:https://www.sympy.org/en/index.html
  4. Sympy 维基百科:https://en.wikipedia.org/wiki/SymPy
  5. GreatX’s Blog:数值 Python:符号计算:https://vlight.me/2018/04/01/Numerical-Python-Symbolic-Computing/
  6. SymPy 符号计算-让Python帮我们推公式:https://zhuanlan.zhihu.com/p/83822118
  7. Python 科学计算利器—SymPy库:https://www.jianshu.com/p/339c91ae9f41

 

本科学数学专业是一个很好的选择吗?

知乎问题:https://www.zhihu.com/question/319574112

选专业这件事情其实是因人而异的,很难对所有的人给出一个标准的答案,肯定是基于每个人的不同条件来给出完全不同的建议。这是一个千人千面的问题,而不是一个数学题,通常只有一个标准的解答。

就个人这几年的经验来看,无论选择什么专业,都会有这个专业的优势和劣势,好比科技是一把双刃剑,专业也是同样的道理。在这种情况下,我们就需要分析数学专业究竟有哪些优势,有哪些劣势。只有这样,我们才能够根据自身的情况来具体分析,然后决定最终是否选择数学专业。

数学专业的劣势

于是,我们来看一下数学专业的劣势有哪些:

  • 理论知识太多;
  • 实用技能偏少;
  • 转行需要时间。

下面来逐一解读以上几点。我们可以先看一下数学系的常见课表:

  • 第一年:数学分析,高等代数,解析几何,C++等;
  • 第二年:常微分方程,离散数学,复分析,概率论,数值计算,抽象代数等;
  • 第三年:实分析,泛函分析,偏微分方程,拓扑学,微分几何,偏微分方程数值解,随机过程,数理统计等。

从课表上面来看,基本上可以确定几个结论。首先,数学专业作为基础学科,其特点就是理论知识偏多,而学习到的技能偏少,毕竟所学的内容都是理论型,培养的学生都是理论型选手。因此直接导致的结果就是数学系的学生掌握了一堆理论,但是却没有办法把它们直接转化成生产力。在实战中,总不能就靠一门 C++ 来谋求工作吧。其次,既然数学系传授给学生的实用的技能偏少,那么数学系的学生在需要转行的话,就肯定要补充新的技能。在从理论派走向实战派的过程中,不仅要找好自己的前进方向,还需要花费一定的时间和精力去转行。在这里需要澄清一点,转行并不是轻轻松松,而是需要花费时间,勇气和精力的。

如果不想继续从事数学科研的话,其实还是建议数学系的学生可以早一点进入公司去实习或者工作,至少在公司能够体验一下与学术界完全不同的人生。人生总是有多种可能性的,其实可以在本科或者硕士阶段多去体验一下人生。与数学系不同的是,对于计算机或者工程类专业的学生而言,到了本科一定的阶段,都会从事某个项目或者大作业,这种时候他们就会在边看边学中得到一定的成长,实践能力的训练其实比数学系的人会早很多。其实数学系也有实践,只不过延后了许多,一般只有到了博士生的阶段才需要进行科研的训练和动手的操作,在本科和硕士阶段一般是不需要的,因为现代数学的发展已经不是大部分硕士生能够完成的了,当然优秀的人总是有的。

数学专业的优势

在讲了数学专业的劣势之后,也需要强调一下数学专业的优势,其优势包括:

  • 底层通用技能;
  • 技能不易淘汰;
  • 逻辑思维能力;
  • 转行就业面广。

众所周知,无论是在学术界,还是工业界,数学基本上就是一切的基础。如果计算机行业没有数学,那就是XX计算机学院与XX培训班的区别;如果金融行业没有数学,那就是文艺复兴公司与XX小银行的区别。虽然我们不能够一味的拔高数学在各个行业的作用,但是很多行业还是离不开数学的。这就是所谓的底层通用技能,无论是计算机,金融还是其他领域,都离不开数学的支持。

除此之外,再次回到那张数学系的常见课表:

  • 第一年:数学分析,高等代数,解析几何,C++等;
  • 第二年:常微分方程,离散数学,复分析,概率论,数值计算,抽象代数等;
  • 第三年:实分析,泛函分析,偏微分方程,拓扑学,微分几何,偏微分方程数值解,随机过程,数理统计等。

当年读本科的时候是这张课表,其实过了十年,也是这张课表。本科的数学课程基本上集中在20世纪初之前的数学内容,最多到了20世纪中期。而微积分的发展时间则更加久远了。对于数学系的教育而言,很难做到跳过数学分析,高等代数的教育直接进入实分析和泛函分析。就算老师能够教,学生也听不懂啊,还是只能够从基础一步一步开始。而工业界用到的数学,通常也就是数学本科三年级的所有课程就能够包括了。很难用到很多研究生方面的知识,甚至很多时候也就用用微积分和线性代数,概率论就足够了。因此,一旦学会了这些课程,则是终身受益的知识,因为数学的另外一个特点就是永恒性。无论个人发生什么,学校发生什么,世界发生了什么,数学定理就是数学定理,一旦被证明且确定了证明是正确的,那就是永恒的。个人会死亡,学校也有可能走向没落,世界也有可能发生变化,但是数学定理就像一个永恒的石头永远放在那里。

除了以上两点,数学是最能够培养学生逻辑思维的学科。在以上的本科生课程里面,几乎所有东西都是从几个公理出发,然后通过严格的证明,一点一点地得到最终的结论,并且构建出整个数学大厦。在本科教育里面,数学系的学生只要认真学习,通常来说,逻辑思维能力和数学推导能力都会得到一个很大的提升。并且在后续的学习或者工作中,数学的烙印都会深深地印在身上。

最后,数学的就业面其实是相对宽广很多,主要包括:

  • 科研工作者:数学界,金融界,经济界,计算机方向等;
  • 计算机行业;
  • 金融行业;
  • 教育培训行业;
  • 其他行业。

除了可以继续从事本专业之外,其他方向无论是金融还是计算机都可以转。

结论: 整体来看,其实如果自身条件 OK 的话,并且也愿意在本科期间选择数学专业的话。其实选择数学专业是一个不错的选择。在数学系本科这几年可以根据自身的情况来继续选择合适自己的发展方向,并且在研究生或者工作的时候选择一个适合自己的舞台。

Riemann Zeta 函数(二)

在上一篇文章里面,我们已经给出了 Riemann Zeta 函数的定义,

\zeta(s) = \sum_{n=1}^{\infty} \frac{1}{n^s}.

其定义域是 [1,\infty)\subseteq\mathbb{R}. 根据级数与定积分的等价关系可以得到:

  1. s = 1 时,\zeta(1) = \infty;
  2. s>1 时,\zeta(s)<\infty.

本文将会重点讲两个内容:

  1. 如何把 Riemann Zeta 函数从 [1,\infty)\subseteq \mathbb{R} 上延拓到 \{s\in \mathbb{C}: \Re(s)>0\} 上;
  2. Riemann Zeta 函数在 \{s\in\mathbb{C}: \Re(s)\geq 1\} 上没有零点。

Riemann Zeta 函数定义域的延拓

如果想把 Riemann Zeta 函数的定义域从 [1,\infty)\subseteq \mathbb{R} 延拓到更大的区域 \{s\in\mathbb{C}:\Re(s)>0\} 上,就需要给出 Riemann Zeta 函数在 \{s\in \mathbb{C}: \Re(s)>0\} 上的定义。而且在原始的定义域 [1,\infty)\subseteq\mathbb{R} 上面,新的函数的取值必须与原函数的取值保持一致。

首先,我们将会在 [1,\infty)\subseteq \mathbb{R} 上面证明如下恒等式:

\zeta(s) = \frac{s}{s-1} - s\int_{1}^{\infty}\frac{\{x\}}{x^{s+1}}dx.

证明:当 s=1 时,上述等式显然成立,两侧都是 \infty.

\frac{s}{s-1}-s\int_{1}^{\infty}\frac{\{x\}}{x^{s+1}}dx

= \frac{s}{s-1} - s\sum_{n=1}^{\infty}\int_{n}^{n+1}\frac{\{x\}}{x^{s+1}}dx 

= \frac{s}{s-1} - s\sum_{n=1}^{\infty}\int_{n}^{n+1}\frac{x-n}{x^{s+1}}dx 

= \frac{s}{s-1} - s\sum_{n=1}^{\infty}\bigg(\int_{n}^{n+1}\frac{1}{x^{s}}dx - \int_{n}^{n+1}\frac{n}{x^{s+1}}dx\bigg)

= \frac{s}{s-1} - s\int_{1}^{\infty}\frac{1}{x^{s}}dx + \sum_{n=1}^{\infty}n\cdot\int_{n}^{n+1}\frac{s}{x^{s+1}}dx

= \sum_{n=1}^{\infty}n\cdot\bigg(\frac{1}{n^{s}}-\frac{1}{(n+1)^{s}}\bigg)

= \sum_{n=1}^{\infty}\bigg(\frac{1}{n^{s-1}}-\frac{1}{(n+1)^{s-1}} + \frac{1}{(n+1)^{s}}\bigg)

= \sum_{n=1}^{\infty}\frac{1}{n^{s}}.

从右式的表达式

\frac{s}{s-1} - s \int_{1}^{\infty}\frac{\{x\}}{x^{s+1}}dx

可以看出 \zeta(s) 可以延拓到 \{s \in\mathbb{C}:\Re(s)>0\} 上。而且右侧的函数在 \{s\in\mathbb{C}:\Re(s)>0,s\neq 1\} 是解析的,并且 s=1 是该函数的一个极点。进一步的分析可以得到,我们得到一个关于 (s-1)\zeta(s) 的解析函数,而且 \lim_{s\rightarrow 1}(s-1)\zeta(s)=1. 综上所述:

  1. Riemann Zeta 函数可以延拓到 \{s\in\mathbb{C}:\Re(s)>0\} 上;
  2. Riemann Zeta 函数在 \{s\in\mathbb{C}:\Re(s)>0, s\neq 1\} 上是解析的;s=1 是 Riemann Zeta 函数的极点。

 

Riemann Zeta 函数的非零区域

著名的 Riemann 猜想说的是 \zeta(s) 函数的所有非平凡零点都在直线 \{s\in\mathbb{C}:\Re(s)=1/2\} 上。因此,数学家首先要找出的就是 Riemann Zeta 函数的非零区域。而本篇文章将会证明 Riemann Zeta 函数在 \{s\in\mathbb{C}:\Re(s)\geq 1\} 上面没有零点。

\Re(s)>1 区域

首先,我们要证明当 \Re(s)>1 时,\zeta(s)\neq 0.

在这里,就需要使用一个重要的恒等式:当 \Re(s)>1 时,

\zeta(s) =\sum_{n=1}^{\infty}\frac{1}{n^{s}}

= \prod_{p}\bigg(1+\frac{1}{p^{s}}+\frac{1}{p^{2s}}+\cdots\bigg)

= \prod_{n=1}^{\infty}\bigg(1-\frac{1}{p_{n}^{s}}\bigg)^{-1},

其中这里的 p 表示所有的素数相乘,而 p_{n} 表示第 n 个素数。

下面我们证明:

\bigg|1-\frac{1}{p_{n}^{s}}\bigg|^{-1}\geq 1-\frac{1}{p_{n}^{\sigma}-1} .

事实上,令 s = \sigma + i t,,当 \sigma=\Re(s)>1 时,我们有

\bigg|1-\frac{1}{p_{n}^{s}}\bigg|^{-1} = \bigg(1+\frac{1}{p_{n}^{s}}+\frac{1}{p_{n}^{2s}}+\cdots\bigg)

\geq 1-\frac{1}{|p_{n}^{s}|}- \frac{1}{|p_{n}^{2s}|} -\cdots

= 1- \frac{1}{p_{n}^{\sigma}} - \frac{1}{p_{n}^{2\sigma}} -\cdots

= 1- \frac{1}{p_{n}^{\sigma}-1}.

因此,

|\zeta(s)| \geq \prod_{n=1}^{\infty}\bigg|1-\frac{1}{p_{n}^{s}}\bigg|^{-1} \geq\prod_{n=1}^{\infty}\bigg(1-\frac{1}{p_{n}^{\sigma}-1}\bigg).

同时,

\lim_{n\rightarrow \infty} \bigg(1- \frac{1}{p_{n}^{\sigma}-1}\bigg) = 1 ,

1-\frac{1}{p_{n+1}^{\sigma}-1} \geq 1- \frac{1}{p_{n}^{\sigma}-1} ,

\sum_{n=1}^{\infty}\frac{1}{p_{n}^{\sigma}}\leq \sum_{n=1}^{\infty}\frac{1}{n^{\sigma}}<\infty when \sigma>1.

所以,当 \Re(s)>1 时,\zeta(s) \neq 0.

\Re(s) =1 直线

Claim 1. 下面我们将会证明恒等式:对于 \sigma >1, \text{ } t\in\mathbb{R},

\Re(\ln\zeta(\sigma + it)) = \sum_{n=2}^{\infty}\frac{\Lambda(n)}{n^{\sigma}\ln(n)}\cos(t\ln(n)) ,

其中当 n 形如 p^{\alpha}, p 是素数,\alpha \geq 1. \Lambda(n) = \ln(p). 而对于其余的 n, \Lambda(n)=0.

事实上,根据 Euler 公式,

\zeta(s) = \prod_{p}\bigg(1-\frac{1}{p^{s}}\bigg)^{-1}.

s = \sigma + it, 可以得到

\ln\zeta(s) = -\sum_{p}\ln\bigg(1-\frac{1}{p^{s}}\bigg)

= \sum_{p}\sum_{\alpha=1}^{\infty}\frac{1}{\alpha p^{\alpha s}}

= \sum_{p}\sum_{\alpha=1}^{\infty}\frac{1}{\alpha p^{\alpha\sigma}}\cdot p^{-i\alpha t}

= \sum_{p}\sum_{\alpha = 1}^{\infty}\frac{1}{\alpha p^{\alpha\sigma}}\cdot e^{-i\alpha t \ln p}

进一步,

\Re(\ln\zeta(s)) = \sum_{p}\sum_{\alpha =1}^{\infty}\frac{1}{\alpha p^{\alpha\sigma}}\cos(\alpha t \ln p)

并且右侧等于

RHS = \sum_{n=2}^{\infty}\frac{\Lambda(n)}{n^{\sigma}\ln(n)}\cos(t\ln(n))

= \sum_{p}\sum_{\alpha = 1}^{\infty} \frac{\ln(p)}{p^{\alpha\sigma}\ln(p^{\alpha})}\cos(t\ln(p^{\alpha}))

= \sum_{p}\sum_{\alpha = 1}^{\infty}\frac{1}{\alpha p^{\alpha\sigma}}\cos(\alpha t\ln p).

所以,恒等式成立,Claim 1 证明完毕。

Claim 2.

\Re(3\ln\zeta(\sigma) + 4\ln\zeta(\sigma+it) + \ln\zeta(\sigma+2it))\geq 0,

其中 \sigma>1, t\in\mathbb{R}. 换句话说

|\zeta(\sigma)^{3}\zeta(\sigma+it)^{4}\zeta(\sigma+2it)|\geq 1.

事实上,

从三角函数的性质可以得到:

3+4\cos(\theta)+\cos(2\theta) = 3 + 4\cos(\theta)+2\cos^{2}(\theta)-1

= 2(\cos(\theta)-1)^{2}\geq 0,

所以,从 Claim 1 可以得到

\Re(3\ln\zeta(\sigma) + 4\ln\zeta(\sigma+it) + \ln\zeta(\sigma+2it))

= \sum_{n=2}^{\infty} \frac{\Lambda(n)}{n^{\sigma}\ln(n)} \cdot ( 3 + 4\cos(t\ln(n)) + \cos(2t\ln(n))) \geq 0.

进一步地,使用 \Re(\ln(z)) = \ln(|z|) 可以得到

0\leq 3\ln|\zeta(\sigma)| + 4\ln|\zeta(\sigma+it)| + \ln|\zeta(\sigma+2it)|

= \ln|\zeta(\sigma)^{3}\zeta(\sigma+it)^{4}\zeta(\sigma+2it)|,

可以推导出 |\zeta(\sigma)^{3}\zeta(\sigma+it)^{4}\zeta(\sigma+2it)|\geq 1. 因此 Claim 2 证明完毕。

Claim 3. \zeta(1+it)\neq 0 对于所有的 \{t\in\mathbb{R}: t\neq 0\} 成立。

反证法:假设 \zeta(s)s=\sigma + it (t\neq 0) 存在阶数为 m 的零点。也就是说:

\lim_{\sigma\rightarrow 1^{+}} \frac{\zeta(\sigma+it)}{(\sigma+it-1)^{m}}=c\neq 0, 其中 m\geq 1.

从 Riemann Zeta 函数的延拓可以知道,\lim_{\sigma\rightarrow 1^{+}}(\sigma -1)\zeta(\sigma) = 1. 并且 \zeta(s)\{s\in\mathbb{C}:\Re(s)>0, s\neq 1\} 上是解析函数。

从 Claim 2 可以得到:

|(\sigma-1)^{3}\zeta(\sigma)^{3}(\sigma+it-1)^{-4m}\zeta(\sigma+it)^{4}\zeta(\sigma+2it)|

\geq |\sigma-1|^{3}|\sigma-1+it|^{-4m}

\geq |\sigma-1|^{3}\cdot |\sigma-1|^{-4m}

= \frac{1}{|\sigma-1|^{4m-3}}.

\sigma\rightarrow 1^{+}, 可以得到左侧趋近于一个有限的值,但是右侧趋近于无穷,所以得到矛盾。也就是说当 t\neq 0 时, \zeta(1+it)\neq 0 成立。

根据之前的知识,s= 1\zeta(s) 的极点,所以我们得到了本篇文章的主要结论:\zeta(s)\{s\in\mathbb{C}:\Re(s)\geq 1\} 上面没有零点。

 

总结

本篇文章从 Riemann Zeta 函数的延拓开始,证明了 Riemann Zeta 函数在 \{s\in\mathbb{C}:\Re(s)\geq 1\} 上没有零点。在下一篇文章中,笔者将会证明在 \Re(s)=1 附近一个“狭长”的区域上,Riemann Zeta 函数没有零点。

 

从调和级数到 RIEMANN ZETA 函数(一)

Riemann Zeta 函数

Riemann Zeta 函数(Riemann zeta function),\zeta(s),是一个关于复数 s 的方程。在复平面上,当复数 s 的实数部分 \sigma=\Re s >1 时,\zeta(s) 就是如下的级数形式:

\zeta(s) = \sum_{n=1}^{\infty}\frac{1}{n^{s}}.

调和级数的概念与性质

既然提到了级数,首先让我们来回顾一下级数的定义是什么?

级数的定义:在数学中,一个有穷或者无穷的序列 (x_{0},x_{1},x_{2},...) 的形式和 S = x_{0}+x_{1}+x_{2}+... 称为级数,里面的每一项都称为级数的通项。

级数收敛的定义:令 S_{n}=x_{0}+...+x_{n},如果存在有限的 S 使得 \lim_{n\rightarrow \infty}S_{n}=S,那么就称该级数收敛。否则,该级数就称为发散级数。

然后下面我们来研究一下调和级数的基本性质。调和级数的表达式写出来十分简单,那就是 Riemann Zeta 函数在 s=1 的取值,i.e.

\zeta(1) = \sum_{n=1}^{+\infty}\frac{1}{n}.

提到级数的收敛或发散,就必须要提到关于级数收敛的等价定理(Cauchy 判别法),那就是:级数 S_{n} 收敛当且仅当对任意的 \epsilon>0,存在 N 使得对于任意的 m, n>N 都有 |S_{m}-S_{n}|<\epsilon.

既然是等价定理,那么就可以使用 Cauchy 判别法来判断调和级数是否收敛。

Method 1.

S_{n}=\sum_{k=1}^{n}\frac{1}{k},

直接通过计算得到

|S_{2n}-S_{n}|=\frac{1}{n+1}+...+\frac{1}{2n}>\frac{1}{2n}+...+\frac{1}{2n}=\frac{1}{2},

说明该级数是不收敛的,也就是调和级数是发散的。

除了基于 Cauchy 收敛准则的证明之外,能否写出判断调和级数发散的其他方法呢?答案是肯定的。以下有一种使用初等数学方法就能够解释调和级数发散的方法。

Method 2.

\sum_{n=1}^{+\infty}\frac{1}{n}

=1+\frac{1}{2}+(\frac{1}{3}+\frac{1}{4})+(\frac{1}{5}+\frac{1}{6}+\frac{1}{7}+\frac{1}{8})+...

>1+\frac{1}{2}+(\frac{1}{4}+\frac{1}{4})+(\frac{1}{8}+\frac{1}{8}+\frac{1}{8}+\frac{1}{8})+...

=1+\frac{1}{2}+\frac{1}{2}+\frac{1}{2}+...=+\infty.

既然都提到了高等数学,那么当然不能仅仅局限于使用初等数学的技巧来解决问题。而且如果只是用初等数学的方法,在拓展性方面就会受到极大的限制。

Method 3. 调和级数的发散可以通过定积分的技巧来进行解决。

HarmonicSeries

1+\frac{1}{2}+...+\frac{1}{n}

>\int_{1}^{2}\frac{1}{x}dx + \int_{2}^{3}\frac{1}{x}dx+...+\int_{n}^{n+1}\frac{1}{x}dx

=\int_{1}^{n+1}\frac{1}{x}dx=\ln(n+1)

因此,\sum_{n=1}^{\infty}\frac{1}{n}=+\infty.

从上面的定积分的方法可以预计出调和级数的量级大约是对数的量级,那么能否精确的估计出来呢?例如下面这个问题:

问题:\lim_{n\rightarrow +\infty}\frac{\sum_{k=1}^{n}\frac{1}{k}}{\ln(n)}=?

通过 L’Hospital 法则可知:\lim_{x\rightarrow 0}x/\ln(1+x)=1.

通过 Stolz 定理可知:

\lim_{n\rightarrow +\infty}\frac{\sum_{k=1}^{n}\frac{1}{k}}{\ln(n)}

= \lim_{n\rightarrow +\infty}\frac{\frac{1}{n}}{\ln(n/(n-1))}

= \lim_{x\rightarrow 0}\frac{x}{\ln(1+x)}=1

除此之外,我们同样可以证明

\lim_{n\rightarrow+\infty}(1+\frac{1}{2}+...+\frac{1}{n}-\ln(n))

这个极限是存在并且有限的。

调和级数的推广

那么,如果在考虑 \zeta(2) 也就是级数

\zeta(2) = \sum_{n=1}^{\infty}\frac{1}{n^{2}}

是否收敛的时候,能否用到以上类似的技巧呢?首先,确实也存在各种各样的初等数学技巧,例如:

Method 1.

\sum_{n=1}^{+\infty}\frac{1}{n^{2}}<1+\sum_{n=2}^{+\infty}\frac{1}{n(n-1)}=1+\sum_{n=2}^{+\infty}(\frac{1}{n-1}-\frac{1}{n})=2.

Method 2. 使用数学归纳法。也就是要证明:

\sum_{k=1}^{n}1/k^{2}\leq 2-\frac{1}{n}.

n=1 的时候,公式是正确的。假设 n 的时候是正确的,那么我们有\sum_{k=1}^{n}1/k^{2}\leq 2-\frac{1}{n}。计算可得:

\sum_{k=1}^{n+1}\frac{1}{k^{2}}

<2-\frac{1}{n}+\frac{1}{(n+1)^{2}}

= 2- \frac{1}{n+1}-\frac{1}{n(n+1)^{2}}

\leq 2-\frac{1}{n+1}.

因此,不等式正确,所以 \sum_{n=1}^{+\infty}1/n^{2} 收敛。

其次,在判断调和级数发散的时候,使用的定积分的方法同样可以应用在这个场景下。

Method 3.

1+\frac{1}{2^{2}}+...+\frac{1}{n^{2}}

<1+\int_{1}^{2}\frac{1}{x^{2}}dx+...+\int_{n-1}^{n}\frac{1}{x^{2}}dx

=1+\int_{1}^{n}\frac{1}{x^{2}}dx=1+1-\frac{1}{n}<2.

那么这个是针对次数等于2的情况,对于一般的情形,

\zeta(s)=\sum_{n=1}^{+\infty}\frac{1}{n^{s}},\sigma = \Re(s)>1.

使用定积分的技术,同样可以证明对于任意的 \sigma = \Re(s)>1,都有 \zeta(s) 是收敛的。但是 \zeta(1) 是发散的。

Riemann Zeta 函数中某些点的取值

除此之外,既然 \zeta(s)\sigma = \Re(s)>1 的时候收敛,能否计算出某些函数的特殊值呢?答案是肯定的,例如,我们可以使用 Fourier 级数来计算出 \zeta(2), \zeta(4), \zeta(6),... 的取值。首先,我们回顾一下 Fourier 级数的一些性质:

假设 f(x) 是一个关于 2\pi 的周期函数, i.e. f(x)=f(x+2\pi) 对于所有的 x \in \mathbb{R} 都成立。那么函数 f(x) 的 Fourier 级数就定义为

a_{0}+\sum_{n=1}^{\infty} (a_{n} \cos(nx) +b_{n} \sin(nx)),

其中,a_{0}= \frac{1}{2\pi} \int_{-\pi}^{\pi} f(x) dx,

a_{n}= \frac{1}{\pi} \int_{-\pi}^{\pi} f(x) \cos(nx) dx n\geq 1,

b_{n}= \frac{1}{\pi} \int_{-\pi}^{\pi} f(x) \sin(nx) dx n\geq 1,

定理 1. 如果 f(x) 在区间 (-\pi, \pi) 上满足 Lipschitz 条件,那么

f(x) =a_{0}+\sum_{n=1}^{\infty} (a_{n} \cos(nx) +b_{n} \sin(nx)).

定理 2. Parseval’s 恒等式.

\frac{1}{\pi} \int_{-\pi}^{\pi} |f(x)|^{2} dx= 2a_{0}^{2}+ \sum_{n=1}^{\infty} (a_{n}^{2}+b_{n}^{2}).

下面我们就来证明下列恒等式:

\sum_{n=1}^{\infty} \frac{1}{(2n-1)^{2}}=\frac{\pi^{2}}{8}

\sum_{n=1}^{\infty} \frac{1}{n^{2}}=\frac{\pi^{2}}{6}

\sum_{n=1}^{\infty} \frac{1}{(2n-1)^{4}}=\frac{\pi^{4}}{96}

\sum_{n=1}^{\infty} \frac{1}{n^{4}}=\frac{\pi^{4}}{90}

证明:

选择在区间 (-\pi, \pi) 上的函数 f(x)=|x|,并且该函数是关于 2\pi 的周期函数。

使用 a_{n}b_{n} 的公式,我们可以得到函数 f(x)=|x| 的 Fourier 级数是

\frac{\pi}{2} + \sum_{n=1}^{\infty} \frac{2((-1)^{n}-1)}{\pi} \cdot \frac{cos(nx)}{n^{2}}

从定理1, 令 x=0, 可以得到

0= \frac{\pi}{2} + \sum_{n=1}^{\infty} \frac{2((-1)^{n}-1)}{n^{2} \pi} = \frac{\pi}{2} + \sum_{m=1}^{\infty} \frac{-4}{(2m-1)^{2}\pi} = \frac{\pi}{2} - \frac{4}{\pi} \sum_{m=1}^{\infty} \frac{1}{(2m-1)^{2}}

因此,\sum_{n=1}^{\infty} \frac{1}{(2n-1)^{2}}=\frac{\pi^{2}}{8} .

假设 S=\sum_{n=1}^{\infty} \frac{1}{n^{2}} , 可以得到

S=\sum_{odd} \frac{1}{n^{2}} + \sum_{even} \frac{1}{n^{2}} = \frac{\pi^{2}}{8} + \frac{1}{4} S .

因此 S=\frac{\pi^{2}}{6} .

从 Parserval’s 恒等式,我们知道

\frac{2\pi^{2}}{3}= \frac{1}{\pi} \int_{-\pi}^{\pi} x^{2}dx = 2\cdot (\frac{\pi}{2})^{2} + \sum_{n=1}^{\infty} \frac{4((-1)^{n}-1)^{2}}{\pi^{2}\cdot n^{4}} = \frac{\pi^{2}}{2} + \sum_{m=1}^{\infty} \frac{16}{\pi^{2} (2m-1)^{4}}

因此 \sum_{n=1}^{\infty} \frac{1}{(2n-1)^{4}} = \frac{\pi^{4}}{96} .

假设 S=\sum_{n=1}^{\infty} \frac{1}{n^{4}} , 得到

S=\sum_{odd} \frac{1}{n^{4}} + \sum_{even} \frac{1}{n^{4}} = \frac{\pi^{4}}{96} + \frac{1}{16} S

因此, S=\frac{\pi^{4}}{90} .

总结

本篇文章从调和级数的发散性开始,介绍了判断调和级数是否收敛的几种方法。进一步考虑了其他级数的收敛性,并通过 Fourier 级数的方法计算出了部分 Riemann Zeta 函数的取值。

Hausdorff dimension of the graphs of the classical Weierstrass functions

In this paper, we obtain the explicit value of the Hausdorff dimension of the graphs of the classical Weierstrass functions, by proving absolute continuity of the SRB measures of the associated solenoidal attractors.

1. Introduction

In Real Analysis, the classical Weierstrass function is

\displaystyle W_{\lambda,b}(x) = \sum\limits_{n=0}^{\infty} \lambda^n \cos(2\pi b^n x)

with {1/b < \lambda < 1}.

Note that the Weierstrass functions have the form

\displaystyle f^{\phi}_{\lambda,b}(x) = \sum\limits_{n=0}^{\infty} \lambda^n \phi(b^n x)

where {\phi} is a {\mathbb{Z}}-periodic {C^2}-function.

Weierstrass (1872) and Hardy (1916) were interested in {W_{\lambda,b}} because they are concrete examples of continuous but nowhere differentiable functions.

Remark 1 The graph of {f^{\phi}_{\lambda,b}} tends to be a “fractal object” because {f^{\phi}_{\lambda,b}} is self-similar in the sense that

\displaystyle f^{\phi}_{\lambda, b}(x) = \phi(x) + \lambda f^{\phi}_{\lambda,b}(bx)

We will come back to this point later.

Remark 2 {f^{\phi}_{\lambda,b}} is a {C^{\alpha}}-function for all {0\leq \alpha < \frac{-\log\lambda}{\log b}}. In fact, for all {x,y\in[0,1]}, we have

\displaystyle \frac{f^{\phi}_{\lambda, b}(x) - f^{\phi}_{\lambda,b}(y)}{|x-y|^{\alpha}} = \sum\limits_{n=0}^{\infty} \lambda^n b^{n\alpha} \left(\frac{\phi(b^n x) - \phi(b^n y)}{|b^n x - b^n y|^{\alpha}}\right),

so that

\displaystyle \frac{f^{\phi}_{\lambda, b}(x) - f^{\phi}_{\lambda,b}(y)}{|x-y|^{\alpha}} \leq \|\phi\|_{C^{\alpha}} \sum\limits_{n=0}^{\infty}(\lambda b^{\alpha})^n:=C(\phi,\alpha,\lambda,b) < \infty

whenever {\lambda b^{\alpha} < 1}, i.e., {\alpha < -\log\lambda/\log b}.

The study of the graphs of {W_{\lambda,b}} as fractal sets started with the work of Besicovitch-Ursell in 1937.

Remark 3 The Hausdorff dimension of the graph of a {C^{\alpha}}-function {f:[0,1]\rightarrow\mathbb{R}}is

\displaystyle \textrm{dim}(\textrm{graph}(f))\leq 2 - \alpha

Indeed, for each {n\in\mathbb{N}}, the Hölder continuity condition

\displaystyle |f(x)-f(y)|\leq C|x-y|^{\alpha}

leads us to the “natural cover” of {G=\textrm{graph}(f)} by the family {(R_{j,n})_{j=1}^n} of rectangles given by

\displaystyle R_{j,n}:=\left[\frac{j-1}{n}, \frac{j}{n}\right] \times \left[f(j/n)-\frac{C}{n^{\alpha}}, f(j/n)+\frac{C}{n^{\alpha}}\right]

Nevertheless, a direct calculation with the family {(R_{j,n})_{j=1}^n} does not give us an appropriate bound on {\textrm{dim}(G)}. In fact, since {\textrm{diam}(R_{j,n})\leq 4C/n^{\alpha}} for each {j=1,\dots, n}, we have

\displaystyle \sum\limits_{j=1}^n\textrm{diam}(R_{j,n})^d\leq n\left(\frac{4C}{n^{\alpha}}\right)^d = (4C)^{1/\alpha} < \infty

for {d=1/\alpha}. Because {n\in\mathbb{N}} is arbitrary, we deduce that {\textrm{dim}(G)\leq 1/\alpha}. Of course, this bound is certainly suboptimal for {\alpha<1/2} (because we know that {\textrm{dim}(G)\leq 2 < 1/\alpha} anyway).Fortunately, we can refine the covering {(R_{j,n})} by taking into account that each rectangle {R_{j,n}} tends to be more vertical than horizontal (i.e., its height {2C/n^{\alpha}} is usually larger than its width {1/n}). More precisely, we can divide each rectangle {R_{j,n}} into {\lfloor n^{1-\alpha}\rfloor} squares, say

\displaystyle R_{j,n} = \bigcup\limits_{k=1}^{\lfloor n^{1-\alpha}\rfloor}Q_{j,n,k},

such that every square {Q_{j,n,k}} has diameter {\leq 2C/n}. In this way, we obtain a covering {(Q_{j,n,k})} of {G} such that

\displaystyle \sum\limits_{j=1}^n\sum\limits_{k=1}^{\lfloor n^{1-\alpha}\rfloor} \textrm{diam}(Q_{j,n,k})^d \leq n\cdot n^{1-\alpha}\cdot\left(\frac{2}{n}\right)^d\leq (2C)^{2-\alpha}<\infty

for {d=2-\alpha}. Since {n\in\mathbb{N}} is arbitrary, we conclude the desired bound

\displaystyle \textrm{dim}(G)\leq 2-\alpha

A long-standing conjecture about the fractal geometry of {W_{\lambda,b}} is:

Conjecture (Mandelbrot 1977): The Hausdorff dimension of the graph of {W_{\lambda,b}} is

\displaystyle 1<\textrm{dim}(\textrm{graph}(W_{\lambda,b})) = 2 + \frac{\log\lambda}{\log b} < 2

Remark 4 In view of remarks 2 and 3, the whole point of Mandelbrot’s conjecture is to establish the lower bound

\displaystyle \textrm{dim}(\textrm{graph}(W_{\lambda,b})) \geq 2 + \frac{\log\lambda}{\log b}

Remark 5 The analog of Mandelbrot conjecture for the box and packing dimensions is known to be true: see, e.g., these papers here and here).

In a recent paper (see here), Shen proved the following result:

Theorem 1 (Shen) For any {b\geq 2} integer and for all {1/b < \lambda < 1}, the Mandelbrot conjecture is true, i.e.,

\displaystyle \textrm{dim}(\textrm{graph}(W_{\lambda,b})) = 2 + \frac{\log\lambda}{\log b}

Remark 6 The techniques employed by Shen also allow him to show that given {\phi:\mathbb{R}\rightarrow\mathbb{R}} a {\mathbb{Z}}-periodic, non-constant, {C^2} function, and given {b\geq 2} integer, there exists {K=K(\phi,b)>1} such that

\displaystyle \textrm{dim}(\textrm{graph}(f^{\phi}_{\lambda,b})) = 2 + \frac{\log\lambda}{\log b}

for all {1/K < \lambda < 1}.

Remark 7 A previous important result towards Mandelbrot’s conjecture was obtained by Barańsky-Barány-Romanowska (in 2014): they proved that for all {b\geq 2} integer, there exists {1/b < \lambda_b < 1} such that

\displaystyle \textrm{dim}(\textrm{graph}(W_{\lambda,b})) = 2 + \frac{\log\lambda}{\log b}

for all {\lambda_b < \lambda < 1}.

The remainder of this post is dedicated to give some ideas of Shen’s proof of Theorem1 by discussing the particular case when {1/b<\lambda<2/b} and {b\in\mathbb{N}} is large.

2. Ledrappier’s dynamical approach

If {b\geq 2} is an integer, then the self-similar function {f^{\phi}_{\lambda,b}} (cf. Remark 1) is also {\mathbb{Z}}-periodic, i.e., {f^{\phi}_{\lambda,b}(x+1) = f^{\phi}_{\lambda,b}(x)} for all {x\in\mathbb{R}}. In particular, if {b\geq 2} is an integer, then {\textrm{graph}(f^{\phi}_{\lambda,b})} is an invariant repeller for the endomorphism {\Phi:\mathbb{R}/\mathbb{Z}\times\mathbb{R}\rightarrow \mathbb{R}/\mathbb{Z}\times\mathbb{R}} given by

\displaystyle \Phi(x,y) = \left(bx\textrm{ mod }1, \frac{y-\phi(x)}{\lambda}\right)

This dynamical characterization of {G = \textrm{graph}(f^{\phi}_{\lambda,b})} led Ledrappier to the following criterion for the validity of Mandelbrot’s conjecture when {b\geq 2} is an integer.

Denote by {\mathcal{A}} the alphabet {\mathcal{A}=\{0,\dots,b-1\}}. The unstable manifolds of {\Phi}through {G} have slopes of the form

\displaystyle (1,-\gamma \cdot s(x,u))

where {\frac{1}{b} < \gamma = \frac{1}{\lambda b} <1}, {x\in\mathbb{R}}, {u\in\mathcal{A}^{\mathbb{N}}}, and

\displaystyle s(x,u):=\sum\limits_{n=0}^{\infty} \gamma^n \phi'\left(\frac{x + u_1 + u_2 b + \dots + u_n b^{n-1}}{b^n}\right)

In this context, the push-forwards {m_x := (u\mapsto s(x,u))_*\mathbb{P}} of the Bernoulli measure {\mathbb{P}} on {\mathcal{A}^{\mathbb{N}}} (induced by the discrete measure assigning weight {1/b} to each letter of the alphabet {\mathcal{A}}) play the role of conditional measures along vertical fibers of the unique Sinai-Ruelle-Bowen (SRB) measure {\theta} of the expanding endomorphism {T:\mathbb{R}/\mathbb{Z}\times\mathbb{R} \rightarrow \mathbb{R}/\mathbb{Z}\times\mathbb{R}},

\displaystyle T(x,y) = (bx\textrm{ mod }1, \gamma y + \psi(x)),

where {\gamma=1/\lambda b} and {\psi(x)=\phi'(x)}. In plain terms, this means that

\displaystyle \theta = \int_{\mathbb{R}/\mathbb{Z}} m_x \, d\textrm{Leb}(x) \ \ \ \ \ (1)

where {\theta} is the unique {T}-invariant probability measure which is absolutely continuous along unstable manifolds (see Tsujii’s paper).

As it was shown by Ledrappier in 1992, the fractal geometry of the conditional measures {m_x} have important consequences for the fractal geometry of the graph {G}:

Theorem 2 (Ledrappier) Suppose that for Lebesgue almost every {x\in\mathbb{R}} the conditional measures {m_x} have dimension {\textrm{dim}(m_x)=1}, i.e.,

\displaystyle \lim\limits_{r\rightarrow 0}\frac{\log m_x(B(z,r))}{\log r} = 1 \textrm{ for } m_x\textrm{-a.e. } z

Then, the graph {G=\textrm{graph}(f^{\phi}_{\lambda,b})} has Hausdorff dimension

\displaystyle \textrm{dim}(G) = 2 + \frac{\log\lambda}{\log b}

Remark 8 Very roughly speaking, the proof of Ledrappier theorem goes as follows. By Remark 4, it suffices to prove that {\textrm{dim}(G)\geq 2 + \frac{\log\lambda}{\log b}}. By Frostman lemma, we need to construct a Borel measure {\nu} supported on {G} such that

\displaystyle \underline{\textrm{dim}}(\nu) := \textrm{ ess }\inf \underline{d}(\nu,x) \geq 2 + \frac{\log\lambda}{\log b}

where {\underline{d}(\nu,x):=\liminf\limits_{r\rightarrow 0}\log \nu(B(x,r))/\log r}. Finally, the main point is that the assumptions in Ledrappier theorem allow to prove that the measure {\mu^{\phi}_{\lambda, b}} given by the lift to {G} of the Lebesgue measure on {[0,1]} via the map {x\mapsto (x,f^{\phi}_{\lambda,b}(x))}satisfies

\displaystyle \underline{\textrm{dim}}(\mu^{\phi}_{\lambda,b}) \geq 2 + \frac{\log\lambda}{\log b}

An interesting consequence of Ledrappier theorem and the equation 1 is the following criterion for Mandelbrot’s conjecture:

Corollary 3 If {\theta} is absolutely continuous with respect to the Lebesgue measure {\textrm{Leb}_{\mathbb{R}^2}}, then

\displaystyle \textrm{dim}(G) = 2 + \frac{\log\lambda}{\log b}

Proof: By (1), the absolute continuity of {\theta} implies that {m_x} is absolutely continuous with respect to {\textrm{Leb}_{\mathbb{R}}} for Lebesgue almost every {x\in\mathbb{R}}.

Since {m_x\ll \textrm{Leb}_{\mathbb{R}}} for almost every {x} implies that {\textrm{dim}(m_x)=1} for almost every {x}, the desired corollary now follows from Ledrappier’s theorem. \Box

3. Tsujii’s theorem

The relevance of Corollary 3 is explained by the fact that Tsujii found an explicittransversality condition implying the absolute continuity of {\theta}.

More precisely, Tsujii firstly introduced the following definition:

Definition 4

  • Given {\varepsilon>0}, {\delta>0} and {x_0\in\mathbb{R}/\mathbb{Z}}, we say that two infinite words {u, v\in\mathcal{A}^{\mathbb{N}}} are {(\varepsilon,\delta)}-transverse at {x_0} if either

    \displaystyle |s(x_0,u)-s(x_0,v)|>\varepsilon

    or

    \displaystyle |s'(x_0,u)-s'(x_0,v)|>\delta

  • Given {q\in\mathbb{N}}, {\varepsilon>0}, {\delta>0} and {x_0\in\mathbb{R}/\mathbb{Z}}, we say that two finite words {k,l\in\mathcal{A}^q} are {(\varepsilon,\delta)}-transverse at {x_0} if {ku}, {lv} are {(\varepsilon,\delta)}-transverse at {x_0}for all pairs of infinite words {u,v\in\mathcal{A}^{\mathbb{N}}}; otherwise, we say that {k} and {l} are{(\varepsilon,\delta)}-tangent at {x_0};
  • {E(q,x_0;\varepsilon,\delta):= \{(k,l)\in\mathcal{A}^q\times\mathcal{A}^q: (k,l) \textrm{ is } (\varepsilon,\delta)\textrm{-tangent at } x_0\}}
  • {E(q,x_0):=\bigcap\limits_{\varepsilon>0}\bigcap\limits_{\delta>0} E(q,x_0;\varepsilon,\delta)};
  • {e(q,x_0):=\max\limits_{k\in\mathcal{A}^q}\#\{l\in\mathcal{A}^q: (k,l)\in E(q,x_0)\}}
  • {e(q):=\max\limits_{x_0\in\mathbb{R}/\mathbb{Z}} e(q,x_0)}.

Next, Tsujii proves the following result:

Theorem 5 (Tsujii) If there exists {q\geq 1} integer such that {e(q)<(\gamma b)^q}, then

\displaystyle \theta\ll\textrm{Leb}_{\mathbb{R}^2}

Remark 9 Intuitively, Tsujii’s theorem says the following. The transversality condition {e(q)<(\gamma b)^q} implies that the majority of strong unstable manifolds {\ell^{uu}}are mutually transverse, so that they almost fill a small neighborhood {U} of some point {x_0} (see the figure below extracted from this paper of Tsujii). Since the SRB measure {\theta} is absolutely continuous along strong unstable manifolds, the fact that the {\ell^{uu}}‘s almost fill {U} implies that {\theta} becomes “comparable” to the restriction of the Lebesgue measure {\textrm{Leb}_{\mathbb{R}^2}} to {U}.

tsujiiacta

Remark 10 In this setting, Barańsky-Barány-Romanowska obtained their main result by showing that, for adequate choices of the parameters {\lambda} and {b}, one has {e(1)=1}. Indeed, once we know that {e(1)=1}, since {1<\gamma b}, they can apply Tsujii’s theorem and Ledrappier’s theorem (or rather Corollary 3) to derive the validity of Mandelbrot’s conjecture for certain parameters {\lambda} and {b}.

For the sake of exposition, we will give just a flavor of the proof of Theorem 1 by sketching the derivation of the following result:

Proposition 6 Let {\phi(x) = \cos(2\pi x)}. If {1/2<\gamma=1/\lambda b <1} and {b\in\mathbb{N}} is sufficiently large, then

\displaystyle e(1)<\gamma b

In particular, by Corollary 3 and Tsujii’s theorem, if {1/2<\gamma=1/\lambda b <1} and {b\in\mathbb{N}} is sufficiently large, then Mandelbrot’s conjecture is valid, i.e.,

\displaystyle \textrm{dim}(W_{\lambda,b}) = 2+\frac{\log\lambda}{\log b}

Remark 11 The proof of Theorem 1 in full generality (i.e., for {b\geq 2} integer and {1/b<\lambda<1}) requires the introduction of a modified version of Tsujii’s transversality condition: roughly speaking, Shen defines a function {\sigma(q)\leq e(q)}(inspired from Peter-Paul inequality) and he proves

  • (a) a variant of Proposition 6: if {b\geq 2} integer and {1/b<\lambda<1}, then {\sigma(q)<(\gamma b)^q} for some integer {q};
  • (b) a variant of Tsujii’s theorem: if {\sigma(q)<(\gamma b)^q} for some integer {q}, then {\theta\ll\textrm{Leb}_{\mathbb{R}^2}}.

See Sections 2, 3, 4 and 5 of Shen’s paper for more details.

We start the (sketch of) proof of Proposition 6 by recalling that the slopes of unstable manifolds are given by

\displaystyle s(x,u):=-2\pi\sum\limits_{n=0}^{\infty} \gamma^n \sin\left(2\pi\frac{x + u_1 + u_2 b + \dots + u_n b^{n-1}}{b^n}\right)

for {x\in\mathbb{R}}, {u\in\mathcal{A}^{\mathbb{N}}}, so that

\displaystyle s'(x,u)=-4\pi^2\sum\limits_{n=0}^{\infty} \left(\frac{\gamma}{b}\right)^n \cos\left(2\pi\frac{x + u_1 + u_2 b + \dots + u_n b^{n-1}}{b^n}\right)

Remark 12 Since {\gamma/b < \gamma}, the series defining {s'(x,u)} converges faster than the series defining {s(x,u)}.

By studying the first term of the expansion of {s(x,u)} and {s'(x,u)} (while treating the remaining terms as a “small error term”), it is possible to show that if {(k,l)\in E(1,x_0)}, then

\displaystyle \left|\sin\left(2\pi\frac{x_0+k}{b}\right) - \sin\left(2\pi\frac{x_0+l}{b}\right)\right| \leq\frac{2\gamma}{1-\gamma} \ \ \ \ \ (2)

and

\displaystyle \left|\cos\left(2\pi\frac{x_0+k}{b}\right) - \cos\left(2\pi\frac{x_0+l}{b}\right)\right| \leq \frac{2\gamma}{b-\gamma} \ \ \ \ \ (3)

(cf. Lemma 3.2 in Shen’s paper).

Using these estimates, we can find an upper bound for {e(1)} as follows. Take {x_0\in\mathbb{R}/\mathbb{Z}} with {e(1)=e(1,x_0)}, and let {k\in\mathcal{A}} be such that {(k,l_1),\dots,(k,l_{e(1)})\in E(1,x_0)} distinct elements listed in such a way that

\displaystyle \sin(2\pi x_i)\leq \sin(2\pi x_{i+1})

for all {i=1,\dots,e(1)-1}, where {x_i:=(x_0+l_i)/b}.

From (3), we see that

\displaystyle \left|\cos\left(2\pi x_i\right) - \cos\left(2\pi x_{i+1}\right)\right| \leq \frac{4\gamma}{b-\gamma}

for all {i=1,\dots,e(1)-1}.

Since

\displaystyle (\cos(2\pi x_i)-\cos(2\pi x_{i+1}))^2 + (\sin(2\pi x_i)-\sin(2\pi x_{i+1}))^2 = 4\sin^2(\pi(x_i-x_{i+1}))\geq 4\sin^2(\pi/b),

it follows that

\displaystyle |\sin(2\pi x_i)-\sin(2\pi x_{i+1})|\geq \sqrt{4\sin^2\left(\frac{\pi}{b}\right) - \left(\frac{4\gamma}{b-\gamma}\right)^2} \ \ \ \ \ (4)

Now, we observe that

\displaystyle \sqrt{4\sin^2\left(\frac{\pi}{b}\right) - \left(\frac{4\gamma}{b-\gamma}\right)^2} > \frac{4}{b} \ \ \ \ \ (5)

for {b} large enough. Indeed, this happens because

  • {\sqrt{z^2-w^2}>2(z-w)} if {z+w>4(z-w)};
  • {z+w>4(z-w)} if {z/w:=u < 5/3};
  • {\frac{2\sin(\frac{\pi}{b})}{\frac{4\gamma}{b-\gamma}}\rightarrow \frac{2\pi}{4\gamma} (< \frac{5}{3})} as {b\rightarrow\infty}, and {2\sin(\frac{\pi}{b}) - \frac{4\gamma}{b-\gamma} \rightarrow (2\pi-4\gamma)\frac{1}{b} (>\frac{2}{b})} as {b\rightarrow\infty} (here we used {\gamma<1}).

By combining (4) and (5), we deduce that

\displaystyle |\sin(2\pi x_i)-\sin(2\pi x_{i+1})| > 4/b

for all {i=1,\dots, e(1)-1}.

Since {-1\leq\sin(2\pi x_1)\leq\sin(2\pi x_2)\leq\dots\leq\sin(2\pi x_{e(1)})\leq 1}, the previous estimate implies that

\displaystyle \frac{4}{b}(e(1)-1)<\sum\limits_{i=1}^{e(1)-1}(\sin(2\pi x_{i+1}) - \sin(2\pi x_i)) = \sin(2\pi x_{e(1)}) - \sin(2\pi x_1)\leq 2,

i.e.,

\displaystyle e(1)<1+\frac{b}{2}

Thus, it follows from our assumptions ({\gamma>1/2}, {b} large) that

\displaystyle e(1)<1+\frac{b}{2}<\gamma b

This completes the (sketch of) proof of Proposition 6 (and our discussion of Shen’s talk).