牛顿法解方程-使用浏览器(Geogebra可选)做实验

提起牛顿法解方程,大一以上的同学都听说过。但是如果要求做实验写代码,不少同学可能会心生畏惧,会不会很难很麻烦。事实上,核心代码只有一行而已。本文演示的实验只需要浏览器,可选Geogebra。

1.计算机问题求解的三类方法

裘宗燕老师把计算机问题求解的方法归结为三类,https://zhuanlan.zhihu.com/p/58943925

。例如就解二元一次方程

a * x * x + b * x - c = 0

而言,求解的方法可以从以下角度分为三类。

第一类方法是公式,

x1=(-b +(b^2-4*a*c)^1/2)/ (2*a)
x1=(-b -(b^2-4*a*c)^1/2)/ (2*a)

(在解存在的情况下,)我们把a,b,c代入公式就能求出方程的根。这种方法适用于 我们对于方程的规律全部都了解。对于计算机学生而言,这个公式相当于领域知识或者先验知识。

如果我们对方程的全局知识并不了解,比如五次及五次以上的方程已证明没有求根公式,那么就需要第二类方法,计算机的核心思路——搜索。我们可以想办法穷举/枚举解空间中的所有可能性,把这些x逐一代入方程a * x * x + b * x – c。如果哪个x代入方程以后得到的结果是0,那么当前的x就是一个解。

搜索问题需要考虑以下问题。1.是否可以穷举。有理数集可穷举,实数集不可穷举。二次方程有些解并非有理数,所以严格地说,第二类方法搜索不能得到方程的根。不过,在工程和技术领域,我们可以不那么严格,并不一定得到精确的实数解,而是可以允许在一定误差内(比如,方程代入x以后求得的值与0间的差的绝对值)的就算命中。2.如何穷举,按什么顺序枚举。从小大到?3.如何以更快的速度穷举。二分查找?

第二类方法搜索,需要了解问题的局部规律。方法的应用正是对局数规律在全局上的推广。

第三类方法甚至连局数规律也不必知道,只需要了解个别(也许很多)问题的实例和解的实例。你可能已经知道了,接下来是采用统计学的或神经网络或深度学习或机器学习的手段,把个例推广到一般(并且在一定范围内检验)。

牛顿法是第二类搜索方法,基于……一系列原理,简而言之,对牛顿的信仰,搜索速度飞快,可以采用。

2. 牛顿法 算法步骤的几何解释

假设我们要解的方程是y=(fx),要求误差小于0.01。

抛物线在平面直角坐标系里是这样的。

image2假设我们初始猜测的值是10。

此时误差不小于0.0.1。

过 (10, f(x))点做切线,红色的。

image3

以 切线与横轴交点 作为新的x1,过 (x1, f(x1))点做切线,绿色的。

image4

此时误差不小于0.0.1。

放大一下,是这样的。

以 切线与横轴交点 作为新的x2,过 (x2, f(x2))点做切线,蓝色的。

此时误差不小于0.0.1。

放大一下,是这样的。

image5

以 切线与横轴交点 作为新的x2,过 (x2, f(x2))点做切线,棕色的。

此时误差不小于0.0.1。

放大一下,是这样的。

image6

以 切线与横轴交点 作为新的x2,过 (x2, f(x2))点做切线,青色的。看起来比青色暗一些,是因为在这个分辨下与原方程的黑线重叠了。

此时误差不小于0.0.1。

放大一下,是这样的。

image7

青色的线与方程黑色的线几乎重叠,放大一下。

黑色的线是方程,青色的线是最后一次的切线。

此时的误差为F1E1=0.0014,小于0.01。搜索到了符合要求的解。

image8

总结,每个步骤都有一些共同的特征。

步骤1 猜测x_k,得到(x_k, f(x_k));

步骤2 过点(x_k, f(x_k))做方程的切线,交于(x_k+1, 0),这样得到x_k+1;

步骤3 检测f(x_k+1) 与 f(0) 间的差是否符合要求。

 

缩小看全局,可以看到 从T1(10, f(x))开始,切线 红-绿-蓝-棕-青 的倾斜过程,以及切线与横轴的交点迅速逼近方程的根。

