When Matlab Meets Newton's Method

Newton’s method is a quite useful and fast method for finding approximate roots of equations. A recent assignment for our Scientific Computing Software course (Matlab) was to write a program to find approximate solutions to equations, which involved Newton’s method. Our goal was to achieve the following: the user inputs an equation, and the script automatically finds the root. This seems like a fairly simple iterative loop program, but due to the specific nature of Matlab itself, it presented several difficulties.

Matlab was born for numerical computation (especially matrix operations); therefore, it is not particularly adept at handling symbolic computation. This poses difficulties for our programming. A quick search online reveals that most Matlab Newton’s method programs online require the user to input both the equation and its derivative function simultaneously. This is clearly inconvenient because Matlab itself possesses differentiation capabilities. Let us analyze where the difficulties lie.

The most basic functionality we want to implement is defining a function, then being able to calculate specific function values based on that function, automatically finding the derivative of that function, and then calculating the derivative values. These seemingly basic functions are difficult to reconcile in Matlab because the definition of a "function" in Matlab is very broad. An M-file with a specific function is called a "function," and an operational expression f(x) might also be a function. Clearly, the latter can be differentiated, while the former cannot. Consequently, Matlab takes a "one-size-fits-all" approach—it cannot differentiate a function handle directly!!

So, how does Matlab’s differentiation function work? In fact, it performs string operations on functional expressions. For example:

diff(x^2)

In this case, the content inside the parentheses is a symbolic expression, not a function handle. What do functions in Matlab look like? A simple function definition in Matlab can be written as:

f=@(x) x^2

or

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

This way, entering f(2) yields 2^2. Functions defined like this cannot be differentiated directly. So how do we reconcile the contradiction between the two? There is a little trick:

diff(f(x))

Note that here you cannot change f(x) to f; f is a function handle, while f(x) is the value of that function when the independent variable is x! This point is quite important. Of course, the result of such differentiation is also a symbolic expression; we need to redefine it as a function:

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

With this, we can easily write the following script:

function f=jfc(g)

syms x
p=0.000001;  % Define precision
n=0;  % Calculate number of iterations

s=inline(g,'x');  % Define the function
ds=inline(diff(s(x)),'x');   % Define the derivative function

x=1;  % Initial value
e=s(x); 

% Iteration process below
while abs(e)>p
    n=n+1;
    x=x-s(x)/ds(x);
    e=s(x);
end

% Output
[x n e]

This is a simple program and is not yet perfect (for instance, initial values and other conditions have not been handled well). Usage:

jfc(x^2-2)

We have another approach, which is to keep the entire process purely in symbolic form without defining any functions. However, in this case, we cannot evaluate the function normally. How do we evaluate it? There is a clever trick—the "substitution function" in Matlab. For example, if f(x)=x^2, to find f(5):

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

This will output f=25. subs is a magical function; it is equivalent to replacing x in x^2 with 5. This is a method for evaluating functions! Thus, 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

After reading this article, readers who are familiar with Matlab might feel that this is not a particularly remarkable trick. Indeed, it seems very simple now, but I also wonder why similar programs do not appear online. I dare not jump to conclusions. Perhaps it is because I have just started, and I am proceeding from the most basic principles. Of course, there is another small reason: when the built-in functions of a program are not enough, you write new ones yourself; when you cannot adapt to the world, let the world adapt to you.

When reposting, please include the original address of this article: https://kexue.fm/archives/1995

For more detailed reposting matters, please refer to: Scientific Space FAQ