Matlab: символьное вычисление членов посл-ти A006206

[maximkarimov]

GAA в сообщении #1458183 писал(а):

maximkarimov, можно проще в MuPAD
Matlab R2013b::MuPAD

A006206 := proc(n)
  local sum, d;
 begin
     sum := 0;
     for d in numlib::divisors(n) do 
       sum := sum + numlib::moebius(n/d)*(numlib::fibonacci(d+1)+numlib::fibonacci(d-1))
     end;
    sum/n;
 end;

A006206(89)
    44720842390302450 

Суперрр! А можно ли это превратить в m-функцию (чтоб ее можно было в скрипт встроить и там вызывать)?

[wrest]

maximkarimov
Калькулятор PARI/GP избавлен от потерь точности при вычислениях с рациональными (и соответственно целыми) числами. Гляньте, на досуге. У нас и тема есть: «интерактивный курс: введение в программирование на PARI/GP»
Код для A006206 (взят из статьи oeis) очень короткий:
a(n)=if(n<1, 0, sumdiv(n, d, moebius(n/d)*(fibonacci(d-1)+fibonacci(d+1)))/n)

[ShMaxG]

maximkarimov
Если без MuPad, можете попробовать использовать функцию vpa.

[GAA]

maximkarimov, я просто код для Maple из A006206 подправил под MuPad. У меня есть возможность пробовать на достаточно старой версии Matlab. В R2013b вроде бы получается длиннее (нужно ещё символьную m-функцию для fibonacci сделать), чем у Вас.

maximkarimov в сообщении #1458153 писал(а):

Проблема в том, что все вместе не получается - ошибка выскакивает на самой последней строчке:

Используется синтаксис Matlab M

a= s/n;
 

Крутил и так и эдак - не получается. Потому и прошу помощи.

Вам нужно символьное деление. При n=89 сумма s=980154972736918050

>> sym('3980154972736918050')/sym('89')
ans = 44720842390302450 

Т.е. если s и n — символьные, то в указанном Вами месте ошибки не должно быть.
Приведите, пожалуйста, полностью текст m-функций (с использованием символьных вычислений) и сообщение об ошибке.

[maximkarimov]

GAA в сообщении #1458201 писал(а):

maximkarimov, я просто код для Maple из A006206 подправил под MuPad. У меня есть возможность пробовать на достаточно старой версии Matlab. В R2013b вроде бы получается длиннее (нужно ещё символьную m-функцию для fibonacci сделать), чем у Вас.

maximkarimov в сообщении #1458153 писал(а):

Проблема в том, что все вместе не получается - ошибка выскакивает на самой последней строчке:

Используется синтаксис Matlab M

a= s/n;
 

Крутил и так и эдак - не получается. Потому и прошу помощи.

Вам нужно символьное деление. При n=89 сумма s=980154972736918050

>> sym('3980154972736918050')/sym('89')
ans = 44720842390302450 

Т.е. если s и n — символьные, то в указанном Вами месте ошибки не должно быть.
Приведите, пожалуйста, полностью текст m-функций (с использованием символьных вычислений) и сообщение об ошибке.

Вот код:

Используется синтаксис Matlab M

n=89;
n=sym(n);
s=0;
for i=1:n
    if mod(n,i)==0
        s=s+moebius(n/i)*(fibonacci(i+1)+fibonacci(i-1)); 
    end
end
a=s/n
 

Вот он же в виде функции:

Используется синтаксис Matlab M

function [a] = calcType1Long (n)
    n=sym(n);
    s=0;
    for i=1:n
        if mod(n,i)==0
            s=s+moebius(n/i)*(fibonacci(i+1)+fibonacci(i-1)); 
        end
    end
    a=s/n;
end
 

Дает такой ответ:

Используется синтаксис Matlab M

>> calcType1Long(89)
 
ans =
 
3980154972736925972/89 ~= 3980154972736918051/89
 

