牛顿法是求方程近似根的一个相当有用而且快捷的方法,我们最近科学计算软件课程(Matlab)的一个作业就是编写求方程近似解的程序,其中涉及到牛顿法。我们要实现的目标是,用户输入一道方程,脚本就自动求出根来。这看起来是一个挺简单的循环迭代程序,但是由于Matlab本身的特殊性,却产生了不少困难。

Matlab是为了数值计算(尤其是矩阵运算)而生的,因此它并不擅长处理符号计算。这就给我们编程带来了困难。在网上随便一搜,就可以发现,网上的Matlab牛顿法程序都是要求用户同时输入方程及其导函数,这显然是不方便的,因为Matlab本身就具备了求导功能。下面我们来分析一下困难在哪里。

我们要实现的最基本功能是定义一个函数,然后可以根据该函数求具体的函数值,并且自动求该函数的导数,接着求导数值。这些看起来很基本的功能在Matlab中却很难调和,因为Matlab的“函数”定义很广,一个具有特定功能的M文件叫“函数”,一个运算式$f(x)$也可能是一个函数,显然后者是可以求导的,前者却不行,所以Matlab一刀砍——不能对函数求导!!

那Matlab的求导功能是怎样的呢?事实上,它是对函数式的字符串操作的,比如

diff(x^2) 

其中括号里边是一个字符串式子,而不是函数。那Matlab中的函数是怎样的呢?Matlab中简单定义一个函数,可以写

f=@(x) x^2 

f=inline(x^2,'x')

这样输入f(2)求出$2^2$了。这样定义的函数是不能求导的。那么怎么调和两者的矛盾呢?里边有个小技巧:

diff(f(x)) 

注意这里不能将f(x)改为f,f是一个函数,而f(x)是该函数自变量为x时的值!这一点相当重要。当然,这样求导的结果也是一个字符串,我们要把它重新定义为函数:

df=inline(diff(f(x)),'x') 

这样我们就可以轻松写出下面的脚本:

function f=jfc(g)

syms x
p=0.000001;  %定义精度
n=0;  %计算迭代次数

s=inline(g,'x');  %定义函数
ds=inline(diff(s(x)),'x');   %定义导函数

x=1;  %初始值
e=s(x); 

%下面是迭代过程
while abs(e)>p
    n=n+1;
    x=x-s(x)/ds(x);
    e=s(x);
end

%输出
[x n e]

这是一个简单的程序,还不完善(比如初值等等情况还没有处理好)。
用法:

jfc(x^2-2)

我们还有另外一个思路,就是整个过程纯粹用字符串形式,而不定义任意函数。但是这样我们就不能根据函数求值了,怎么求值呢?有一个巧妙的技巧——Matlab中的“替换函数”。比如$f(x)=x^2$,求$f(2)$:

f=subs(x^2,'5',x)

这将输出f=25。Subs是一个神奇的函数,它等价于将$x^2$中的x替换成2,这是求函数值的方法!于是我们也可写出以下版本的程序:

function f=jfc(g)

p=0.000001;
x=1;
n=0;

ds=diff(s);
e=subs(s,'x',x);

while abs(e)>p
    n=n+1;
    ds0=subs(ds,'x',x);
    x=x-e/ds0;
    e=subs(s,'x',x);
end

[x n p]

结论

看完本文,比较了解Matlab的读者会感觉这是一个并没有什么了不起的技巧。的确,现在看来这真的很简单,可是我也想不到为什么网上并没有出现类似的程序。我不敢妄下定论。也许是因为我才刚刚入门,一些从最基本的原理出发吧。当然,还有一个小原因,就是当程序现成的功能不够用时,就自己动手编写新的;当你不能适应世界时,就让世界适应你吧。

转载到请包括本文地址:https://kexue.fm/archives/1995

更详细的转载事宜请参考:《科学空间FAQ》

如果您还有什么疑惑或建议,欢迎在下方评论区继续讨论。

如果您觉得本文还不错,欢迎分享/打赏本文。打赏并非要从中获得收益,而是希望知道科学空间获得了多少读者的真心关注。当然,如果你无视它,也不会影响你的阅读。再次表示欢迎和感谢!

如果您需要引用本文,请参考:

苏剑林. (May. 22, 2013). 《当Matlab遇上牛顿法 》[Blog post]. Retrieved from https://kexue.fm/archives/1995

@online{kexuefm-1995,
        title={当Matlab遇上牛顿法},
        author={苏剑林},
        year={2013},
        month={May},
        url={\url{https://kexue.fm/archives/1995}},
}