[SICP] 用流来解微分方程

win_hate

SICP 的 3.5.4 节以解微分方程为例，来解释如何使用“延迟参数”。抛开延迟参数不谈，解微分方程也是很有意思的。下面是本人对该节内容的一点评论。


ps: 本来这个帖子应该是今天早上发的。但不知为什么，在发送的时候 CU 说我只是个"游客"，要求我登录。登录后文字就全都没有了。郁闷中......

win_hate

给定微分方程 y'=f(y)，初始值 y(x_0)，求 y(x).

把区间 [x0, x] 平均分成 n 份，得到序列 x_0, x_1, ..., x_n=x。 如果能计算出相应的 y 在这些点上的导数值 y'(x_0), y'(x_1), ..., y'(x_n)，那么根据牛顿---莱布尼茨公式，我们有:

y(x_n) － y(x_0) ～ [y'(x_0)+y'(x_1)+...+y'(x_{n-1})]dx                              

所以 y(x_n) = y(x_0) + [y'(x_0)+y'(x_1)+...+y'(x_{n-1})]dx

但是，并不能直接从 x_i 得到 y'(x_i)，因为 y'(x)=f(y(x))。为求得 y'(x_i)，我们需要先求出 y(x_i)，而要想求出 y(x_i)，就需要先求出 y'(x_0), y'(x_2), ...., y'(x_{i-1})。 最终归结到初值 y(x_0) 上。

也就是说，y 和 y' 序列是互递归的，可以用如下序列来示意：

y(x_0)-> y'(x_0) -> y(x_1) -> y'(x_1) -> ...... -> y(x_n) -> y'(x_n) ......

显然，这可以看成一个反馈系统，所以用流来描述是非常自然的。具体的代码 SICP 上有，这里就不多说了。我们说点书上没有的。


从递推的观点看，y(x_n)=y(x_0)+[f(y(x_0))+f(y(x_1)+...+f(y(x_{n-1})]dx 可以表示为

y(x_n)=y(x_{n-1})+f(y(x_{n-1})dx

如果把 y(x_n) 记成 g(n)，则 g(n)=g(n-1)+f(g(n-1))dx。这是个简单明了的迭代关系，数学上称为 Euler 向前公式。它是最简单(陋)的求解一阶微分方程的数值方法，在其基础上可以衍生出更精确，但也更复杂的龙格---库塔(Runge-Kutta)方法。

SICP 上给出了如下的例子: y'=y, y(0)=1, 求 y(1).

我们知道这个 y 是 e^x，且 e^1=e.  用上面的公式进行计算对比是很有意思的。

由于 y'(x_{n-1})=y(x_{n-1})，有

y(x_n)=y(x_{n-1})+y(x_(n-1))dx，即： y(x_n)=y(x_{n-1}) (1+dx)

容易得到 y(x_n)=y(x_0) (1+dx)^n

如果令 x_0=0, x_n = 1, dx=0.001，有 n=1000。所以

e=y(1)=(1+0.001)^1000= 2.7169239322358924573830881219475771890

更一般地，有

y(1)=(1+1/n)^n，当 n->＋inf 时，有 y(1)=e。

这个极限在此时此地跳出来，一点也不意外。

[ 本帖最后由 win_hate 于 2009-3-13 18:26 编辑 ]