免费注册 查看新帖 |

Chinaunix

  平台 论坛 博客 文库
最近访问板块 发新帖
查看: 3798 | 回复: 1

[SICP] 用流来解微分方程 [复制链接]

论坛徽章:
0
发表于 2009-03-13 18:22 |显示全部楼层
SICP 的 3.5.4 节以解微分方程为例,来解释如何使用“延迟参数”。抛开延迟参数不谈,解微分方程也是很有意思的。下面是本人对该节内容的一点评论。


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

论坛徽章:
0
发表于 2009-03-13 18:22 |显示全部楼层
给定微分方程 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 编辑 ]
您需要登录后才可以回帖 登录 | 注册

本版积分规则 发表回复

  

北京盛拓优讯信息技术有限公司. 版权所有 京ICP备16024965号-6 北京市公安局海淀分局网监中心备案编号:11010802020122 niuxiaotong@pcpop.com 17352615567
未成年举报专区
中国互联网协会会员  联系我们:huangweiwei@itpub.net
感谢所有关心和支持过ChinaUnix的朋友们 转载本站内容请注明原作者名及出处

清除 Cookies - ChinaUnix - Archiver - WAP - TOP