Умножьте две переменные в Matlab с помощью vpa — высокая точность

Я хочу быть уверенным, что две переменные, a и b, умножаются с высокой точностью, т. е. произведение c из a и b с произвольной точностью (в моем случае 50 правильных десятичных цифр):

a = vpa(10/3,50);  
b = vpa(7/13,50);  
c = eval(vpa(vpa(a,50)*vpa(b,50),50));      % I basically want to do just c = a*b;

что дает мне

a = 3.3333333333333333333333333333333333333333333333333
b = 0.53846153846153846153846153846153846153846153846154
c = 23.333333333333333333333333333333

Тестирование

d = eval(vpa(c*13,50))

дает

d = 23.333333333333333333333333333333333333335292490585

что показывает, что умножение для получения c не было выполнено с 50 значащими цифрами.

Что здесь не так, но, что более важно, как мне получить правильный результат для a*b и для других операций, таких как exp?


matlab numeric symbolic-math arbitrary-precision
person Futurist    schedule 07.06.2015    source источник

Ответы (1)


arrow_upward
1
arrow_downward

Во-первых, следует использовать vpa('7/13',50) или vpa(7,50)/13, чтобы избежать возможности потери точности из-за того, что 7/13 вычисляется с плавающей запятой двойной точности (я полагаю, что vpa, как и sym, пытается угадать обычные константы и рациональные дроби, но на это не стоит полагаться).

Проблема в том, что хотя a и b хранятся как 50-значные значения переменной точности, ваше умножение по-прежнему выполняется в соответствии со значением по умолчанию digits (32). Второй аргумент vpa указывает только точность переменной, а не какие-либо последующие операции над ней или с ней (документация в этом отношении не особенно полезна).

Один из способов выполнить то, что вы хотите, будет:

old_digits = digits(50);
a = vpa('10/3')
b = vpa('7/13')
c = a*b
d = 13*c
digits(old_digits);

Другим было бы использование точных символьных выражений для всей математики (потенциально более дорого), а затем преобразование результата в 50-значную переменную точность в конце:

a = sym('10/3')
b = sym('7/13')
c = a*b
d = vpa(13*c,50)

Оба метода возвращают 23.333333333333333333333333333333333333333333333333 вместо d.

person horchler    schedule 08.06.2015
comment
Хорошо спасибо; просто сомнительный вопрос: как я могу быть уверен, что в более сложном вычислении vpa действительно вычисляет 50 десятичных цифр, а не просто показывает их, внутренне представляя их в виде рациональных дробей или других символьных выражений? Есть ли простое численное доказательство? - person Futurist; 10.06.2015
comment
@Futurist: Если ваш код правильный, он должен. Ваш вопрос может быть задан о любых вычислениях. В конце концов, единственный способ доказать это самому себе — заранее знать ответ. С практической точки зрения вы можете попытаться выполнить те же вычисления символически и/или использовать большее значение для digits, чтобы проверить, как изменится результат. Или вы можете выполнить вычисления в другом программном обеспечении (я часто проверяю свои результаты в Mathematica или Wolfram Alpha) — в некоторых случаях разные алгоритмы могут обнаруживать (или не обнаруживать) разные решения или выполнять арифметические действия по-разному. - person horchler; 11.06.2015
comment
Вы также можете попробовать ограничить свое решение, используя интервальную арифметику — см. INTLAB Toolbox для Matlab. - person horchler; 11.06.2015