科学空间|Scientific Spaces

牛顿法是求方程近似根的一个相当有用而且快捷的方法，我们最近科学计算软件课程（Matlab）的一个作业就是编写求方程近似解的程序，其中涉及到牛顿法。我们要实现的目标是，用户输入一道方程，脚本就自动求出根来。这看起来是一个挺简单的循环迭代程序，但是由于Matlab本身的特殊性，却产生了不少困难。

Matlab是为了数值计算（尤其是矩阵运算）而生的，因此它并不擅长处理符号计算。这就给我们编程带来了困难。在网上随便一搜，就可以发现，网上的Matlab牛顿法程序都是要求用户同时输入方程及其导函数，这显然是不方便的，因为Matlab本身就具备了求导功能。下面我们来分析一下困难在哪里。

我们要实现的最基本功能是定义一个函数，然后可以根据该函数求具体的函数值，并且自动求该函数的导数，接着求导数值。这些看起来很基本的功能在Matlab中却很难调和，因为Matlab的“函数”定义很广，一个具有特定功能的M文件叫“函数”，一个运算式$f(x)$也可能是一个函数，显然后者是可以求导的，前者却不行，所以Matlab一刀砍——不能对函数求导！！

点击阅读全文...

进入期末了，很多同学都开始复习了，这学期我选的几门课到现在还不是很熟悉，本想也在趁着这段时间好好看看。偏生五天前我在浏览数学研发论坛的编程擂台时看到了这样的一道题目：

设对于给定的$L$，方程
$$\frac{1}{x}+\frac{1}{y}=\frac{1}{n}$$
满足$0 < x < y \leq L$的正整数解共有$f(L)$种情况。比如$f(6)=1,f(12)=3,f(1000)=1069$，求$f(10^{12})$。

这道题目的来源是Project Euler的第454题：Diophantine reciprocals III（丢潘图倒数方程），题目简短易懂，但又不失深度，正符合我对理想题目的定义。而且最近在学习Python学习得不亦乐乎，看到这道题目就跃跃欲试。于是乎，我的五天时间就没有了，而且过程中几乎耗尽了我现在懂的所有编程技巧。由于不断地测试运行，我的电脑发热量比平时大了几倍，真是辛苦了我的电脑。最后的代码，自我感觉已经是我目前写的最精彩的代码了。在此与大家共享和共勉~

上述表达式是分式，不利于编程，由于$n=\frac{xy}{x+y}$，于是上述题目也等价于求$(x+y)|xy$（意思是$x+y$整除$xy$）的整数解。

点击阅读全文...

Mathieu方程

在文章《有质动力：倒立单摆的稳定性》中，我们分析了通过高频低幅振荡来使得倒立单摆稳定的可能性，并且得出了运动方程
$$l\ddot{\theta}+[h_0 \omega^2 \cos(\omega t)-g]\sin\theta=0$$

由此对单摆频率的下限提出了要求$\omega \gg \sqrt{\frac{g}{h_0}}$。然而，那个下限只不过是必要的，却不是充分的。如果要完整地分析该单摆的运动方程，最理想的方法当然是写出上述常微分方程的解析解。不过很遗憾，我们并没有办法做到这一点。我们只能够采取各种近似方法来求解。近似方法一般指数值计算方法，然后笔者偏爱的是解析方法，也就是说，即使是近似解，也希望能够求出近似的解析解。

点击阅读全文...

在对数学或物理进行事后分析，往往会发现一些奇怪的现象，也有可能得到一些更为深刻有趣的结果。比如本文所要谈及的傅里叶变换，可以由一种“异想天开”的思路得来。

洛朗展式

我们知道，在原点处形态良好的函数，可以展开为泰勒级数
$$f(x)=\sum_{n=0}^{\infty}a_n x^n$$
我们发现，上面的幂都是正的，为什么不能包含$x$的负数次幂呢？比如$\frac{\sin z}{z^2}$展开为
$$\frac{1}{z}-\frac{z}{6}+\frac{z^3}{120}\dots$$
显然也是一件合理的事情。于是，结合复变函数，我们得到解析函数的洛朗展式
$$f(z)=\sum_{n=-\infty}^{+\infty}a_n z^n$$
这是函数的双边展开。其中

