0

I am writing a simple newton method 

x_(n+1) = x_n - f(x_n) / f_prime(x_n)

to find the roots (can be a real number or a complex number) of a quadratic function: 

f(x) = a*x*x + b*x + c

(a, b, c are given constants and are all real numbers). I know Newton method will fail if the start point or some iteration point in the loop has a zero derivative. I want to use a if statement inside my for/while loop to avoid this situation. Does Julia have something like stop 0 syntax in Fortran ?

The generic Newton's Method root-finding code:

function newton_root_finding(f, f_diff, x0, rtol=1e-8, atol=1e-8)
    f_x0 = f(x0)
    f_diff_x0 = f_diff(x0)
    x1 = x0 - f_x0 / f_diff_x0
    f_diff_x1 = f_diff(x1)

    @assert abs(f_diff_x0) > atol + rtol * abs(f_diff_x0) "Zero derivative. No solution found."

    while abs(f_x0) > atol + rtol * (abs(f_x0))
        x0 = x1
        f_x0 = f(x0)
        f_diff_x0 = f_diff(x0)
        x1 = x0 - f_x0 / f_diff_x0
    end

    return x1
end 

function quadratic_func(x)
    a = 1.0
    b = 0.0 
    c = 2.0
    return a*x*x + b*x + c
end

function quadratic_func_diff(x)
    a = 1.0
    b = 0.0 
    c = 2.0   
    return 2.0*a*x + 1.0*b + 0.0*c

end

newton_root_finding(quadratic_func, quadratic_func_diff, 1.0 + 0.5im)

In the above code I used a @assert macro to make it happens, but I don't want to use any macro. I want to use a if statement inside my while loop to halt it. Another thing I've noticed is that if I change to @assert abs(f_diff_x0) != 0 this test will be ignored. Is that because of some round-off errors that "zero derivative" doesn't exactly equal to 0?

phipsgabler
  • 20,535
  • 4
  • 40
  • 60
Tack_Tau
  • 628
  • 1
  • 6
  • 12
  • 2
    using @assert is exactly like using an if statement: ``` julia> @macroexpand @assert True == False "Wrong" :(if True == False nothing else (Base.throw)((Base.AssertionError)("Wrong")) end) ``` – Thomas Jalabert Oct 24 '19 at 13:53

2 Answers2

2

To address your first question, you can use break to exit the while loop, like

function test()
    i = 0
    while true
        i += 1
        if i > 10
            break
        end
    end
    return i
end

As to your second question, when comparing floating point numbers it is often better to use isapprox (provide an atol if you compare against zero) instead of == or !=.

Lufu
  • 436
  • 4
  • 11
2

The way to exit from the inside of a loop in general is a break statement; a return fulfills the same purpose, because it just exits the whole function.

For the comparisons you can use Base.isapprox(x, y; atol=atol, rtol=rtol). It's documentation starts with:

Inexact equality comparison: true if norm(x-y) <= max(atol, rtol*max(norm(x), norm(y))).

norm falls back to abs for numbers. And I think you might have a bug in both comparisons, always comparing the value at x0 to itself.

As for the breaking on zero derivatives, an @assert is, I think, appropriate here: if you get zero derivative, you don't stop iteration and return the result, but you throw an error to signify an infeasible condition. I'd thus write your function as follows:

function newton_root_finding(f, ∂f, x0, rtol=1e-8, atol=1e-8)
    x_old = x0
    y_old = f(x0)

    while true
        df_old = ∂f(x_old)
        @assert !isapprox(df_old, 0, rtol=rtol, atol=atol) "Zero derivative. No solution found."

        x_new = x_old - y_old / df_old
        y_new = f(x_new)

        isapprox(y_old, y_new, rtol=rtol, atol=atol) && return x_new

        x_old, y_old = x_new, y_new
    end
end

This returns 3.357392012620626e-26 + 1.4142135623730951im on your test case, approximately sqrt(2)im.

phipsgabler
  • 20,535
  • 4
  • 40
  • 60