When Matlab Meets Newton's Method

当Matlab遇上牛顿法 苏剑林 · 2013-05-22 · Mathematics equation software approximation

Newton's method is a rather useful and fast way of finding approximate roots of an equation. A recent assignment in our Scientific Computing Software course (Matlab) was to write a program for finding approximate solutions of equations, which involved Newton's method. The goal we were asked to achieve was: the user types in an equation, and the script automatically finds the root. This looks like a fairly simple looping/iterative program, but because of some quirks specific to Matlab, it turned out to involve quite a few difficulties.

Matlab was built for numerical computation (especially matrix operations), so it isn't particularly good at symbolic computation. This is exactly what makes things difficult for us here. A quick search online shows that all the Matlab Newton's-method programs out there require the user to supply both the equation and its derivative — which is obviously inconvenient, since Matlab itself already has a differentiation feature. Let's look at where the difficulty actually lies.

The most basic functionality we need to implement is: define a function, then be able to evaluate that function at specific points, and also automatically compute its derivative and evaluate the derivative. These seemingly basic requirements turn out to be hard to reconcile in Matlab, because Matlab's notion of "function" is very broad — an M-file implementing some specific piece of functionality is called a "function", but an expression like $f(x)$ can also be called a "function". Clearly the latter can be differentiated, while the former cannot. So Matlab draws a hard line here — you cannot differentiate a function! more

So how does Matlab's differentiation feature actually work? In fact, it operates on the string form of an expression, for example

diff(x^2) 

where what's inside the parentheses is a string expression, not a function. So how are functions defined in Matlab? A simple way to define a function in Matlab is to write

f=@(x) x^2 

or

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

and then entering f(2) gives you $2^2$. A function defined this way cannot be differentiated. So how do we reconcile the two? There's a small trick:

diff(f(x)) 

Note that here you cannot change f(x) to just f — f is a function, while f(x) is the value of that function when the independent variable is x! This distinction matters a great deal. Of course, the result of differentiating this way is also a string, so we need to redefine it as a function:

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

With that, we can easily write the following script:

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]

This is a simple program, and it's not fully polished yet (for instance, the handling of the initial value and other edge cases still needs work).

Usage:

jfc(x^2-2)

There's also another approach: carry out the entire process purely in terms of strings, without defining any functions at all. But then how do we evaluate things numerically? There's a neat trick for this — Matlab's "substitution function". For example, given $f(x)=x^2$, to find $f(2)$:

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

This will output f=25. Subs is a wonderfully handy function: it's equivalent to substituting x=2 into $x^2$, which is exactly how you evaluate a function! So we can also write the following version of the program:

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]

Conclusion

Readers who are already fairly familiar with Matlab will, after reading this, probably feel that this isn't really much of a trick. Indeed, looking back now it really is quite simple, and yet I can't quite figure out why nothing like this seems to appear anywhere online. I won't jump to conclusions — maybe it's just because I'm still a beginner, starting from the most basic principles. But there's also a small underlying reason: when the ready-made functionality of a program isn't enough, you write new functionality yourself; when you can't adapt to the world, make the world adapt to you.