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(随便什么值都可以),运行如下代码:
(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 ( (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分形,都存在多个不动点。
这里面牵涉到的东西用数学分析起来,太复杂了,恐怕连现在的数学家也没完全弄清楚。
下面有一个视频介绍:
好了,暂时介绍到这里。本人水平有限,文中错误请大家多多指正。
那么,你还有一些一两行就能干成大事的代码呢?不妨分享给大家,谢谢!
**** Hidden Message *****

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最小为多少?
_$ (defunf(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
_$ (f3 5 7 2 4 5)
89
_$
页: [1]
查看完整版本: 不动点迭代法--从一行代码干大事说起