点击阅读全文...

BoJone很喜欢Python，也很喜欢数论，所以就喜欢利用Python玩数论了。平时也喜欢自己动手写一些数论函数，毕竟Python支持大整数高精度运算，这点是非常好的；但是，在很多实际应用中，还是希望能有一个现成的数论函数库来调用。之前尝试过数学研发网的HugeCalc库，但是由于各种不熟悉不了了之。后来论坛上的无心老兄推荐了PARI/GP，小试一下，居然在Python上成功调用了。以后再也不用担心Python上的数论计算问题了，呵呵~

点击阅读全文...

ODE的坐标变换

熟悉理论力学的读者应该能够领略到变分法在变换坐标系中的作用。比如，如果要将下面的平面二体问题方程
$$\left\{\begin{aligned}\frac{d^2 x}{dt^t}=\frac{-\mu x}{(x^2+y^2)^{3/2}}\\
\frac{d^2 y}{dt^t}=\frac{-\mu y}{(x^2+y^2)^{3/2}}\end{aligned}\right.\tag{1}$$
变换到极坐标系下，如果直接代入计算，将会是一道十分繁琐的计算题。但是，我们知道，上述方程只不过是作用量
$$S=\int \left[\frac{1}{2}\left(\dot{x}^2+\dot{y}^2\right)+\frac{\mu}{\sqrt{x^2+y^2}}\right]dt\tag{2}$$
变分之后的拉格朗日方程，那么我们就可以直接对作用量进行坐标变换。而由于作用量一般只涉及到了一阶导数，因此作用量的变换一般来说比较简单。比如，很容易写出，$(2)$在极坐标下的形式为
$$S=\int \left[\frac{1}{2}\left(\dot{r}^2+r^2\dot{\theta}^2\right)+\frac{\mu}{r}\right]dt\tag{3}$$
对$(3)$进行变分，得到的拉格朗日方程为
$$\left\{\begin{aligned}&\ddot{r}=r\dot{\theta}^2-\frac{\mu}{r^2}\\
&\frac{d}{dt}\left(r^2\dot{\theta}\right)=0\end{aligned}\right.\tag{4}$$
就这样完成了坐标系的变换。如果想直接代入$(1)$暴力计算，那么请参考《方程与宇宙》:二体问题的来来去去(一)

点击阅读全文...

实变函数中有一个勒贝格控制收敛定理，一般认为它是判断积分和取极限可交换的很好用的方法。勒贝格控制收敛定理是说，如果定义在集合$E$上的函数列$\left\{f_n(x)\right\}$满足$|f_n(x)|\leq F(x)$，而$F(x)$在$E$上可积，那么积分和取极限就可以交换，即
$$\lim_{n\to\infty}\left(\int_E f_n (x)dx\right)=\int_E \left(\lim_{n\to\infty}f_n (x)\right)dx$$
本文不打算谈该定理的证明，只是谈谈该定理的应用相关的话题。首先，请有兴趣的读者，做做以下题目：
$$\lim_{n\to\infty}\left(\int_0^1 \frac{n^2 x}{1+n^4 x^4}dx\right)$$

点击阅读全文...

斯特灵近似，或者称斯特灵公式，最开始是作为阶乘的近似提出
$$n!\sim \sqrt{2\pi n}\left(\frac{n}{e}\right)^n$$
符号$\sim$意味着
$$\lim_{n\to\infty}\frac{\sqrt{2\pi n}\left(\frac{n}{e}\right)^n}{n!}=1$$
将斯特灵公式进一步提高精度，就得到所谓的斯特灵级数
$$n!=\sqrt{2\pi n}\left(\frac{n}{e}\right)^n\left(1+\frac{1}{12n}+\frac{1}{288n^2}\dots\right)$$
很遗憾，这个是渐近级数。

相关资料有：
https://zh.wikipedia.org/zh-cn/斯特灵公式

https://en.wikipedia.org/wiki/Stirling%27s_approximation

本文将会谈到斯特灵公式及其渐近级数的一个改进的推导，并解释渐近级数为什么渐近。

点击阅读全文...