Документация

Использовать символьные вычисления, когда это возможно

Рекомендуется выполнять вычисления символически, поскольку точные символьные вычисления не подвержены ошибкам округления. Например, стандартные математические константы имеют собственные символьные представления в символьных математических Toolbox™:

pi
sym(pi)

ans =
    3.1416

ans =
pi

Избегайте ненужного использования числовых аппроксимаций. Число с плавающей запятой аппроксимирует константу; это не сама константа. С помощью этой аппроксимации можно получить неверные результаты. Например, heaviside специальная функция возвращает различные результаты для синуса sym(pi) и синус числового приближения pi:

heaviside(sin(sym(pi)))
heaviside(sin(pi))

ans =
1/2

ans =
     1

Выполнение расчетов с повышенной точностью

Гипотеза Римана утверждает, что все нетривиальные нули дзета-функции Римана (z) имеют одинаковую действительную часть ℜ (z) = 1/2. Чтобы найти возможные нули функции Зета, постройте график её абсолютного значения | (1/2 + iy) |. На следующем графике показаны первые три нетривиальных корня функции Зета | (1/2 + iy) |.

syms y
fplot(abs(zeta(1/2 + i*y)), [0 30])

Использование числового решателя vpasolve для аппроксимации первых трех нулей этой функции Зета:

vpasolve(zeta(1/2 + i*y), y, 15)
vpasolve(zeta(1/2 + i*y), y, 20)
vpasolve(zeta(1/2 + i*y), y, 25)

ans =
14.134725141734693790457251983562
 
ans =
21.022039638771554992628479593897
 
ans =
25.010857580145688763213790992563

Теперь, рассмотрите ту же функцию, но немного увеличьте реальную часть, $$\zeta \frac{\mathrm{}}{\mathrm{}}(10000000012000000000+iy$$). Согласно гипотезе Римана, эта функция не имеет нуля для какого-либо действительного значения y. Если вы используетеvpasolve с 10 значащими десятичными цифрами решатель находит следующий (несуществующий) ноль функции Зета:

old = digits;
digits(10)
vpasolve(zeta(1000000001/2000000000 + i*y), y, 15)

ans =
14.13472514

Увеличение числа цифр показывает, что результат неверен. У функции Дзэты $$\zeta \frac{\mathrm{}}{\mathrm{}}$$(10000000012000000000+iy) нет ноля ни для какой реальной стоимости 14 <y <15:

digits(15)
vpasolve(zeta(1000000001/2000000000 + i*y), y, 15)
digits(old)

ans =
14.1347251417347 + 0.000000000499989207306345i

Для дальнейших вычислений восстановите количество цифр по умолчанию:

digits(old)

Сравнение символьных и числовых результатов

Функции Бесселя с полуцелыми индексами возвращают точные символьные выражения. Аппроксимация этих выражений числами с плавающей запятой может привести к очень нестабильным результатам. Например, точное символическое выражение для следующей функции Бесселя:

B = besselj(53/2, sym(pi))

B =
(351*2^(1/2)*(119409675/pi^4 - 20300/pi^2 - 315241542000/pi^6...
 + 445475704038750/pi^8 - 366812794263762000/pi^10 +...
 182947881139051297500/pi^12 - 55720697512636766610000/pi^14...
 + 10174148683695239020903125/pi^16 - 1060253389142977540073062500/pi^18...
 + 57306695683177936040949028125/pi^20 - 1331871030107060331702688875000/pi^22...
 + 8490677816932509614604641578125/pi^24 + 1))/pi^2

Использовать vpa для аппроксимации этого выражения с 10-значной точностью:

vpa(B, 10)

ans =
-2854.225191

Теперь вызовите функцию Бесселя с параметром с плавающей запятой. Значительная разница между этими двумя приближениями указывает на то, что один или оба результата неверны:

besselj(53/2, pi)

ans =
   6.9001e-23

Увеличение числовой рабочей точности для получения более точной аппроксимации для B:

vpa(B, 50)

ans =
0.000000000000000000000069001456069172842068862232841396473796597233761161

Постройте график функции или выражения

Печать результатов поможет распознать неправильные аппроксимации. Например, числовое приближение следующей функции Бесселя возвращает:

B = besselj(53/2, sym(pi));
vpa(B, 10)

ans =
-2854.225191

Постройте график этой функции Бесселя для значений x вокруг 53/2. График функции показывает, что аппроксимация неверна:

syms x
fplot(besselj(x, sym(pi)), [26 27])