Целое значение никак не получается вывести(((

-- 27.04.2020, 14:41 --

ShMaxG в сообщении #1458199 писал(а):

maximkarimov
Если без MuPad, можете попробовать использовать функцию vpa.

Пробовал - не получается.

[GAA]

Поставьте точку останова на строчке a=s/n; запустите вычисление calcType1Long(89). Как перейдёте на a=s/n; посмотрите значения s и n и что за объекты.
Если значения не те, или объекты не символьные, то надо смотреть выше по тексту функции, аналогичным способом: ставить точку останова и проверять.

-- Mon 27.04.2020 12:52:19 --

Так Вы быстро найдёте место с ошибкой без посторонней помощи.

[maximkarimov]

GAA в сообщении #1458212 писал(а):

Поставьте точку останова на строчке a=s/n; запустите вычисление calcType1Long(89). Как перейдёте на a=s/n; посмотрите значения s и n и что за объекты.
Если значения не те, или объекты не символьные, то надо смотреть выше по тексту функции, аналогичным способом: ставить точку останова и проверять.

-- Mon 27.04.2020 12:52:19 --

Так Вы быстро найдёте место с ошибкой без посторонней помощи.

Все переменные отображаются в Worckspace как символьные.

[GAA]

GAA в сообщении #1458212 писал(а):

посмотрите значения s и n

[maximkarimov]

GAA в сообщении #1458216 писал(а):

GAA в сообщении #1458212 писал(а):

посмотрите значения s и n

Используется синтаксис Matlab M

>> s
 
s =
 
3980154972736925972 ~= 3980154972736918051
 
>> n
 
n =
 
89
 
 

s почему-то считается приблизительно.

-- 27.04.2020, 15:38 --

GAA в сообщении #1458216 писал(а):

GAA в сообщении #1458212 писал(а):

посмотрите значения s и n

Покажите как у Вас реализовано вычисление функции Мебиуса, плиз. Может у меня с ней проблема? Потому как при замене ее на функцию Эйлера точно такой же код работает прекрасно - считает члены последовательности А000358.

[GAA]

Я не реализовывал, а воспользовался имеющейся. R2013b:

Используется синтаксис Matlab M

 for i=1:60, A(i) = feval(symengine,'numlib::moebius', i); end;
>> A
A =

[ 1, -1, -1, 0, -1, 1, -1, 0, 0, 1, -1, 0, -1, 1, 1, 0, -1, 0, -1, 0, 1, 1, -1, 0, 0, 1, 0, 0, -1, -1, -1, 0, 1, 1, 1, 0, -1, 1, 1, 0, -1, -1, -1, 0, 0, 1, -1, 0, 0, 0, 1, 0, -1, 0, 1, 0, 1, 1, -1, 0]
 

-- Mon 27.04.2020 14:05:38 --

Используется синтаксис Matlab M

>> feval(symengine,'numlib::moebius', 89)
ans =
-1

-- Mon 27.04.2020 14:12:44 --

Вашу m-функцию подправить можно для R2013b так:

Используется синтаксис Matlab M

function [a] = calcType1Long (n)
    s=0;
    i= sym('i');
    for i=1:n
     if mod(n,i)==0
        s=s+feval(symengine,'numlib::moebius', n/i)*(feval(symengine,'numlib::fibonacci', i+1)+feval(symengine,'numlib::fibonacci', i-1)); 
     end
    end
    a=s/n;
end

Тогда

Используется синтаксис Matlab M

>> calcType1Long(89)
ans =
44720842390302450

[maximkarimov]

УРА!!!
Проблема была именно в "кривой" реализации функции Мебиуса, которую я скачал с матлабовского файлообменника. Значения она выдавала правильные, но там использовалась функция floor и потому на символьных вычислениях выскакивала проблема.
Переписал ее сам (она ведь простенькая!) и все получилось:

Используется синтаксис Matlab M

>> calcType1Long(89)
 
ans =
 
44720842390302450
 

Спасибо всем!

[maximkarimov]

GAA в сообщении #1458237 писал(а):

Вашу m-функцию подправить можно для R2013b так:

Используется синтаксис Matlab M

function [a] = calcType1Long (n)
    s=0;
    i= sym('i');
    for i=1:n
     if mod(n,i)==0
        s=s+feval(symengine,'numlib::moebius', n/i)*(feval(symengine,'numlib::fibonacci', i+1)+feval(symengine,'numlib::fibonacci', i-1)); 
     end
    end
    a=s/n;
end

Тогда

Используется синтаксис Matlab M

>> calcType1Long(89)
ans =
44720842390302450

С Вашей помощью нашел ошибку - "кривая" реализация функции Мебиуса, поэтому исправил именно ее, а код оставил прежним:

Используется синтаксис Matlab M

function [a] = calcType1Long (n)
    n=sym(n);
    s=0;
    for i=1:n
        if mod(n,i)==0
            s=s+moebius(n/i)*(fibonacci(i+1)+fibonacci(i-1)); 
        end
    end
    a=s/n;
end
 

Результат тот же:

Используется синтаксис Matlab M

>> calcType1Long(89)
 
ans =
 
44720842390302450
 

Еще раз спасибо!

[GAA]

В R2013b приведенную в начальном сообщении функцию moebius красиво и просто переделать под символьные вычисления у меня не получилось. (В частности, при стандартном определении (function mu = moebius(n)... ) тип входного аргумента будет double array.)

На форум заглядывают участники и с новыми версиями Matlab. Поэтому, если у Вас получилось её довести до рабочего состояния, то выкладывайте текст (функции moebius). Возможно, работающий вариант будет интересен участникам.

[maximkarimov]

GAA в сообщении #1458277 писал(а):

На форум заглядывают участники и с новыми версиями Matlab. Поэтому, если у Вас получилось её довести до рабочего состояния, то выкладывайте текст (функции moebius). Возможно, работающий вариант будет интересен участникам.

Вот:

Используется синтаксис Matlab M

function [a] = moebius(n)
if n==1
    a=1;
else   
    v=factor(n);
    v1=unique(v);
    a=-1;
    x=0;
    for i=1:size(v1,2)
      if mod(n,v1(i)^2)~=0 
          x=x+1;
      end
    end
    if x~=size(v1,2)
      a=0;
    elseif mod(size(v1,2),2)==0
      a=1;
    end
end
end
 

[GAA]

Обозначим эту [вторую] модификацию через moebius2.

Используется синтаксис Matlab M

>> for i=1:200, A(i) = moebius(i); end
>> for i=1:200, B(i) = moebius2(i); end
>> any(A-B)
ans =
     0 

Один и тот же результат возвращают.

Такой вариант с текстом moebius из начального сообщения ветки

Используется синтаксис Matlab M

function [a] = calcType1Long (n)
    n = sym(n);
    s=0;
   for i=1:n
     if mod(n,i)==0
        s=s+moebius(n/i)*(feval(symengine,'numlib::fibonacci', i+1)+feval(symengine,'numlib::fibonacci', i-1)); 
     end
    end
    a=s/n;
end

в R2013b даёт

Используется синтаксис Matlab M

>> calcType1Long(89)
ans =
44720842390302450