Ошибка Mathematica CForm / FortranForm?

Question

Я пытаюсь вывести очень сложную матрицу (~ 1,3 МБ в виде обычного текста) из Mathematica для использования в программе на Фортране.Когда я делаю это (через Splice), результирующая матрица отключается на ~ 2%, когда переменные получают числовые значения.Это проблема, так как должно быть собственное значение, которое точно равно нулю, и состав собственных векторов должен быть точно правильным.

Я сделал все обычные должные проверки в отношении точности, правильных переменных, правильныхкод диагонализации и т. д. И все сводилось либо к тому, что сам Фортран не смог справиться с такой большой матрицей, либо к Mathematica испортили вывод FortranForm.

Поэтому я заставил Mathematica дать мне CForm матрицы и попробовал,Кроме того, было ~ 2% от того, что должно было быть, что более поразительно, это было то же самое (с точностью до машины), что и матрица FortranForm!

Кто-нибудь сталкивался с такой проблемой?У вас есть идеи, что может вызвать это?Я боюсь пройти через 25000 строк кода Fortran, отформатированного в Mathematica, чтобы выяснить это.

РЕДАКТИРОВАТЬ: Матрица, о которой идет речь, сложная, а не большая.Это всего лишь 6x6, но каждый элемент индивидуально алгебраически очень беспорядочный, включая тригонометрические функции, логарифмы и различные корни и степени.

Открытый текст элемента (1,1) нашей матрицы,код C и код Fortran .Значения параметров в здравом уме: 0 <лямбда, каппа, Y *** <1;все остальные между 100 и 1000. </p>

Answer 1

Глядя на выражения для первого элемента, это почти наверняка проблема из-за арифметики с плавающей запятой конечной точности. Мне кажется, есть два способа попытаться подойти к этой проблеме. Первый заключается в использовании большей точности в коде Си. Возможно, вы могли бы использовать препроцессор макросов (или просто найти и заменить), чтобы изменить все объявления в коде C с double на long double. В зависимости от того, какой компилятор вы используете, вы получите от 64 бит до 80 или 128 (если вы не используете компилятор Microsoft , который рассматривает длинные двойные числа как стандартные двойные). Если вы опубликуете, какой компилятор C вы используете, я могу помочь посмотреть, какие изменения вам нужно будет сделать, чтобы пойти по этому пути. В самом конце вы можете привести к удвоению, если вам нужно передать его Mathematica или другой программе на Си, которая просто использует double.

Кроме того, компиляторы C могут иметь разные режимы с плавающей запятой. Возможно, стоит проверить, установлен ли режим для точного, строгого или быстрого режима. Если это быстрый или быстрый математике, то это может способствовать вашим проблемам с округлением.

Другой путь заключается в использовании Mathematica для перестановки членов каждого элемента с целью оптимизации точности с плавающей запятой. В то время как в целом (a * b) * c = (c * b) * a в стране арифметики с плавающей запятой это не так. Лучшая эвристика, о которой я слышал, это:

1. Добавление: добавьте термины от наименьшего к наибольшему, чтобы минимизировать ошибки округления
2. Вычитание: избегайте или задерживайте вычитание значений, которые по значению очень близки друг к другу, чтобы сохранить как можно больше значащих цифр
3. Умножение / деление: избегайте умножения / задержки умножения на очень малые числа или деления на очень большие числа по причинам значимых цифр.

На практике это означает, что факторинг этих выражений - ваш лучший выбор, особенно проблемные. если a1 ~ = a2, то выражение (a1-a2) (b-d) должно иметь больше точности / значащих цифр, чем (a1-a2) * b- (a1-a2) * d, потому что у вас на одну операцию умножения меньше.

Answer 2

Может ли быть следующая простая вещь: Mathematica довольно умна в числовой оценке, поэтому, например, эти два выражения

N[(10^18+100)-10^18]
N[10^18+100]-N[10^18]

приводят к совершенно разным выводам

Out[3]= 100.
Out[4]= 128. 

с Cили в Фортране вы не получаете няни, и выражение будет оцениваться как во второй строке.
Может быть, вы могли бы посмотреть на одну запись матрицы и попытаться оценить ее в Mathematica, как в C / Fortran, и посмотреть, если этосоответствует тому, что вы видите?

РЕДАКТИРОВАТЬ:
Кажется, что простой тест для проблем с промежуточной точностью должен использовать Compile, который отбрасывает все проверки точности, так что

In[10]:= Compile[{a,b},a+100-b][10^18,10^18]
Out[10]= 128.

Не могли бы вы проверить,скомпилированная версия вашей матрицы согласна с C / Fortran?