Как остановить цикл в Julia и одновременно распечатать ErrorMsg без использования макросов?

Question

Я пишу простой метод Ньютона

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

, чтобы найти корни (могут быть действительным или комплексным числом) квадратичной функции:

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

(a, b, c даны постоянные и являются действительными числами). Я знаю, что метод Ньютона потерпит неудачу, если начальная точка или некоторая итерационная точка в цикле имеет нулевую производную. Я хочу использовать оператор if внутри моего цикла for / while, чтобы избежать этой ситуации. Есть ли у Джулии что-то вроде stop 0 синтаксиса в Фортране?

Общий код поиска корня метода Ньютона:

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)

В приведенном выше коде я использовал макрос @assert, чтобы сделать егобывает, но я не хочу использовать какой-либо макрос. Я хочу использовать оператор if внутри моего цикла while, чтобы остановить его. Еще одна вещь, которую я заметил, это то, что если я перейду на @assert abs(f_diff_x0) != 0, этот тест будет проигнорирован. Это из-за некоторых ошибок округления, что «нулевая производная» точно не равна 0?

Answer 1

Чтобы ответить на ваш первый вопрос, вы можете использовать break для выхода из цикла while, например

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

Что касается вашего второго вопроса, при сравнении чисел с плавающей запятой часто лучше использовать isapprox (укажите atol, если сравнивать с нулем) вместо == или !=.

Answer 2

Способ выхода из внутренней части цикла в общем случае - оператор break;return выполняет ту же цель, потому что он просто выходит из всей функции.

Для сравнения вы можете использовать Base.isapprox(x, y; atol=atol, rtol=rtol). Документация начинается с:

Неточное сравнение на равенство: true, если norm(x-y) <= max(atol, rtol*max(norm(x), norm(y))).

norm возвращается к abs для чисел. И я думаю, что вы могли бы иметь ошибку в обоих сравнениях, всегда сравнивая значение в x0 с самим собой.

Что касается разрыва по нулевым производным, я думаю, здесь @assert уместно: есливы получаете нулевую производную, вы не останавливаете итерацию и не возвращаете результат, но вы выдаваете ошибку, чтобы обозначить недопустимое условие. Таким образом, я бы написал вашу функцию следующим образом:

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 

Это возвращает 3.357392012620626e-26 + 1.4142135623730951im в вашем тестовом примере, приблизительно sqrt(2)im.