image9

3. 递归公式推导

 

因为每个步骤有共同特征,因此可以抽象为递归或迭代过程。

我们根据几何解释求出递归公式。

图示为几何解释,其中(x_k, f(x_k))的高度为 f(x_k)。

image10

我们根据导数定义得到公式,其中f(x)和f'(x)根据方程由程序员手动求出。

公式在此!

image11

如果不会求导,可以借助符号运算工具,比如 Mathmatica,或者https://www.wolframalpha.com/。如果像我一样访问不到,再如 https://www.derivative-calculator.net/https://zh.numberempire.com/derivativecalculator.php

image12

或者

image13

3. JavaScript代码

为什么想起来用 JavaScrtip实现呢。因为正在重读 SICP(计算机程序的构造和解释,JavasScrtip版本),发现这个开发环境真是友好便利。本文即重现SICP中的例题,再次推荐SICP。

在浏览器里按F12,然后Ctrl+o(相信用菜单打开文件也一样)打开 newton.js 文件。还可以不编译调试,随手改随手运行。

代码很短,核心部分只有一行。以下分段解释代码。

// In Geogebra
// f: y=2 x^(2)+3 x-4
// NSolve(0=2 x^(2)+3 x-4)
// = {x = -2.35, x = 0.85}
a = 2;
b = 3;
c = -4;
epsilon = 0.01;

以上,要解的方程是y=2 x^(2)+3 x-4。为增加通用性(和趣味性),a,b,c三个参数可以修改。Epsilon是误差的最大值。

function f(x)
{
    return a*x*x+b*x+c;
}
function df(x) // derived function
{
    return 2*a*x+b;
}

方程,以及方程的切线。其中f(x)是方程,df(x)是方程的切线/导数。

以下就是核心部分了。

function guess(x)
{
    return Math.abs(f(x)) <= epsilon ? x : guess( x - f(x)/df(x));
}

这一行就是核心

Math.abs(f(x)) <= epsilon ? x : guess( x - f(x)/df(x) )

如果误差小于预期,那么当前的x就是方程的根;

如果误差大于预期,那么继续猜下去,猜测的x_k+1<= x_k – f(x)/df(x),即上面 递归公式推导中的“公式在此!” 手写公式中的f'(x)就是代码中的df(x)。

guess(0); //0.8517236753856471
guess (-1); //-2.350915564067371

这两行是测试运行,分别由左边和右边逼近得到两个不同的根。

4.验证

求解的根代入方程,可以确认f(x)误差小于预期。

也可以通过Geogebra算一下。在左边输入,求出方程的近似解。

image14

4.递归有多快

程序员会思前想后很多因素,也许在某个具体的问题中根本用不上,但是万一墨菲法则发作了呢。比如,牛顿法收敛有多快,递归需要多少次,会不会暴栈。

道理上,收敛非常快,非常~是多快。我们看一下感性认识。修改上述js代码,把guess改成下面这样。

function guess(x)
{  
    console.debug(x);
    return Math.abs(f(x)) <= epsilon ? x : guess( x - f(x)/df(x)
);
}

加了一行调试,把浏览器按以后F12的debug设置为允许显示。

image15

猜测一下,以10作为初始值,guess执行6次。

image16

如果猜测得糟糕一些呢,比如10000?guess执行16次。

image17猜测的初始值为原来的1000倍,guess增加了10次。

收敛有多快呢,我看把x系列放到excel里看一下。

无标题

猜得再离谱一些,1000000000。33次。33次这种量级,对于计算机而言,就是瞬间。

33

而且我们观察上面两图可以看出,最初的收敛非常陡峭,也就是说,极快地逼近目标值。如果容忍的误差大一些,会有更好地表现。

当然,牛顿法有适用范围,或者说有不适用的范围。然后,没有一种工程方法是完美的,只有不完美的工程师把某种工程方法应用在不适合的问题上。所以,牛顿法(以及Geogebra,SICP,JavaScriipt,Excel等等)真是强大啊。

微信图片_20220814182439

Leave a Reply

Your email address will not be published. Required fields are marked *