Как вычислить функцию Brjuno в Maxima CAS?

#function #maxima #numerical-computing

Вопрос:

Я пытаюсь вычислить реальную функцию 1-Брюно для действительных чисел в x в диапазоне [0, 1]

   G(x):= block(
            [g],
        if (x=0) 
            then g : 0
            else g : float(1/x - floor(1/x)),
        return( g)
        )$
    
   bet(j,x) := block(
        [r],
        if (j=-1) 
            then r:1
            else r : product (G(x)^i, i, 0, j),
        return( r)
        )$
    
    A(i,x) := bet(i-1,x)*log(1/G(x)^i)$
    B(x):= sum(A(i,x), i, 0,10)$

Есть численные ошибки :

                         0
    expt: undefined: 0.0

вызвано функцией для некоторых значений x : 1/2,1/5,1/10,1/3,1/9,1/4,0

Как я могу решить эту проблему ?

1. Попробуй simp:false; tellsimp(0^0, 1); tellsimp(0.0^0, 1); simp:true; .

2. @RobertDodier expt: не определено: от 0 до отрицательного показателя.

3. Хорошо, это означает, что в какой-то момент у вас есть что-то/0. Есть ли ограничение, которое применяется в этом случае? или какой-то другой особый случай? Возможно, вам потребуется выяснить, является ли «что-то» само по себе нулем или ненулевым. Кстати, я попробовал несколько примеров, и кажется, что «от 0 до отрицательного показателя» встречается для «хороших» значений, таких как 1/2, но для случайных значений я не получаю ошибку.

4. @RobertDodier Я использую просто plot2d(B(x),[x,0,1]); тогда ошибок нет, но диаграмма выглядит иначе, чем в статье. Поэтому я перешел на дискретный набор точек ( рациональные углы) в диапазоне [0,1] ,

Ответ №1:

Это создает сюжет, похожий на бумагу. Я использовал уравнение из википедии и аппроксимацию непрерывной дробью:

 cflength: 10 $
B(x):= if length(x)=1 then 0 else (local(x), x[1]: 0, -log(cf2num(x)) cf2num(x)*B(rest(x))) $
cf2num(e):=ev(cfdisrep(e), numer, infeval) $
x0: makelist(i*sqrt(2)/1421, i, 1, 1000), numer $
x: map(cf, x0) $
y: map(B, x) $
draw2d(points(x0, y)) $

1. Идеальный. Я хотел бы поместить изображение в общий доступ. Это ваш код, так что йопу может это сделать. Если у вас нет времени, я могу сделать это с информацией об авторе. Дайте мне знать

2. @Adam Не стесняйтесь публиковать это. Спасибо вам за то, что вы это сделали.

3. @slitvinv Готово : en.wikipedia.org/wiki/Brjuno_number#Brjuno_function

4. commons.wikimedia.org/wiki/File:Brjuno_function.png

5. Я немного изменил код. Результат хороший, так что я надеюсь, что все в порядке

Ответ №2:

Функция Брюно имеет «логарифмические особенности в рациональных точках», поэтому с помощью простого »

 G(x):= block(

    [g],
    if (x=0) 
        then g : 0
        else g : float(1/x - floor(1/x)),
    return( g)
    )$
    

beta(j,x) := block(

    [r],
    if (j=-1) 
        then r:1
        else r : product (G(x)^i, i, 0, j),
    return( r)
    )$
    
    
A(i,x) := beta(i-1,x)*log(1/G(x)^i)$
B(x):= sum(A(i,x), i, 0,10)$
plot2d(B(x),[x,0,1], [y, -1, 10]);

дает сюжет :

Но диаграмма не совпадает с рисунком 1 ( слева ) в документе.