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

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

Выполнение вычислений символически рекомендуется, потому что точные символьные вычисления не подвержены ошибкам округления. Например, стандартные математические константы имеют свои собственные символьные представления в Symbolic Math 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 \left(\frac{1000000001}{2000000000}+iy\right)$$. Согласно Римановой гипотезе, эта функция не имеет нуля ни для какого действительного значения y. Если вы используете vpasolve с 10 значительными десятичными цифрами, решатель находит следующий (несуществующий) нуль Дзета-функции:

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

ans =
14.13472514

Увеличение числа цифр показывает, что результат является неправильным. Дзета-функция $$\zeta \left(\frac{1000000001}{2000000000}+iy\right)$$ не имеет нуля ни для какого действительного значения 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])