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

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

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

pi
sym(pi)

ans =
    3.1416

ans =
pi

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

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

ans =
1/2

ans =
     1

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

Гипотеза Римана утверждает, что все нетривиальные нули дзета-функции Римана (z) имеют одинаковую действительную часть ℜ (z) = 1/2. Чтобы найти возможные нули функции Zeta, постройте график ее абсолютного значения |ζ (1/2 + iy) |. Следующий график показывает первые три нетривиальных корня функции Zeta |ζ (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 значащих десятичных цифр решатель находит следующий (несуществующий) нуль функции Zeta:

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

ans =
14.13472514

Увеличение количества цифр показывает, что результат неправильен. Функция Zeta $$\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])