вычисление корней определителя в Матлаб

[volchenok]

Здравствуйте, уважаемые участники форума. Мне нужна помощь по быстрому вычислению корней определителя. Задача такая: дана матрица большого размера, которая является функцией комплексного аргумента. Нужно найти такие значения этого аргумента, чтоб ее определитель обращался в ноль. Задачу я решаю методом Мюллера, так как аналитически тут нельзя ничего найти, матрица большая, все нелинейно зависит от искомых корней. Ну и в процессе решения (как в принципе любого численного поиска корней) мне нужно вычислять значение определителя в некоторых точках и так получается что это значение огромно. Матлаб воспринимает это как бесконечность. Вот в этом вопрос и состоит. Естественным образом я попытался делить значение определителя на его же значение в начальном приближении, но то что в знаменателе то же считается матлабом бесконечностью. Я тогда внес знаменатель под знак определителя (т.е. перемножал две матрицы под знаком определителя) и в итоге нашел все корни. Но цена этого ооооочень большое время счета и именно в месте перемножения матриц. Найденные корни я пытался брать в качестве начальных приближений, но и они без нормировки на значение определителя в начальной точке, не обращали в ноль определитель и он все еще считался матлабом бесконечностью. В общем, найти корни у меня получилось только путем домножения матрицы на обратнуюю в начальном приближении и последующем взятием определителя, но это приводит огромной трате ресурсов. Подскажите пожалуйста другие способы выхода из этой ситуации. Я уверен, проблема эта не нова и должны существовать уже проверенные методы борьбы с ней. Спасибо.

[B@R5uk]

В чём проблема поделить все элементы на константу?

[volchenok]

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

[B@R5uk]

volchenok в сообщении #1274100 писал(а):

она либо маленькая, что не влияет на огромность определителя, либо большая, настолько, что матлаб уже считает начальное приближение корнем

Сдаётся мне, что не правильно ваша программа написана. Поскольку у вас алгоритм Мюллера, то на одну и ту же константу надо делить три определителя. Потом, вы пишите:

volchenok в сообщении #1274014 писал(а):

Задача такая: дана матрица большого размера, которая является функцией комплексного аргумента.

У вас в процессе вычисления есть комплексные числа или же сам аргумент, по которому ищется ноль, является комплексным числом? В последнем случае задача эффективно становится двумерной, ответ — континуальным множеством (фигура на плоскости), а алгоритм Мюллера — не применим.

И да, было бы здорово увидеть вашу программу. А то мы о каких-то воздушных материях разговариваем. Если же программа — секрет, то попытайтесь оценить значение определителя через максимальный по модулю элемент матрицы. Его надо возвести в степень размера матрицы. Это число должно быть меньше realmax, который, если верить Гуглу, равен 1.797693134862316e+308. Вычисляете из этого числа корень степени размера матрицы и сравниваете с максимальным элементом. Если он того же порядка, то вероятно выйдет переполнение при вычислении определителя. Матрицу можно поделить на этот максимальный элемент, тогда результат точно будет не больше единицы.

-- 12.12.2017, 11:13 --

volchenok в сообщении #1274014 писал(а):

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

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

[volchenok]

Цитата:

Сдаётся мне, что не правильно ваша программа написана. Поскольку у вас алгоритм Мюллера, то на одну и ту же константу надо делить три определителя

Безусловно, я считаю три определителя, вернее один, но в трех разных точках, метод Мюллера, он же метода парабол, а парабола определяется значениями в трех точках. Этот метод чуть ли не единственный, который позволяет находить комплексные корни, на Фортране библиотека именно его реализовала. В других своих задачах я проверял этот метод и в его правильности сомнений нет. Конечно в данном случае комплексные и корни и значение самого определителя, но это ничему не противоречит.

Цитата:

И да, было бы здорово увидеть вашу программу. А то мы о каких-то воздушных материях разговариваем. Если же программа — секрет, то попытайтесь оценить значение определителя через максимальный по модулю элемент матрицы. Его надо возвести в степень размера матрицы. Это число должно быть меньше realmax, который, если верить Гуглу, равен 1.797693134862316e+308. Вычисляете из этого числа корень степени размера матрицы и сравниваете с максимальным элементом. Если он того же порядка, то вероятно выйдет переполнение при вычислении определителя. Матрицу можно поделить на этот максимальный элемент, тогда результат точно будет не больше единицы.

Дело не в секретности, вернее не сколько в ней, сколько в большой сложности кода. Если оставить толкьо саму матрицу и код вызова ее определителя, то останется набор непонятных подфункций которые прийдется обьяснять, сути вопроса это все равно не изменит. Задача, мне кажется, без кода и так понятна: есть матрица, с комплексным аргументов, нужно найти корни ее определителя. Решается методом парабол (Мюллера). Хотя в любом методе нужно считать значение функци в каких-то точках. Относительно максимального числа Вы правы, Гугл не врет, ну вот в моей матрице максимальные элемент порядка 1000, а если проделать процедуру, описанную вами, для сравнения, то получается число порядка единицы. Хотя я и так выше написал, что проблема в большой величине определителя. При делении на максимальный элемент выдало ожидаемый результат, что корнем является начальное приближение. Я уже пробовал делить на меньшие значения матрицу и получал такой же ответ, при 250 где-то происходит переход на очень большие значения определителя.

