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

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

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

copy

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

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

copy

(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次吧。改进的代码如下：

copy

(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)
下面两行代码，如果把初始化撇开不算的话，真正的代码只有一行：

copy

(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站的这个：
四、完善
我上面的关于不动点迭代的程序虽然能达到目的，但是还是很粗糙，下面我对它稍微进行了一下改造：

copy

;;;利用不动点迭代法求解非线性方程
(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 ( (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:00

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

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

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

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

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

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

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

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

高飞

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

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

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

敬佩高飞鸟老师

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

鸟叔真是神奇！超级赞，利用鸟叔这个迭代，我写了个已知弦长弧长求半径飞快
(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:00

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

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

发表回复

精选推荐