当Matlab遇上牛顿法
By 苏剑林 | 2013-05-22 | 58468位读者 |牛顿法是求方程近似根的一个相当有用而且快捷的方法,我们最近科学计算软件课程(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://spaces.ac.cn/archives/1995
更详细的转载事宜请参考:《科学空间FAQ》
如果您还有什么疑惑或建议,欢迎在下方评论区继续讨论。
如果您觉得本文还不错,欢迎分享/打赏本文。打赏并非要从中获得收益,而是希望知道科学空间获得了多少读者的真心关注。当然,如果你无视它,也不会影响你的阅读。再次表示欢迎和感谢!
如果您需要引用本文,请参考:
苏剑林. (May. 22, 2013). 《当Matlab遇上牛顿法 》[Blog post]. Retrieved from https://spaces.ac.cn/archives/1995
@online{kexuefm-1995,
title={当Matlab遇上牛顿法},
author={苏剑林},
year={2013},
month={May},
url={\url{https://spaces.ac.cn/archives/1995}},
}
May 23rd, 2013
这是个好方法,我以前也困惑过。
不过以后我要装matlab2012来跟上你的代码啦
这不是2012版专有的函数,我在上课时的2011版也可以实现。
June 1st, 2013
我发现过一个类似牛顿法的迭代法,没做什么深入研究,楼主有兴趣吗?
September 7th, 2013
要不试试mathematica【哪有不合适关键词啊】?可以求导求积分。
可以,但是数值编程还是matlab比较简单,尤其是对于数组的处理,以及代码上的简洁性~
July 17th, 2015
听起来不错
不懂matlab,下半年大二开始学.
May 3rd, 2018
我系省外读紧书嘅广州人,写得几好啊。多谢晒~