Цитата:

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

К сожалению все не так просто. Я уже писал, что вычислить напрямую ни значения определителя, ни значения обратной матрицы я не могу - я выхожу за рамки матлаба. Если б я смог вычислить значение определителя обратной матрицы, тогда вопроса бы не было потому что вопрос заключается в избавлении от умножения матриц, которым я пользуюсь, чтоб посчитать определитель. Я заменяю отношение двух определителей на определитель произведение исходной матрицы и обратной в начальном приближении.

Но все равно спасибо Вам за советы. Я буду рад выслушать еще, так как у меня тоже пока идей не появилось, а вопрос еще актуален.

[B@R5uk]

volchenok в сообщении #1274397 писал(а):

При делении на максимальный элемент выдало ожидаемый результат, что корнем является начальное приближение.

Так и проверку окончания работы тоже поправить надо!

volchenok в сообщении #1274397 писал(а):

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

Какой размер матрицы? Какова максимальная величина у модулей элементов матрицы? Это мне для справки. А по проблеме можно сделать следующее. Из каждого столбца вынесите положительный множитель так, чтобы значения элементов в каждом столбце не превышали по модулю единицу. Если верить вам, то произведение этих множителей превысит 308-ю степень десятки. Поэтому, их надо не перемножать, а прологарифмировать и сложить. Получится расширение представления числа с плавающей точкой, вернее расширение экспоненциальной части этого представления. При желании даже можно заморочиться и выносить в качестве множителей только степени двойки.

[volchenok]

Цитата:

Так и проверку окончания работы тоже поправить надо!

Как таковой, на данном этапе, проверки нет. Я пока вручную останавливаю счет, когда значения модуля функции становятся достаточно маленькими и когда модуль разности корней на соседних итерациях тоже малы. Но это пока не суть, так как хочется чтоб расчет хотя бы начался. То что матлаб считает корнем начальное приближение при достаточно больших константах деления видно, когда определитель уже в начальном приближении сразу ноль (а не какая-то малая величина).

Цитата:

Какой размер матрицы? Какова максимальная величина у модулей элементов матрицы? Это мне для справки. А по проблеме можно сделать следующее. Из каждого столбца вынесите положительный множитель так, чтобы значения элементов в каждом столбце не превышали по модулю единицу. Если верить вам, то произведение этих множителей превысит 308-ю степень десятки. Поэтому, их надо не перемножать, а прологарифмировать и сложить. Получится расширение представления числа с плавающей точкой, вернее расширение экспоненциальной части этого представления. При желании даже можно заморочиться и выносить в качестве множителей только степени двойки.

Размер матрицы 2000 на 2000 , но в будущем планируется еще большим. Максимальная величина (по модулю), как я уже писал, равна тысячи. Минимальная равна нулю. Матрица, к сожалению, не имеет особой структуры, не содержит много нулевых элементов.
Относительно совета, то правильно ли я вас понял, что нужно вынести с каждого столбца максимальный (по модулю) элемент, потом их прологарифмировать (натуральным логарифмом), сложить и в расчетах делить матрицу на полученную сумму?
Еще раз спасибо, что уделяете время.

[B@R5uk]

volchenok в сообщении #1274507 писал(а):

и в расчетах делить матрицу на полученную сумму

Нет. Сумма логарифмов — это показатель степени, на которую надо умножить значение определителя получившейся матрицы (с элементами не большими единицы), чтобы получить то число, которое выходит за пределы возможностей матлаба. Поскольку это число выходит за пределы, то, разумеется, умножать на эту степень не надо. Но во всех последующих вычислениях со значением определителя этот показатель, очевидно, надо учитывать. Фактически, вы получите число, представленное в виде $A\cdot\exp(B)$, или в виде $A\cdot10^B$, или в виде $A\cdot2^B$. Операции над числами, представленными в таком виде придётся написать вручную. Контроль за точностью вычислений тоже желательно осуществить. Потому что $1 + eps = 1$, не смотря на то, что $eps \ne 0$.

В связи с существованием вопроса точности вычислений, я крайне рекомендую вам в качестве множителя из столбцов выносить степени двойки, тогда двоичных логарифм этой степени всегда будет целым числом (на практике, можно составить таблицу). Умножение/деление на степени двойки никогда не приводят к потере точности в типах с плавающей точкой, даже когда десятичное представление имеет вид типа $1,0633824e+37$. В то время как последовательность операций логарифмирования/экспонирования над числами с плавающей точкой точность снижают весьма заметно. Это может быть весьма критично в алгоритмах поиска какого-либо значения последовательными приближениями, когда надо вычитать близкие по величине числа.

[volchenok]

Спасибо, буду пробовать реализовать операции с такими числами, надеюсь расчеты быстрее с ними будут, чем перемножение матриц. Если будут у Вас какие-то другие мысли по этому поводу, то буду рад выслушать. О своих мыслях и результатах я тоже напишу. Спасибо Вам.