不动点迭代法--从一行代码干大事说起

highflybird · 发表于 2022-7-26 19:59

本帖最后由 highflybird 于 2022-7-26 20:16 编辑

一、计算器上的连按
cos(cos(cos(cos(cos(....)))))是多少？

我还记得以前在中学的时候，曾经对计算器着迷。有一天，发现了一个有趣的事情，对一个数进行余弦计算，不管你开始输入一个什么值，连按几次cos键后，就会发现数字再也不变了，(如果你设置的是弧度模式，可能要几十次，我真的这样测试过！）但当时对其中的缘故并未深究。直到很多年后，我才知道，这个就是不动点迭代。
二、不动点迭代
更有个定理描述这一类的迭代：巴拿赫不动点定理。具体的我也也不讲了，我也没那个水平讲清楚。
这个不动点迭代的最后得到的结果，其实就是解方程f(x)=x。对于余弦的这个迭代，就是求解非线性方程 cos(x)=x。这个迭代收敛比较慢。如果采用角度制的话，就是解方程cos(x*pi/180)=x，这个收敛很快，所以几次按键就可以得到很精确的解。
自从有了编程后，利用不动点迭代法，就可以快速解方程。
我给你展示如何用LISP代码求解 cos(x)=x，代码简洁到可以只用一行：
假设x我已经赋初始值1（随便什么值都可以），运行如下代码：

(repeat 100 (setq x (cos x)))

就这样，我们得到了LISP 允许范围的高精度解：

(rtos x 2 20) "0.7390851332151607"

从这里也显示了LISP代码是如何高效简洁，从这方面来说，它胜过了很多的编程语言。
在这里说明一点的是，实际迭代要不了100次，大概在81次就已经得到解了。
三、提高收敛速度
当然不同的函数，收敛的速度可能会天差地别。要提高收敛速度，还是需要储备一点数学知识的。
譬如对于这个求解cos(x)=x的，我们需要把两种方法结合在一起，就可以极大程度地改善收敛速度。
首先想到的自然是伟大的牛顿法：
把它稍微变化一下，等于是求解 f(x)=x-cos(x)的零点。求导得到，f'(x)=1+sin(x)
依据牛顿法思想，我们很容易得到另外一个不动点迭代：
g(x)=x-(x-cos(x))/(1+sin(x))
这次再也不需要100次的迭代了，实际需要5次就可以了，我这里多加1次吧。改进的代码如下：

(repeat 6 (setq x (- x (/ (- x (cos x)) (1+ (sin x))))))

我们再来看一个例子：
求解： f(x)=x^5-5*x-2=0.
一般一元五次方程无根式解。所以这里我们不能指望有什么求解公式了，只好采用数值求解了。
利用不动点迭代和牛顿法: g(x)=x-f(x)/f'(x) 即 g(x)=x-(x^5-5*x-2)/(5*x-5)
下面两行代码，如果把初始化撇开不算的话，真正的代码只有一行：

(setq x 0.5) ;先初始化x
(repeat 10 (setq x (- x (/ (- (* x x x x x) (* 5 x) 2) (- (* 5 x x x x) 5)))))

这样就得到了20位精度的解。
再次感叹LISP语言的伟大！
当然还有很多方法可以提高收敛速度。有很多工程数学或者数值计算的书都有提到不动点迭代法。读者如果要细究，请不妨参阅这些书籍。
我这里推荐一本：
《工程数学模型及数值计算方法》刘小华编写
网上的视频也有很多关于不动点的，譬如B站的这个：

四、完善
我上面的关于不动点迭代的程序虽然能达到目的，但是还是很粗糙，下面我对它稍微进行了一下改造：

;;;利用不动点迭代法求解非线性方程
(defun c:bdd (/ EPS EXPR I X X0 X1 X2)
(setq expr (getstring "\n输入字符串表达式:" ))
(CAL:Expr2Func expr 'func '(x))
(initget 1)
(setq x1 (getreal "\n下:"))
(initget 1)
(setq x2 (getreal "\n上:"))
(setq i 0)
(setq eps 1e-14)
(if (> x1 x2)
(setq x0 x2)
(setq x0 x1)
)
(setq x (func x0))
(while (and (< i 100) (> (abs (- x x0)) eps))
(setq x0 x)
(setq x (func x0))
(setq i (1+ i))
)
(princ "\n方程的解是: ")
(princ (rtos x0 2 20))
(princ "\t方程的值是: ")
(princ (rtos x 2 20))
(princ "\t迭代次数: ")
(princ i)
(princ)
)

此段程序需要加载我的表达式求值的子程序，包含在附件里面。

当然，读者还可以按照自己的要求进一步完善，利用它，基本能解决你遇到的大多数问题。
五、延伸
不动点迭代不仅可以对实函数迭代，还可以进行复函数迭代。还有其它形式的不动点迭代。

也许你可能会发现，对于某些函数，利用不动点迭代法，可能根本不收敛，为什么呢？可能要涉及到高深的数学知识了。
另外，不动点迭代可能不只是一个不动点，可以由一个不动点，跳到另一个不动点，对某一类迭代来说，可能有多个不动点。
譬如分形，典型的例子如Mandelbrot分形，Julia分形，都存在多个不动点。
这里面牵涉到的东西用数学分析起来，太复杂了，恐怕连现在的数学家也没完全弄清楚。
下面有一个视频介绍：

好了，暂时介绍到这里。本人水平有限，文中错误请大家多多指正。
那么，你还有一些一两行就能干成大事的代码呢？不妨分享给大家，谢谢！

tigcat · 发表于 2022-7-26 21:34

计算器要感谢牛顿和泰勒,是他们提供了理论支持.

自贡黄明儒 · 发表于 2022-7-27 07:11

大师数学很好，与lisp没有关系吧？

highflybird · 发表于 2022-7-27 08:12

自贡黄明儒发表于 2022-7-27 07:11
大师数学很好，与lisp没有关系吧？

与LISP没多大关系，只不过用这个简单而有效的工具，结合不动点理论，就可以变成解决实际问题的强有力手段。文末也期待各位网友分享看起来简单，实际能完成复杂工作的代码。

20060510412 · 发表于 2022-7-27 08:50

高飞鸟老师重出江湖，实在是幸甚。

ynhh · 发表于 2022-7-27 08:50

高飞

czb203 · 发表于 2022-7-27 08:58

高飞鸟老师重出江湖，实在是幸甚。

zgzzsn · 发表于 2022-7-29 06:13

敬佩高飞鸟老师

纵横八方 · 发表于 2022-7-29 20:24

本帖最后由纵横八方于 2022-7-30 14:42 编辑

鸟叔真是神奇！超级赞，利用鸟叔这个迭代，我写了个已知弦长弧长求半径飞快
(SETQ L (getdist "指定弧长") D (getdist "指定弦长"))
;(SETQ L 1026.45691 D 1000);举例
C:BDD
输入此表达式：
X-(L*SIN(X)-X*D)/(L*cos(X)-D)
圆心角马上出来 R=(/ L (*2 圆心角))就的出来了

大家测试一下
命令: BDD
输入字符串表达式:X-(L*SIN(X)-X*D)/(L*cos(x)-D)
下:1000
上:1
方程的解是: 0.3947911107667987 方程的值是: 0.3947911107667998 迭代次数: 7

mahuan1279 · 发表于 2022-8-1 11:02

本帖最后由 mahuan1279 于 2022-9-28 17:04 编辑

中国剩余定理。正整数a,b,c两两互质（a<b<c），一正整数N除以a余x，除以b余y，除以c余z，求N最小为多少？
_$ (defun  f(a b c x y z)
(defun f1 (M n)
(setq sum M)
(while (/= (setq i (rem sum n)) 1)
   (setq sum (+ sum M))
   )
(/ sum M)
)
(setq N (rem  (+ (* x b c (f1 (* b c) a))  (* y a c (f1 (* a c) b)) (* z a b (f1 (* a b) c)))  (* a b c))))
F
_$ (f  3 5 7 2 4 5)
89
_$

账号		自动登录	找回密码
密码			注册

[【高飞鸟】] 不动点迭代法--从一行代码干大事说起

本帖子中包含更多资源

评分

评分