CC BY
21ВЫСОКОМОЛЕКУЛЯРНЫЕ СОЕДИНЕНИЯ, Серия Б, 2005, том 47, № 1, с. 121-143
- ДИСКУССИИ
УДК 541.64:539(199+3)
К ВОПРОСУ ОБ УСТОЙЧИВОМ ОПРЕДЕЛЕНИИ РЕЛАКСАЦИОННОГО СПЕКТРА ИЗ ДАННЫХ ПО МЕХАНИЧЕСКОЙ РЕЛАКСАЦИИ ПОЛИМЕРОВ1
© 2005 г. В. А. Дубовицкий, В. И. Иржак
Институт проблем химической физики Российской академии наук 142432 Черноголовка Московской обл., пр. Ак. Семенова, I Поступила в редакцию 13.05.2003 г.
Принята в печать 27.06.2004 г.
В работе предложен новый формализм для описания релаксационных кривых, основанный на понятии интеграла Стилтьеса относительно функции распределения экспоненциальных мод по временам релаксации. Этот подход позволяет единообразно описать все возможные разновидности релаксационного спектра, включая дискретный и непрерывный типы, начальный мгновенный спад кинетической кривой и ее стационарное плато. Важным обстоятельством является наличие взаимно однозначного непрерывного соответствия между монотонной функцией распределения по временам релаксации и временной или частотной зависимостью релаксационного модуля упругости. Благодаря этому обратная задача восстановления распределения по временам релаксации является математически корректной в классическом смысле и не требует привлечения техники регуляризации. Описан эффективный алгоритм численного решения. Приведены примеры обработки модельных кинетических данных для модифицированной цепи Рауза, а также экспериментальных данных по частотной зависимости модуля упругости линейных полимеров. Таким образом, предложен способ нахождения релаксационного спектра из данных по механической релаксации полимеров, точность которого лимитирована лишь адекватностью математической модели, а также величиной систематической ошибки и полнотой этих данных.
ВВЕДЕНИЕ
В настоящее время признанным способом описания вязкоупругих свойств полимеров в линейной области их механического поведения является построение распределения максвелловых экспоненциальных мод ехр(-//т) по временам релаксации т, т.е. релаксационного спектра (PC). С одной стороны, предполагается, что PC коррелирует со строением макромолекул и их окружением [1]. С другой стороны, через него при помощи простых интегральных преобразований выражаются важнейшие непосредственно измеряемые характеристики полимерных материалов: кинетические функции напряжения и деформации, соответствующие частотные модули, мгновенный модуль упругости, ньютоновская вязкость и т.д. [2]. Зная PC, можно в широких пределах с доста-
1 Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований и регионального фонда р2001 Подмосковье (коды проектов 01-03-97001 и 02-01-00538).
E-mail: dubv@icp.ac.ru (Дубовицкий Владимир Абрамович).
точной точностью описать поведение полимера при произвольном законе деформирования или нагружения. Именно поэтому так важно иметь возможность определять РС конкретного материала, исходя из данных ограниченного эксперимента. И это тем более важно, что современная экспериментальная техника достигла высокого совершенства, позволяя определять сотни значений с относительной точностью до 10"4 в диапазоне более 14 декад по времени или частоте.
Исторически первым является представление РС в виде конечного набора бесконечно узких линий с определенными временами релаксации и ве-сами-интенсивностями. В непрерывном формализме РС описывается в терминах функции плотности р(т) распределения интенсивности мод по временам релаксации.
В литературе [3-6] неоднократно обсуждалось, что принятые методы нахождения РС грешат изначально предопределенным произволом, выявляющимся при попытках приложения моде-
лей к интерпретации конкретных экспериментальных данных. Так, в дискретном варианте, реализуя различные процедуры аппроксимации, легко получить неограниченно много разных комплектов времен релаксации и весов дискретных мод, которые описывают с точностью до шума экспериментальные кривые, однако эти наборы не имеют видимой корреляции между собой и структурой полимера. При непрерывном способе описания возникает некорректно поставленная задача восстановления функции плотности р(т), как решения линейного интегрального уравнения первого рода. Такие уравнения традиционно решают методами регуляризации, при этом произвол восстановления плотности маскируется способом организации процесса регуляризации и дискретизации интегрального уравнения (выбор окна, сетки и т.д.). Вместе с тем замечено [6], что, несмотря на видимый произвол описания РС, интегральные его характеристики зачастую ведут себя согласованно. Сложилась странная ситуация, которая отчетливо сформулирована в работе [4]: "...выбор способа описания РС и метода его вычисления определяет персональное предпочтение, а не объективная реальность". В недавнем обзоре [6] пессимистическая картина неуловимости РС сформулирована еще яснее, сделан вывод: "РС - не более чем способ аппроксимации в пределах точности измерений, удобный для различных операций в рамках теории вязкоупругости сплошных сред" и приведены многие ссылки подобного рода.
Уместно подчеркнуть, что в цитируемой дискуссии смешаны две совершенно независимые проблемы. Одна относится к физической адекватности математических моделей, описывающих взаимосвязи вязкоупругих характеристик и основанных на понятии РС. Другая состоит в математически корректной записи этих представлений. Понятно, что нелогично критиковать и отвергать удобное фундаментальное понятие о РС, исходя только из наблюдающихся в вычислительной практике ненадежных результатов его восстановления в рамках неустойчивых форм описания.
Итак, явно назрел вопрос о том, существует ли вообще устойчивый способ описания РС, и не связана ли наблюдающаяся неустойчивость с неполнотой традиционных дискретной и непрерывных математических моделей?
Разрешение отмеченных противоречий может быть получено на основе развиваемого одним из авторов настоящей статьи общего метода оптимальных интегральных представлений (ОИП) корректной переформулировки и решения обратных задач спектроскопии. Ранее этот подход применялся для иных приложений [7, 8]. Логически метод тесно связан с математическим аппаратом проблемы моментов [9] и теории Шоке [10]. Применительно к данной работе изложение очень упрощается благодаря тесной связи обсуждаемых интегральных преобразований с преобразованием Фурье. Можно сказать, что задачи восстановления РС являются удобным полигоном для демонстрации общей теории ОИП. Этому и посвящена настоящая работа.
Особенностью статьи является параллельное изложение восстановления РС по кинетическим и частотным функциям. Статья построена таким образом, что в ее первой части дана только идея устойчивого описания и восстановления РС с использованием интеграла Стилтьеса, и поэтому с точки зрения математики эта часть грешит недостаточной точностью изложения. Предложенный способ проиллюстрирован рядом примеров.
Поскольку рассматриваемые вопросы являются относительно сложными и тонкими, во второй части статьи приведены математически более корректные и точные формулировки (Приложение). Значительно более деликатные вопросы аналитических асимптотик ошибки восстановления РС в статье не затронуты.
ПРИМЕНЕНИЕ ИНТЕГРАЛА СТИЛТЬЕСА ДЛЯ ОПИСАНИЯ И ВОССТАНОВЛЕНИЯ РЕЛАКСАЦИОННОГО СПЕКТРА
Экспериментальные данные по механической релаксации обычно получают в виде зависимости кинетического модуля g(t) от времени или соответствующего комплексного модуля С(со) от частоты. Используя представления [1, 2] о РС, эти функции описывают в дискретном варианте как
N
8(0 = Р~+ X ехР(~'/х*)Р*
Т-^ «СОТл
С(со) = у --
г ^ 1 + 10)Т*
к = 1 *
Здесь {хк, рк}, к = 1,..., N есть набор N времен релаксации и весов соответствующих мод Максвелла, а р„- вес инертной составляющей. В непрерывном варианте
g(t) = р„ + |ехр(-г/х)р(х)</х
(2)
G(t0) = p~ + iTT^p(t)£/x'
о
где р(т) - некоторая интегрируемая функция, имеющая смысл плотности распределения экспоненциальных мод по времени релаксации. При этом по своему физическому смыслу числа хк, рк, р(х), неотрицательны. Кинетическая и частотная функции связаны друг с другом через преобразование Фурье формулой G((о) = /coF[g(i)]((0), т.е.
G( о)) = pM + /Mjexp(-i7M)(g(i)-pJift (3)
о
Представления (1) и (2) допускают естественное объединение на языке понятия функции распределения и конструкции интеграла Стилтьеса [11,12].
Введем функцию распределения Р(т), обладающую следующими свойствами: она определена на расширенной неотрицательной полуоси /?+ = = [0, непрерывна слева при 0 < х < и удовлетворяет равенству Р(0) = 0. Мы можем выразить теперь функции g(t), G(cо) при помощи интеграла Стилтьеса относительно Р(х):
g{t) = Jexp (-tlx)dP(x),
(4)
©о
G(ö» = fr^b^(x)
J 1 + iCOX
о
Дискретному случаю (1) на языке интеграла Стилтьеса соответствует ступенчатая функция распределения, принимающая для конечных х > 0 значения
N
Р(Х).= £ptA(T-Tt),
к = 1
где h - единичная функция Хевисайда, а для х = <=°
Р(х) = рм + l р^. Непрерывному варианту (2)
соответствует функция распределения, принимающая для конечных х > 0 значения Р(х) =
= ßp (ВД, а для х = - Р(Х) = +
Итак, функция Р{х) является здесь первообразной для плотности р(х).
Подчеркнем, что интегральное представление Стилтьеса отнюдь не сводится к формальному объединению (1) и (2). Оно имеет смысл для произвольной интегрирующей функции ограниченной вариации Р(х), которая в общем случае однозначно разлагается в сумму трех компонент: компоненты скачков, являющейся суммой счетного семейства ступенчатых функций, "гладкой" абсолютно непрерывной компоненты и непрерывной, но недифференцируемой сингулярной компоненты. Функция распределения Р(х) задана на расширенной положительной полуоси, она может иметь разрывы в конечных неотрицательных точках и на бесконечности. Случай скачка Р при х = 0 соответствует вкладу стадии мгновенного спада упругости, а скачок на бесконечности соответствует наличию инертной компоненты. При аккуратной формулировке при наличии у Р(х) упомянутых разрывов в нуле и на бесконечности придется немного уточнить подынтегральные ядра в формулах (4). Это сделано в Приложении. Заметим также, что введение в определение функции распределения условия Р(0) = 0 означает явное предположение о конечности интегрального PC, что равносильно требованию конечности мгновенного модуля упругости g(0) = G(°°).
Особо ценным для наших целей свойством является однозначная и непрерывная обратимость соотношений связи (4) вязкоупругих функций g(t), G(со) с функцией распределения Р(х). В Приложении приведены характеризующие это свойство теоремы. Понятие непрерывности формулируется на языке поточечной (слабой) сходимости монотонных функций. Наглядно слабая сходимость последовательности функций распределения Рп(х) к функции Р(х) означает геометрическую сходимость графиков с дополнительным условием сходимости значений на бесконечности: limP„(°o) = Р(°°). При этом оказывается, что для
п
обеспечения свойства непрерывной обратимости вовсе не надо знать значения g{t), G(co) при всех /,
со, а достаточно отслеживать указанные функции на "представительных" множествах времен и частот, имеющих в простейшем варианте вид
Т = {0}и[г\Г], О. = {oo}u[co\ö>"],
где ?", со', со" - границы "экспериментальных" интервалов. Особую роль в обращении (4) играет мгновенный модуль упругости, т.е. g(0) = G(°°).
Свойства единственности и непрерывности в сущности означают, что уравнения (4) порождают корректно поставленные обратные задачи на нахождение функции распределения Р(т) по данным g(t) или G(co). Но следует специально оговорить, что обратимая непрерывная связь Р(х) g(t), G(co) гарантирована не для произвольных мыслимых функций g(t), G(со), а лишь для тех, для которых уравнения (4) разрешимы в классе монотонных Р(х). Именно из-за "внешнего" ограничения разрешимости уравнения, связывающие Р(т) с g(t) и G(cо), не следует интерпретировать как традиционные линейные операторные уравнения первого рода. Эти стеснительные ограничения являются источником технических трудностей при решении, но именно они же обеспечивают корректность соответствующих обратных задач!
Чтобы получить способ устойчивого восстановления Р{х) по приближенным g и G, не связанный условием строгой разрешимости соотношений (4), введем естественное (и общепринятое) понятие обобщенного решения как решения задачи минимизации функционала ошибок. В простейшем случае квадратичного функционала эта процедура состоит в решении задач
Jv(i) jexp(-i/x)</P(x)-g(i)
dt
min
(5)
JV(w)J
fCOX
1 + /сох
dP(T)-G( со)
0
d<a
min, (6)
где v(t), д((0) - некоторые неотрицательные весовые функции, определяющие функционал ошибок, а Г, й - области определения экспериментально измеренных g(t), G(со). Минимум ищется в классе всех функций распределения. Решение Р(х) задач (5) и (6), коль скоро оно существует, называется оптимальным интегральным представлением соответствующего кинетического или ча-
стотного модуля. При постановке задач (5), (6) можно без заметного усложнения задать разнообразные ограничивающие условия, учитывающие априорную информацию об искомом PC. Например, указать область возможного сосредоточения PC, исключить мгновенный скачок упругости, задать мажоранту PC и т.п. Выбор весовых функций вполне свободен и определяется априорным знанием о распределении ошибок в экспериментальных значениях. В частности, постоянные v(f), |i(co) разумно применять для относительно узких (две-три декады) диапазонов t, со. Выбор в качестве весов v(/) = г1, ц(со) = со-1 дает равномерный в логарифмическом масштабе функционал ошибок, который естественно использовать для обработки данных с широким диапазоном изменения t, со (более четырех декад). Для ОИГТ доказаны свойства существования, единственности и устойчивости. При этом оказалось, что запас устойчивости неожиданно велик: ОИП устойчивы как к интегрально малым погрешностям функций g(t), G(со), так и к белому шуму конечной амплитуды, наложенному на них. Все эти результаты справедливы и для более общего, чем в задаче (6), комбинированного функционала ошибок, позволяющего восстанавливать Р(х) по вещественной и мнимой части G(со) (т.е. модулю сохранения и потери) или их сочетанию.
При практической численной реализации метода имеем конечное множество времен и частот измерения {г;}, {со¡},j = 1, ...,п функций g(t), G(со). Введя фиксированную сеть узлов {хг}, / = 1, ..., m, покрывающую область возможного сосредоточения PC, дискретизируем непрерывные задачи (5) и (6). При этом искомая функция распределения аппроксимируется ступенчатой комбинацией
Р(х) = ,Л(Х-Х(),
где р = {р,), / = 1,..., т есть набор неизвестных неотрицательных коэффициентов. После применения к внешнему интегралу выражений (5), (6) некоторой линейной квадратурной формулы получим в качестве дискретных аналогов
7 = 1
V
£ехр(-г,/х,)р*-,?(0)
\к :
min, р > 0 (7)
1>
j = 1
Si
к= 1
+ * со ,-т
чч
ç>k-G((Oj)
min, р > 0 (8)
Здесь Vj, Д, - некоторый набор положительных весовых коэффициентов, определяемых совокупностями узлов {¿у}, {со,} и использованной квадратурной схемой. Например, постоянные весовые коэффициенты получатся для формулы прямоугольников и равномерно распределенных узлов при v(0 = 1, |х(со) = 1, а также для логарифмически равномерно распределенных узлов при v(t) = t'1, (i(co) = о)-1. Задачи (7), (8) относятся к проблемам квадратичного программирования. Для их решения известен ряд эффективных конечно-шаговых алгоритмов, доступных в стандартных пакетах математических программ (Maple, Matlab) и фортран-библиотеках (IMSL, CXML), в их числе алгоритм "неотрицательного метода наименьших квадратов" NNLS [13]. Определенные трудности при численном решении могут возникнуть из-за весьма высокой размерности (7), (8) (m, п несколько сотен), поскольку для достижения хорошей аппроксимации исходных непрерывных задач (5), (6) дискретными аналогами необходимо располагать примерно 10-30 узлов tp со,, хк на декаду соответствующего диапазона.
Практика решения многочисленных задач на восстановление PC из модельных и реальных данных релаксации показала, что способ размещения экспериментальных и расчетных узлов сам по себе никак не влияет на качество решения, существенна лишь их густота. При различных способах задания узлов вычисляется различный набор весов, но получаемая в итоге ступенчатая функция распределения определяется с точностью, зависящей от густоты сеток узлов и величины ошибки в g(t), G(со). Не имеет значения возможная при m > п неединственность решения (недо-определенный случай), поскольку любой набор параметров р, реализующий минимум (7), (8), дает равноценную аппроксимацию искомого распределения PC. Введение густой избыточной сети х является даже предпочтительным практическим приемом! Эти практические наблюдения согласуются с изложенной в Приложении теорией корректности ОИП. Более подробно особенности численного решения экстремальных задач минимизации функционала ошибок квадратичного и более общего "гельдерова" типа обсуждаются также в Приложении.
Поскольку свойство устойчивости восстанавливаемых по описанному методу функций распределения PC имеет место только в специфическом "слабом" смысле, уместно пояснить те характеристики PC, которые согласованы с такого рода устойчивостью. По определению слабой непрерывности, устойчивыми являются линейные интегральные средние функции распределения,
т.е. выражения f{x)dP{x), где/(х) - некоторая
непрерывная осредняющая функция. Частный случай осреднения - моменты, для которых /(х) = х". Устойчив, как сказано выше, график функции распределения, причем как геометрическая фигура "в целом". Напротив, локальные числовые характеристики, в частности плотность р(х) = = dP(x)/dx, неустойчивы. Не имеет "слабо устойчивой" трактовки традиционная практика (например, работа [3]) конвертирования дискретного PC в функцию плотности путем прямой интерполяции гребенчатого графика набора 8-функций. Для корректного сглаживания произвольных и в частности, дискретных PC следует привлекать интегральное осреднение по тому или иному ядру. Так как структурные особенности PC зачастую распределены достаточно равномерно в логарифмическом масштабе по х, разумно вычислять сглаженную плотность на основе операции
мультипликативной свертки р£(х) = с£ (x/E)dP(
где - произвольное сглаживающее ядро,
удовлетворяющее условиям а£ > 0, ое (Qft/J£> = 1,
ое(^) = 0 при - 1| > е. Например, удобно для легкости вычислений взять степенное ядро сЕ(^) = = СЛ(% - 1 + е)+(1 + е - причем п > 0 определяет степень гладкости, а Спе - нормирующая константа. Функция плотности рЕ в той же степени гладкая, что и соответствующее ядро. При £ —► 0 непрерывные распределения Ре(х) =
= слабо сходятся к Р(х) и потому по-
строенные сглаженные распределения сколь угодно точно приближаются к исходному дискретному прототипу.
ПРИМЕРЫ ВОССТАНОВЛЕНИЯ РС
Ниже приведены примеры восстановления РС для модельных и реальных экспериментов. Всю-
ду в разделе значения модулей g, G\ G", времени t, частоты со даны соответственно в Па, с, рад/с.
Дискретный спектр конечной неоднородной цепочки Рауза
Рассмотрим упрощенную модель сокращения макромолекулы из ^-звеньев [14]. Считая, что молекула расположена на прямой, обозначим координаты центров ее звеньев, пронумерованных по возрастанию, через ..., хр. При естественных предположениях (высокая вязкость среды, гуковское взаимодействие ближайших соседей с единым для всей цепи коэффициентом упругости, нулевая масса звеньев) изменение координат опишется системой дифференциальных уравнений
= х2 — х} ^
Spxp = хр_х-хр, sixi = xi_l-2xi + xj + l (1 <i<p)
(s, - положительные константы вязкого сопротивления).
Вычислим функцию изменения во времени расстояния между крайними звеньями g(t) = xp(t) -- Xj(t) в виде разложения по экспонентам. Запишем для этого систему (9) в матричной форме х = S~'Ax, где S = diag{s, }f= ,, причем символ diag() обозначает квадратную диагональную матрицу, заданную соответствующим вектором, а А = = 1 -квадратная трехдиагональная симметричная матрица: а„ + [ = 1 ,ап = арр = 1, аи = -2 при 1 < i <р. Решение линейной системы с постоянными коэффициентами (9) выражается через матричный экспоненциал
x(t) = ехр(г5_1Л)л:(0) (10)
Опишем вычисление этого экспоненциала. Пусть D = diag{ A/T!}f= |. Приведем симметричную матрицу D~lAD~l к диагональной форме вещественной ортогональной матрицей U е Rpxp, т.е.
D~lAD~l = UTAU, (И)
где А = diag{X,jf= ,, а - набор собственных чисел D-'AD"1. Тогда Л = DUTAUD и так как exp(fA) = diag{f^,}f= ,,
exp(/S_1A) = D'^diag^Jf^i/D (12)
Для численной реализации матричной факторизации (11) используем быстро сходящийся QR-алгоритм со сдвигами [13]. Поскольку А является отрицательно полуопределенной матрицей ранга р - 1, одно из чисел Xt нулевое, а прочие отрицательны. Переставим спектр и столбцы U так, чтобы Хр = 0. Формулы (10) и (12) приводят нас к искомому экспоненциальному разложению р
8(0 = Ip4expHXt)
k = l (13)
р
Рk = ^ xj(®)djUkj(ukp/dp - ukl/dl), j = 1
причем так как limg(r) = рм = 0, последнее слага-
( —>
емое в выражении (13) нулевое, поэтому
р-1
80) = 5>ехр(-г/тк). хк = -1/Хк (14) к= 1
Если начальные координаты звеньев растут с постоянным шагом (эквидистантны), то коэффициенты разложения (14) неотрицательны, но для иных начальных условий они произвольного знака. "Времена релаксации" в разложении (14) -числа хк = -1/Хк, положительность их вытекает из отрицательности Хк (к < р).
Рассмотрим модельный пример восстановления спектра релаксации по функции релаксации длины эквидистантной цепочки для р = 6, содержащей неоднородные звенья. Зададим набор констант сопротивления и начальных условий: si = i, хi (0) = (г - 1 )/(р - 1), г = 1,..., р. При подготовке модельных "экспериментальных" данных после численной диагонализации матрицы £)_|А£Н была рассчитана функция длины g(t), а по формуле (1) ее частотный модуль G(co). Далее в расчетных узлах i, со наложен шум размаха £ = 10~2 от максимума соответствующей амплитуды, сформированный программой - генератором равномерно распределенных случайных чисел. Вычисленные коэффициенты р(, т,, а также соответствующие
моменты = dP(X) для к = 0, 1, 2 приведены в
табл. 1.
Таблица 1 (пример 1). Дискретный РС и его моменты в модифицированной модели Рауза: точный РС и результаты восстановления при е = 10~2
Точный PC PC восстановлен no g(t) PC восстановлен no G(co)
Pi P, Pi
0.5586 0.2287 х 10"1 0.4281 x 10"1 0.5865 x 10"3 0.58131 0.26375 x 10"'
0.1068 х 101 0.1409 х 10"1 0.9439 0.1910 x 10-1 0.59489 0.14317 x 10"2
0.1671 х 101 0.3533 x 10"1 0.9660 0.4656 x 10"1 0.16426 x 101 0.13763 x 10"1
0.3165 х 101 0.9069 x 10"2 0.1760 x 101 0.1874 x 10"2 0.16810 x 101 0.38617 x 10"1
0.1087 х 102 0.4186 0.1041 x 102 0.4172 0.10655 x 102 0.10467 x 10'1
0.1066 x 102 0.01067 x 10"1 0.10904 x 102 0.40859
0.1000 x 103 0.3831 x 10"2 0.10000 x 103 0.11762 x 10'2
|io = 0.5 , = 4.666, ц2 = 49.67 До = 0.4998, m = 4.907, p2 = 84.82 Ho = 0.5004 , щ = 4.788, (X2 = 61.69
Результаты восстановления из g(t) априорно известного спектра даны вместе с моментами в табл. 1, а также изображены на рис. 1. При имитации эксперимента использовано п = 300 логарифмически равномерных времен и частот в сегментах t € [Ю-4, 103], со е [10-1, 105], а восстановление спектра производилось по схеме (7) на
логарифмически равномерной сетке из т = 400 узлов в сегменте х е [Ю-2,102]. На рис. 1а, 16 видно, что "экспериментальные" значения g(t) приближены с точностью до шума. Из рис. 1в следует, что точный дискретный спектр "рассыпался" на группы близких по х интенсивностей, но, как видно из рис. 1г, ступенчатая первообразная кри-
Па 0.61-
0.4 0.2 0
р, Па 10°Ь
10"1
Ю-2
10
гЗ
(а)
10°
(в) <6
10°
ю2
105 t, с
+ I
о 2
Относительная ошибка (б)
0.01
-0.01
Ч. ; •• •У. '
*.* » *» ♦
- . . •» V. • '
Л Па
0.3
0.1
104 х, с
10° (г)
105 и с
-
- — 1 — 2
_
10°
ю2
104 х, с
Рис. 1. Восстановление PC по кинетической кривой g(t) в примере 1 : а - "экспериментальные" значения (1) и расчетная кривая (2); б - относительная ошибка аппроксимации g(t); в - дискретный PC: рассчитанный PC (Г), исходный PC (2); г - функции распределения PC: рассчитанный (7), исходный (2).
в, Па
0.6 ь
0.4 0.2
0
р, Па 10° I-
10-1
10"
(а)
10° (в)
<Ь
5
сн-
+ .
10°
Относительная ошибка
0.01
-0.01
105 ш, рад/с
+ 1 о 2
Р, Па 0.81-
0.6
0.4
0.2
102
ю4
х, с
(б)
10° (г)
-з£
ю5
со, рад/с
— 1
— 2
10°
10
ю4
X, с
Рис. 2. Восстановление РС по комплексному частотному модулю со) в примере 1: а — "экспериментальные" значения С(со), С'(со) (1) и соответствующие расчетные кривые (2,3); б - относительная ошибка аппроксимации С(со); в - дискретный РС: рассчитанный РС (7), исходный РС (2); г - функции распределения РС: рассчитанный (1), исходный (2).
вая восстановленного распределения близка к истинной. Относительные ошибки в определении моментов |10, р2 равны соответственно 0.0004, 0.0514,0.7070.
В той же табл. 1 и на рис. 2 приведены результаты восстановления по схеме (8) для соответствующего "частотного" эксперимента. Результат этого восстановления в смысле близости графиков функций распределения и их моментов подобен предыдущему. Относительные ошибки в определении моментов Цо, (д.[, ц2 составляют 0.0008, 0.0259, 0.2415.
В табл. 2 приведены также результаты восстановления для того же модельного эксперимента, но с уровнем шума £ = Ю-3. При таком уровне шума набор относительных ошибок в моментах есть 0,0.00042,0.0006 для "кинетического" восстановления и 0, 0.0002,0.0042 для "частотного". Согласие моментов здесь идеально.
Заметим, что поскольку восстановление делалось при т> п, целевой функционал задач минимизации вырожден, и решение не является единственным. Но это не имеет заметного влияния на качество восстановления — все минимизирующие наборы р равно приемлемы!
Непрерывный ВБ\У спектр
Рассмотрим модельный пример восстановления РС по частотному модулю, сгенерированному по непрерывному спектру типа ВБХУ [3].
2
р(т) = Х^'х^^х), (15)
I = 1
где хь х2, х3, п|, п2, с,, с2-параметры, удовлетворяющие соотношениям 0<т3<х2<х1,0<гс1< 1, 1 < п2 < 0, а %[а,ь](т) есть характеристическая функция сегмента [а, Ь], т.е. функция, равная единице при а < х < Ь и нулю в противном случае.
Таблица 2 (пример 1). Дискретный РС и его моменты в модифицированной модели Рауза: результаты восстановления при е = 10
PC восстановлен по g(t) PC восстановлен по G(co)
t, P, P,
0.1000 х 10"1 0.8093 x 10"4 0.60878 0.16995 x 10"1
0.5813 0.7610 x 10"2 0.62300 0.12678 x 10"1
0.5949 0.1872 x 10"1 0.15685 x 101 0.63100 x ÎO'3
0.1464 х 101 0.8293 x 10"2 0.16051 x 101 0.39410 x 10-1
0.1498 x 101 0.3593 x 10-1 0.25469 x 101 0.96683 x 10"2
0.3360 x 101 0.1036 x 10"1 0.26064 x 101 0.96145 x 10"3
0.3438 x 101 0.7915 x 10"3 0.10655 x 102 0.82258 x 10"1
0.1066 x 102 0.5521 x 10-1 0.10904 x 102 0.33736
0.1090 x 102 0.3630 0.89099 x 102 0.76454 x 10"6
0.91180 x 102 0.30438 x 10"4
Цо = 0.5 , щ = 4.665, |x2 = 49.66 Ho = 0.5, (i! = 4.668, ц2 = 49.9
Слагаемые /=1,2 образуют по терминологии [3] стеклообразную и высокоэластическую компоненты PC. В отличие от работы [3] в формулу (15) введен дополнительно положительный параметр х3, чтобы обеспечить конечность мгновенного модуля упругости, т.е. g(0).
Пусть Xj = 0.08, х2 = 0.0094, х3 = 4 х 10"7, щ = = 0.188, п2 = -0.695, су = 13726, с2 = 221.3932. Коэффициенты си с2 выбраны так, чтобы примерно уравнять мгновенные модули упругости вязко-упругой и стеклообразной компонент спектра (15). Подобно предыдущему примеру по формулам (2), (15) были генерированы массивы значений функции g(t), G(со) для логарифмически равномерных сеток {i,}, {со,} при п = 300 на сегментах t € [ 10~W, 102], со е [10°, Ю10] и на них наложен шум размаха е = 0.01 от максимальной амплитуды. Результаты восстановления для m = 600 на логарифмически равномерной сетке в сегменте х е [Ю-8, Ю-1] по схемам (7) и (8) представлены в табл. 3.
Сопоставление моментов восстановленных дискретных PC с априори известными для спектра (15) значениями = 0.8838 х 107, ц, = 0.697 х 103, |х2 = 0.2513 х 102 показывает, что относительные ошибки составляют соответственно 0.0017,0.146, 0.293 для кинетического и 0.00079,0.282,0.487 для частотного варианта восстановления. Графическое изображение результатов восстановления дано в окнах рис. 3 и рис. 4. Видно, что "экспериментальные" значения с точностью до шума аппроксимированы рассчитанными g(t), G(со). На рис. 36
и рис. 46 наложены изображения исходного В8\У и рассчитанного дискретного спектра, а на рис. Зв и рис. 4в приведены соответствующие функции распределения. Отчетливо видна близость графиков функций распределения, тогда как сопоставлять плотность (14) и вычисленный дискретный спектр можно лишь в смысле осреднения последнего по некоторому локальному сглаживающему ядру.
Охарактеризуем в этом примере влияние на качество восстановления точности данных, а также ширину диапазонов соответствующих "экспериментальных" времен и частот. Так, положив уровень шума е = Ю-4, восстановим при тех же сетках (г}, {со}, {х} ступенчатые функции распределения, графики которых изображены на рис. 5. Графики найденных функций распределения стали заметно ближе к соответствующему точному. Моменты распределений приведены в табл. 3, они имеют относительные ошибки 0, 0.007, 0.028 для кинетического и 0,0.023,0.1356 для частотного варианта. Сузим теперь "экспериментальные" диапазоны времен и частот до г е [Ю-8, 1], со е [1, 108]. На том же рис. 5 видно, что графики функций распределения__изменились незначительно. Что касается моментов, то, как следует из табл. 4, моменты нулевого порядка определены точно, а моменты порядка 1, 2 имеют относительную ошибку 0.0086,0.161 для кинетического варианта восстановления и 0.412, 1.884 для частотного варианта. Большее ухудшение качества частотного восстановления связано, вероятно, с тем, что ядро
Таблица 3 (пример 2). Дискретный РС и его моменты в ВБХУ-модели
РС восстановлен по g(t) РС восстановлен по 6(00)
Р; Р/
0.2124 х 10"7 0.2492 х 105 0.4689 х 10"6 0.1258 х 107
0.5661 х 10-6 0.1204 х 106 0.6305 х 10"6 0.1929 х 107
0.5816 х Ю-6 0.3849 х 107 0.6477 х 10"* 0.7023 х 106
0.1617 х 10"5 0.1023 х 106 0.1574 х 10"5 0.2281 х 107
0.1661 х 10"5 0.2544 х Ю7 0.1617 х 10"5 0.4461 х 106
0.4495 х 10"5 0.1101 х 107 0.5727 х Ю-5 0.1270 х 107
0.1726 х 10"4 0.1955 х 106 0.5883 х 10"5 0.6663 х 105
0.1773 х 10"4 0.5629 х 106 0.2029 х 10"4 0.3775 х 106
0.1047 х 10"3 0.2598 х 105 0.2084 х 10"4 0.9262 х 105
0.1076 х 10-3 0.2089 х 106 0.7582 х 10"4 0.2302 х 106
0.2614 х 10"3 0.3144 х 105 0.7789 х 10"4 0.4120 х 105
0.9513 х 10"3 0.5398 х 105 0.3611 х 10"3 0.8709 х 105
0.2675 х 10'1 0.2460 х 105 0.3709 х 10"3 0.8280 х 104
0.4179 х 10~2 0.3696 х 105
0.5838 х 10"1 0.1072 х 105
Мо = 0.8845 х 107, щ = 0.7679 х 103, р2 = 0.1766 х 102 Но = 0.8837 х 107, (X! = 0.8! 92 х 103, р2 = 0.3719 х 102
/сот/(1 + /сот) выходит на стационар с линейной скоростью по со, тогда как "кинетическое" ядро ехр(-г/т) спадает по г экспоненциально. Поэтому при равной точности данных для хорошего воспроизведения моментов запас ширины диапазона частот по сравнению с шириной зоны сосредоточения РС должен быть более значительным (3-4 порядка), в то время как для кинетического вос-
становления достаточно двух порядков. Уменьшение ширины "экспериментального" диапазона можно, согласно изложенной теории, компенсировать ростом числа узлов п и точности данных. В нашем примере при е = 10-4 и указанных "экспериментальных" диапазонах по г и со погрешность в восстановленных моментах р2 будет ~1% для логарифмически равномерной сетки {со} при п = 10000.
Рис. 3. Восстановление РС по кинетической кривой g{t) в примере 2: а — "экспериментальные" значения (/) и расчетная кривая (2); б - рассчитанный дискретный РС (/), плотность р(т) исходного непрерывного ВБ\V-cneKTpa (2); в - функции распределения РС: ступенчатая функция распределения рассчитанного дискретного РС (/), первообразная исходного непрерывного В8\¥-спектра (2).
со, рад/с т, с т, с
Рис. 4. Восстановление РС по комплексному частотному модулю С(м) в примере 2: а - "экспериментальные" значения С'(со), С'(со) (/) и соответствующие расчетные кривые (2,3); б - рассчитанный дискретный РС (/), плотность р(т) исходного непрерывного В8\У-спектра (2); в - функции распределения РС: ступенчатая функция распределения рассчитанного дискретного РС (/), первообразная исходного непрерывного В8\У-спектра (2).
Рис. 5. Результаты восстановления РС для шума е = Ю-4 и различных диапазонов "экспериментальных" времен/частот в примере 2: а - восстановленные дискретные РС для диапазонов г е [Ю-10, 102], со е [10°, Ю10]: "кинетическое" восстановление (У), "частотное" восстановление (2), плотность р(т) исходного ВБХУ-спектра (3); б - функции распределения, соответствующие рисунку а: "кинетическое" восстановление (1), "частотное" восстановление (2), первообразная В^Ш-спектра (3); в - функции распределения для суженных "экспериментальных" диапазонов г е [Ю-8,10°],сое [10°, 108]: "кинетическое" восстановление (/), "частотное" восстановление (2), первообразная В8\¥-спектра (3).
Таблица 4 (пример 2). Моменты восстановленных дискретных РС при е = 10"4 в ВБХУ-модели для различных "экспериментальных" диапазонов времен и частот
Восстановление Ш) Ri
"Кинетическое" для г е [Ю-10 102] 0.8838 х 107 0.6921 х 103 0.2442 х 102
"Частотное" для со € [1,10ю] 0.8838 х 107 0.6808 х 103 0.2172 х 102
"Кинетическое" для г е [Ю-8, 1] 0.8838 х Ю7 0.7030 х 103 0.2918 х 102
"Частотное"для®е [1, 108] 0.8838 х 107 0.9843 х 103 0.7249 х 102
Релаксационный спектр полистирола
Рассмотрим восстановление РС по частотному модулю релаксации напряжения из серии экспериментов, приведенной в работе [3]. Конкретно мы возьмем данные для монодисперсного ПС с М = 3 х 106. Комплексный частотный модуль С(оУ) задан таблично своими действительной и мнимой составляющими С'(со), С"(со). В данном примере мы имеем массив из п = 44 комплексных данных и
используем при расчете по схеме (8) логарифмически равномерную сетку с т = 400 на сегменте те [10"5, 104].
Дискретный спектр, восстановленный по методу ОИП, и его моменты представлены в табл. 5 и на рис. 6. Рисунок 66 показывает, что экспериментальные данные приближены посредством рассчитанной после восстановления теоретической кривой С?(со) с квадратичной ошибкой -1%.
Рис. 6. Восстановление РС монодисперсного полистирола с М = 3 х 106 по комплексному частотному модулю С(со) на основе данных [3]: а — экспериментальные значения С(со), С'(оз) (!) и соответствующие расчетные кривые (2,5); б - относительная ошибка аппроксимации С(со): ошибка ОИП-восстановления (/), ошибка тК-восстановления (2); в - восстановленный дискретный РС: ОИП-спектр (1), ШК-спектр (2); г — функции распределения восстановленных дискретных РС: ОИП-спектр (1), ИИБ-спектр (2).
Таблица 5 (пример 3). Восстановленный по частотному модулю РС полистирола и его моменты (результаты восстановления методом ОИП для двух сегментов времен релаксации т и восстановления методом регуляризации IRIS из [3])
ОИП; сегмент восстановления
[10~5 Ю4] [io-7 105]
т,- Р. Р/ t, 9i
0.10000 х 10"4 0.12287 х 107 0.10000 x 10"6 0.97748 x 108 2.68 x 1(Г5 1.01 x 106
0.10155 XlO"3 0.40375 х 105 0.96600 x 10^ 0.24927 x 106 3.74 x lO"4 9.52 x 104
0.10512 х 10"3 0.16569 х 106 0.30606 x 10-3 0.25727 x 105 2.26 x 10"3 2.53 x 104
0.29679 х Ю-3 0.16334 х 105 0.32051 x 10"3 0.27750 x 105 1.9 x 10"2 9.24 x 103
0.30724 х Ю-3 0.37448 х 105 0.10155 x 10"2 0.32273 x 105 1.29 x lO"1 8.64 x 103
0.99616 х 10-3 0.19573 х 105 0.44437 x 10"2 0.32650 x 104 7.32 x 10"1 1.04 x 104
0.10312 х 10"2 0.13401 х 105 0.12838 x 10"1 0.11399 x 104 4.21 1.6 x 104
0.45649 х 10"2 0.32123 х 104 0.13445 x 10"1 0.80148 x 104 2.38 x 101 2.44 x 104
0.13342 х 10~1 0.91125 х 104 0.16231 0.12164 x 104 1.36 x 102 3.33 x 104
0.16674 0.64814 х 104 0.16997 0.78375 x 104 7.37 x 102 4.34 x 104
0.17261 0.25728 х 104 0.15559 x 101 0.95812 x 104 2.66 x 103 5.33 x 104
0.15263 х 101 0.31653 х 104 0.29679 x 101 0.48952 x 104 7.79 x 103 1.99 x 104
0.15800 х 101 0.65448 х 104 0.78191 x 101 0.11375 x 105
0.29452 х 101 0.46268 х 104 0.31200 x 102 0.22423 x 105
0.77592 х 101 0.11539 xlO5 0.15679 x 103 0.31856 x 105
0.30961 х 102 0.13634 х 105 0.16419 x 103 0.17245 x 104
0.32051 х 102 0.89417 х 104 0.75242 x 103 0.32168 x 105
0.15740 х 103 0.20915 х 105 0.78794 x 103 0.10883 x 105
0.16294 х 103 0.13051 х 105 0.26143 x 103 0.47128 x 105
0.77294 х 103 0.57850 х 104 0.62805 x 104 0.20985 x 105
0.80015 х 103 0.39656 х 105 0.65770 x 104 0.51079 x 104
0.28781 х 104 0.14607 xlO5 0.10000 x 106 0.11244 x 104
0.29794 х 104 0.34765 х 105
0.46714 х 104 0.10166 х 105
0.10000 хЮ5 0.11862 x 105
Цо = 0.1742 х 107, и, = 0.3542 x 109, Lu = 0.983 x 108, x. = 0.4399 x 109, Po = 0.13491 x 107, p! = 0.33397 x 109,
р2 = 0.1867 х 1013 p2 = 0.1264 x 1014 |i2 = 0.16089 x 1013
В табл. 5 дан также дискретный РС, восстановленный по программе IRIS (одна из реализаций метода регуляризации), приведенный в работе [3], и его моменты, а на рис. 66,6в, 6г произведено наложение графиков ОИП- и IRIS-восстановле-ний. Итак, моменты порядка 0, 1, 2 различаются на 22, 5.7, 13.8% при восстановлении РС на сегменте х е [10~5, 104]. Из рис. 66 видно, что хотя ошибка аппроксимации частотных модулей в обоих восстановленных РС близка по порядку, но ОИП-уклонение "случайно", а IRIS-уклонение явно имеет систематические "волны". Что касается сопоставления функций распределения, то из
рис. 6г отчетливо видна параллельность их графиков в зоне х е [Ю-4, 101]. Параллельный сдвиг функций распределения из-за начальных скачков не является в этом случае существенным расхождением восстановления - он просто вызван недостаточной шириной диапазона данных в области высоких частот. Так, на рис. 4а видно, что модули не выходят в диапазоне частот эксперимента на стационар, причем в зоне о) > 104 начинает проявляться стеклообразная составляющая спектра. Для ее достаточно точного воспроизведения необходимо иметь, вероятно, либо значительно более прецизионные данные в имеющемся диапазо-
не частот (ориентировочно точность ~10"5), либо данные с прежней точностью, но в диапазоне до 0) ~ 108. В действительности, исходя из приведенного в работе [3] неполного и зашумленного массива данных С(со), а также не располагая априорной информацией о месте сосредоточения искомого РС, вообще нельзя достоверно судить о значениях его моментов. Так, проведя по схеме (8) ОИП-восстановление на сегменте х е [Ю-7, 105], мы получим РС, также приведенный вместе с моментами в табл. 5. Найденные моменты различаются по сравнению с соответствующими ШК-моментами уже на 7400, 31, 685%! И при этом в среднем диапазоне 10-4 < х < 101 все восстановленные ОИП-распределения почти идентичны.
Таким образом, в среднем диапазоне ОИП- и ШК-восстановления РС близки на языке функций распределения, хотя параметры индивидуальных линий точечных спектров существенно различны. Близость графиков функций распределения объяснима в рамках теории устойчивости интегральных представлений, поскольку авторы работы [3] предъявили также некоторый неотрицательный РС, аппроксимирующий с шумовой точностью экспериментальные данные. Заметим также, что в работе [3] дискретный спектр в конечном счете конвертирован путем интерполяции весов 5 функций в гладкий кусочно-степенной BSW-cпeктp, и эта операция сглаживания несколько ухудшила исходную аппроксимацию.
СРАВНЕНИЕ С ЛИТЕРАТУРНЫМИ ДАННЫМИ
К настоящему времени сформировалась обширная литература, посвященная обратной задаче восстановления РС, рассмотренной в настоящей статье. Заметим, что подобные постановки возникают в различных областях естествознания, не обязательно связанных с реологией полимеров, на которую ориентировано наше изложение. Подробный обзор методов решения содержит работа [15], причем здесь задача трактуется как проблема обращения преобразования Лапласа и ассоциируется с лазерными корреляционными экспериментами. Создано и эксплуатируется множество компьютерных программ, например СОШШ [16], ЕРЯ [17], ОЕЫЕКЕС [18], ЕТ1К1ШЗ [19], в том числе встроенных в стандартное математическое обеспечение серийных спектрометров. Как правило, их алгоритмы основаны на ко-
нечномерной аппроксимации - подгонке под эксперимент спектра заранее предписанной формы (полуобратная задача) или на процессах регуляризации того или иного типа, устраняющих некорректность путем внесения в задачу искусственного возмущения. К этому кругу идей примыкают и комбинированные алгоритмы с максимизацией энтропии РС, а также с предварительным сингулярным разложением матрицы линейной математической модели.
Каково же соотношение представленной работы с известными?
Теория ОИП Стилтьеса как раз и объясняет, почему, когда и в каком смысле непохожие на первый взгляд результаты различных восстановлений РС согласуются по интегральным характеристикам. Согласие возникает тогда и только тогда, когда они достаточно хорошо аппроксимируют данные и удовлетворяют некоторому ограничению принадлежности "выпуклому локально слабо компактному множеству", которое, например, выделяется простейшим условием неотрицательности. В этом случае теория ОИП гарантирует, что независимо от способа получения РС он близок к истине на языке функций распределения.
Отметим еще одно характерное обстоятельство в практике решения подобных задач. Начиная примерно с 80-х годов XX века, практики твердо усвоили, что без наложения дополнительных ограничений типа неотрицательности технология регуляризации интегрального уравнения первого рода в классической манере (например, по Тихонову) бесполезна для большинства реальных обратных задач спектроскопии. Это связано с тем, что в ходе подгонки параметра регуляризации она приводит сразу от сильно осциллирующих к недопустимо уширенным "решениям". Поэтому во многих стандартных пакетах (например, Con-tine) предусмотрена возможность наложения априорного линейного ограничения на решение, выбираемого пользователем, в том числе ограничения неотрицательности. Но минимизируемый функционал (как правило, квадратичный) в любом случае одновременно подвергается процедуре возмущения-регуляризации. Иными словами, регуляризируется изначально корректная (на языке распределений) обратная задача! Вполне естественно, что такая верность идее регуляриза-
ции имеет неприятные методические последствия. В известных пакетах решения обратных задач спектроскопии, и в частности, восстановления РС алгоритмы минимизации квадратичного функционала ошибок изначально привязаны к строго положительному значению параметра регуляризации. При его уменьшении они вырождаются и теряют работоспособность. Таким образом, в рамках стандартных пакетов программ определения РС возможно достичь лишь приближения к минимуму чистых функционалов ошибок и только лишь для квадратичного варианта восстановления.
В демонстрационных примерах статьи также рассмотрен только вариант квадратичной схемы восстановления. Реализация схемы более общего "гельдерова" функционала ошибок может иметь преимущество при обработке обширных массивов прецизионных экспериментальных данных.
Отметим еще одно характерное обстоятельство. В работах [3, 19, 20] обсуждается влияние на РС неполноты экспериментальных данных из-за конечности диапазона частот измерения релаксационных модулей. Такая неполнота рассматривается как непреодолимый барьер для точности методик восстановления. Яснее всего этот барьер виден на процедуре преобразования (конверсии) частотных модулей С(со) = С(со) + ¿С'(со) в кинетическую релаксационную функцию g{t). Например, в работе [19] приведены формулы
s(0 = P- + -J
2 fG'(co) - Р.» . , w
1 sin(coi)fife)
со
(16)
8{t) = р~ + л1
2 rG"(co)
CO
COS((Ot)d(>),
основанные на соотношении (3), а также на обратном преобразовании Фурье [21]. Эти явные формулы устойчивы к интегрально малым в квадратичном смысле ошибкам в выражении (С(со) - р„)/о. Конечно, для применения соотношений (16) нужно еще достаточно точно знать вес инертной составляющей модуля р„, который вносит равноценный с С(со) вклад в ошибку. Но в явных формулах (16) ясно видно существенное прямое влияние диапазона частот. Пусть для просто-
ты рм = 0. В силу теоремы Планшереля [21], справедливо тождество
J rcJ со со
dm
Поэтому при использовании в выражениях (16) конечного частотного диапазона П = [со', со"] интеграл неизбежно искажен на значение выражений
2 тс
' J|G'(co)/co|:
dcо
2тс~'j"|G"(co)/co|2i/co,
Q'
где £2' = (0, со') и (со", - объединение интервалов дополнительных к Q частот. Вклад дополнения "экспериментального" диапазона в ошибку g(t) может быть велик, особенно для сред с выраженной стеклообразной компонентой, и его нельзя устранить, повышая точность измерения G(co) в ограниченном диапазоне частот. В противовес этому ограничению "неполноты данных" Фурье-методики конверсии G(co) —► g(t ) в рамках развиваемого нами ОИП-метода полная для восстановления РС информация содержится в произвольно малой окрестности времен или частот, и неполноту экспериментального диапазона можно компенсировать, по крайней мере в принципе, уменьшением систематической ошибки и шага измерений в любом доступном диапазоне. Информация используется неявно, в силу однозначной связи значений аналитических функций переменных г, со, выраженных интегралами Стилтьеса (4), причем ошибка при ОИП-конверсии стремится к нулю равномерно по t при повышении качества данных, тогда как при Фурье-конверсии (16) точность улучшается всего лишь в среднеквадратическом смысле.
ОБСУЖДЕНИЕ РЕЗУЛЬТАТОВ
Функции g(t), имеющие интегральное представление (4), изучены в литературе: это класс так называемых абсолютно монотонных функций, которые бесконечно дифференцируемы при t > 0 и удовлетворяют при всех целых п> 0 и вещественных t > 0 неравенству (-\ )"g,n>(t) > 0. Соот-
ветствующая формула (4) есть представление Бернштейна [9, 11, 12] абсолютно монотонной функции. Обратно, это условие определяет возможность представления g(t) в виде (4), т.е. является критерием существования и соответственно самой возможности восстановления PC по кинетическим данным. Указанный критерий весьма полезен. Например, отмеченная в работе [24] неудача при определении PC из кинетической функции напряжения "последовательной разгрузки" связана, по-видимому, с несоответствием условий эксперимента и критерия абсолютной монотонности, т.е. в конечном счете с неадекватностью PC-модели релаксации рассматриваемому физическому процессу.
К сказанному надо добавить, что представление PC в виде интегральной функции распределения Р(х) вообще снимает проблему неединственности представления PC. В приведенных в статье примерах, различные, казалось бы, спектры дают практически одну и ту же Pix). По-видимому, это обстоятельство и является причиной того, что эксперимент является нечувствительным к локальному виду PC, как было показано во многих работах и подытожено в работе [6].
То, что моменты PC однозначно характеризуют поведение полимера в терминальной области (вязкость, податливость), хорошо известно [1, 2]. В настоящей работе продемонстрировано, что предлагаемый метод восстановления PC приводит к получению устойчивых значений моментов и вообще любых интегральных характеристик PC. Локальные же характеристики PC, и в частности предельные времена релаксации, неустойчивы и неинформативны. Поэтому факт зависимости максимального значения времени релаксации от способа восстановления PC, на который указано в работе [6], не столь существен, поскольку определяющим релаксационные свойства фактором являются моменты распределения. А они-то обладают высокой степенью устойчивости, как следует из настоящей работы.
Вместе с тем, предлагаемый математический подход предоставляет возможность дискриминировать экспериментальные методы, позволяющие определить PC полимера. Формально задача минимизации функционала ошибок в классе монотонных функций распределения всегда имеет единственное решение, т.е. ОИП. Получаемый в
итоге PC дает, по определению, наилучшее из возможных объяснение данных в рамках линейной модели с неотрицательным PC, имеющим физический смысл и связь со структурой полимера. Но в рамках применимости линейной теории возможны ситуации со знакопеременным разложением вязкоупругих функций по экспоненциальным модам. Признаком реализации этого случая является то, что найденное ОИП не обеспечивает хорошего (с точностью до шума) воспроизведения экспериментальных g(t), G(со). Например, если в первой задаче, рассмотренной выше, начальное состояние цепочки не является эквидистантным, то получить неотрицательный PC, воспроизводящий модельную релаксацию, невозможно, так как некоторые из коэффициентов (14) при экспоненциальных членах отрицательны. Данное обстоятельство указывает на то, что установление PC как характеристики структуры полимера требует совершенно определенной постановки эксперимента [24]. Хорошо известно, что в случае релаксации цепи в режиме деформации и напряжения получается несколько различный PC [14]. Очевидно, это обусловлено тем, что структура полимера проявляется весьма опосредованно в случае перемещения цепи как целого, поэтому PC зависит от характера ее движения.
Заметим, что сам по себе метод ОИП легко обобщить на случай PC произвольного знака. Для этого следует перейти в формулах (4) от интегрирования по монотонной функции распределения Р(х) к интегрированию относительно функции распределения ограниченной вариации и ввести в постановку экстремальных задач (5), (6) априорное предположение принадлежности Р(х) некоторому выпуклому слабо локально компактному множеству - ограничению (М0 в обозначении Приложения). Такое обобщение в статье не сделано лишь потому, что в существующих физических теориях фигурирует именно неотрицательный PC.
ЗАКЛЮЧЕНИЕ
Таким образом, задача восстановления монотонной функции распределения спектра времен релаксации по кинетической кривой или частотному модулю математически корректна. Просто трактовать ее необходимо в рамках метода ОИП, а не традиционно, как задачу подгонки (фиттин-
га) набора параметров точечного спектра или как проблему решения интегрального уравнения первого рода относительно функции плотности непрерывного спектра. Еще раз подчеркнем особенности предлагаемого подхода.
Точность восстановления функции распределения РС и любых ее интегральных средних лимитирована лишь систематической ошибкой и "густотой" экспериментальных данных. Фактором, влияющим на систематическую ошибку, естественно, является и качественная адекватность математической модели, описывающей рассматриваемые релаксационные процессы в терминах экспоненциальных мод Максвелла.
Теоретически методом ОИП возможно восстановление РС с любой точностью по набору кинетических или частотных данных в любых диапазонах. Однако сужение диапазонов требует взамен резкого увеличения точности и числа данных.
Для надежного восстановления "быстрой" части РС при характерной для современной экспериментальной практики точности и полноте данных (и ~ 100, £ ~ 0.01) необходимо, чтобы диапазон эксперимента было шире зоны сосредоточения РС на 1-2 порядка при "кинетическом" восстановлении и на 3-4 порядка - при "частотном". В противном случае функция распределения определяется с точностью до начального скачка, локальная структура которого неинформативна, а также резко падает точность определения моментов.
При решении задачи восстановления РС полезно включать в число его параметров "крайние" мгновенную и инертную компоненты. В этом случае в число данных следует включить значение мгновенного модуля упругости. Но можно также отказаться от восстановления крайних компонент, введя в процедуру восстановления явное априорное ограничение зоны расположения РС либо оценочную мажоранту функции распределения.
Утверждения, характеризующие понятие о РС как формальную фикцию только вследствие неудач его применения для интерпретации экспериментальных данных, в рамках математически некорректного аппарата (при дискретном и не-
прерывном описании РС) являются неосновательными.
Наличие устойчивого математического аппарата ОИП для восстановления РС из динамических экспериментов также не означает, разумеется, что РС есть физическая реальность. Однако предложенный подход делает перспективным и актуальным пересмотр (уже на основе устойчивой трактовки РС в терминах функций распределения) результатов релаксационных экспериментов с целью выяснения рамок применимости несомненно удобного фундаментального понятия о РС и его взаимосвязи со структурой полимеров.
ПРИЛОЖЕНИЕ
ПРЕДСТАВЛЕНИЕ РЕЛАКСАЦИОННЫХ
ФУНКЦИЙ ИНТЕГРАЛОМ СТИЛТЬЕСА
Для начала уточним способ представления релаксационных функций интегралом Стилтьеса. Условимся [7, 8] описывать РС не плотностью или набором параметров дискретных мод, а функцией распределения Р(т), где Р{х) есть интегральная интенсивность РС на полуинтервале [0, х). Тем самым наше рассмотрение заранее ограничено случаем конечного интегрального спектра, что равносильно также конечности значения мгновенного модуля упругости. Функция Р(т) обладает в силу своего определения следующими свойствами: она определена на расширенной неотрицательной полуоси R+ = [0, °о], непрерывна слева при 0 < х < «> и удовлетворяет равенству Р(0) = 0. Обозначим класс всех таких монотонных функций через M(R+). Мы можем выразить функции g(t), G(co) при помощи интеграла в смысле Стилтьеса относительно Р(т) соотношениями (4). В формулах (4) на самом деле приходится интегрировать функции непрерывные внутри R+, но имеющие концевые разрывы первого рода. Поэтому под интегралом понимается не стандартный интеграл Римана-Стилтьеса, определенный для непрерывных функций, а интеграл Лебега-Стилтьеса. Впрочем, для наших целей его можно эквивалентно описать в рамках специально модифицированной конструкции Римана-Стилтьеса. А именно, под интегралом функции/(х) по функ-
Таблица 6. Определение ядер a (i, х) и Ь (со, т)
X а (t, х) b (со, х)
г = 0 0 < i < °° f=oo со = 0 0 < (0 < оо (0 = °°
х = 0 1 0 0 0 0 1
0<Т<оо 1 exp (-i/x) 0 0 /сот/(1 + im) 1
X = ОО 1 1 1 1 1 1
ции Р(х), имеющей ограниченную вариацию на полуоси R+, понимается предел при е —► 0 сумм
S(îlÊ) = /(0)Р(0+)+/(оо)(Р(со)-Р(оо-)) +
"E-1 (П.1)
I = 1
построенных для параметризованной посредством параметра £ > 0 сети пар цепочек чисел л£ =
ле п - 1
= ({<}"= 1, ) таких, что 0 < X! < < х2
1 < < \ < причем lirnx, = О,
£¿0
limxe£ =оо, lim шах (Х;+1 - х^ ) = 0.
eiO «-1 ^0iSiS)l«-2
Этот предел (т.е. интеграл) существует для ограниченных на R+ функций/(х), являющихся непрерывными на интервале 0 < х < °° и имеющими односторонние пределы на концах 0,
Подынтегральные ядра (4) не определены для всех сочетаний конечных и бесконечных неотрицательных значений t, х, со. Поскольку эти особенности вносят ненулевой вклад в интегралы (4) при наличии разрывов Р на при х = 0, <», произведем аккуратное доопределение exp(-i/x), /сох/(1 + /сох) функциями a(t, х), ¿(со, х), положив их равными пределам соответствующих выражений по конечным положительным х при фиксированных t, со:
a(t*, х*) = lim exp(-r*, х)
T -» T*
¿»(со*, x*) = lim /со*х/(1 + /со*х)
T -> T*
Значения таким образом доопределенных ядер сведены в табл. 6. Введя a(t, х), ¿»(со, х), переопределим формулы (4) соотношениями
оо во
g(t) = Ja(t, x)dP(x), G(со) = Jè(co, x)dP(x), (П.2)
о о
имеющими точный смысл при любых 0 < т, t, со < t». При t = t», со = 0 ядра a(t, х), ¿»(со, х) имеют разрыв первого рода при х = что не препятствует их интегрируемости.
Заметим, что за рамки предположения интегральной конечности выходит модель PC критического геля из работы [20], для которого р(х) = г~"_1/Г(п), 0 < п < 1 и соответственно
р (x)dx = оо, g(t ) = г", G(co) = Г( 1 - n)(i(ü)n. Чтобы
использовать представление Стилтьеса в такой ситуации, следует при определении функции Р(х) вместо краевого условия Р(0) = 0 ввести условие закрепления Р(х0) = 0 в произвольной фиксированной положительной точке х0. При этом для х, < х2 разность Р(х2) - Р(хх) имеет смысл интегрального спектра на полуинтервале [xl7 х2). Представления (4) сохранят силу, но в них вместо определенного интеграла Стилтьеса фигурирует несобственный интеграл [11] на открытой полуоси (0, оо). Мы исключили этот случай из основного рассмотрения, поскольку он требует ряда оговорок, притом что ситуация с интегрально не ограниченным PC достаточно экзотична.
ЕДИНСТВЕННОСТЬ И УСТОЙЧИВОСТЬ ИНТЕГРАЛЬНЫХ ПРЕДСТАВЛЕНИЙ
В приложениях при использовании представлений (4) для восстановления Р(х) неизбежно возникает вопрос о единственности функции распределения, соответствующей данным вязкоупругим функциям g(t), G(со), а также об ее устойчивости к возмущениям последних. В формулировке нижеследующих результатов участвуют некоторые множества неотрицательных чисел Т, £2, которые являются абстрактными прототипами совокупности времен и частот измерения значений g(t), G(со). Будем предполагать, что Т, £2 борелевские, обладают конечной и положительной предельной точкой, и сверх того Т содержит 0, а £2 содержит о«. Следовательно, T, £2 должны быть бесконечны. Примером подходящих под этот критерий множеств являются Т= (0| u [г,, i2], £2 = и [со,, со2], где 0 < г, < г2, 0 < со, < со2 - некоторые числа.
Свойство единственности заключается в следующем.
Предложение 1. Пусть Р\(Х), Р2{т) - две функции распределения. Если при всех t е Т выполне-
но #!</) = g2(t) или при всех со е £2 выполнено G,(co) = G2(со), то Р,(х) = Р2(х) при любом 0 < х <
Свойство устойчивости формулируется в терминах слабой сходимости монотонных функций. Напомним это понятие, лежащее в основании аппарата теории вероятностей.
Определение. Говорят, что последовательность монотонных на отрезке [а, ß] функций Рп слабо сходится к монотонной функции Р (общепринятое обозначение - Рп => Р), если Рп{х) —«-—»- Р(х) при х = сх, ß, а также в точках а < х < ß, в которых функция Р непрерывна.
Известно, что монотонная функция имеет не более чем счетное множество точек разрыва, поэтому слабая сходимость влечет поточечную сходимость на плотном подмножестве отрезка мощности континуум, содержащем его концы. Наглядным геометрическим эквивалентом слабой сходимости является стремление линейных замыканий графиков Рп (замыкание вклеиванием отрезков скачков) к линейному замыканию графика Р. Справедливо следующее свойство, близкое принципу непрерывности характеристических функций в теории вероятностей [11].
Предложение 2. Пусть Р, g, G, Рп, gn, G„ in = 1,. 2, ...) есть последовательность из функций распределения и соответствующих кинетических и частотных функций - модулей. Тогда слабая сходимость Рп => Р на расширенной полупрямой R+ = = [0, оо] равносильна каждому из двух условий: при всех t е Т выполнено gn(t) —- g(t); при всех со € £2 выполнено Gn(со) —- G(co).
Условие сходимости gn(t), Gn{со) при t = 0, со = т.е. сходимость мгновенных модулей упругости, является существенным. Если его отбросить, то из сходимости gn{t) —- g(t), G„(со) —► G(co) при t, со > 0 вытекает лишь, что limP„(x) = Р(х) + с ■ где с - произвольная постоянная.
Сформулированные выше свойства единственности и непрерывности означают, что обратные задачи
A[P]=g, PeM(R+) (П.З)
ß[P] = G, Ре M(R+) (П.4)
являются корректно поставленными в классическом смысле Адамара. Здесь А[Р](0 = \ a(t,
x)dP(x), teT и Я[Р](со) = (со, x)dP(x), со е Q - интегральные операторы, а функции g, G принадлежат метрическим пространствам Zg = A[M(R+)], ZG= B[M(R+)] соответственно. При аккуратной математической формулировке обратных задач (П.З), (П.4) следует еще определить метрики в образе и прообразе операторов А, В. В качестве таковых можно выбрать, например, интегральные метрики
dp(PUP2) = |Р,(«)-Р2(оо)| + + j\Pl(x)-P2(x)\/(l+x1)dx
dg(gltg2) = M0)-*2(0)|+J|*,(f)-S2(O|<fc
т
dc(GltG2) = |G,(oo)-G2(OO)| +
+ JlG^co) - G2(co)|JCO,
a
где dt, d(0 обозначают интегрирование по некоторым конечным положительным мерам Радона на
г, а
Хотя операторы А, В линейны, уравнения (П.З), (П.4) не являются традиционными линейными уравнениями первого рода из-за ограничений выпуклости областей определения и значений. Эти стеснительные ограничения влекут трудности при решении, но именно они же обеспечивают непрерывную обратимость операторов A: M(R+) —► Zg, В: M(R+) —► ZG и, следовательно, корректность соответствующих обратных задач.
Известно [11, 21], что из слабой сходимости функций распределения Рп => Р вытекает сходимость интегралов Стилтьеса limjf{x)dPn(x) =
= |f(x)dP(x) для произвольной непрерывной на
R+ функции /(х). Если в качестве / брать срезанные степенные функции /(х) = min{x", С"}, где С > 0 параметр, то из Предложения 2 следует, что интегральные моменты распределений сходятся, когда поточечно сходятся соответствующие кинетические функции g„(t), а сами распределения Рп(х) сосредотрочены в интервале времен релаксации х < С.
ОПТИМАЛЬНЫЕ ИНТЕГРАЛЬНЫЕ ПРЕДСТАВЛЕНИЯ
При экспериментальном измерении g(t), G(со) получаются функции с шумовой составляющей, заведомо нарушающие тонкие дескриптивные условия разрешимости g е Zg, G е ZG задач (П.З), (П.4). Определим поэтому естественное понятие обобщенного решения, рассмотрев экстремальные задачи минимизации функционалов ошибок. Чтобы ввести эти задачи в достаточно гибкой форме, зададим набор их параметров q, Т, Q., v(dt), H(dco), М0. Здесь q > 1 есть числовой параметр функционала, именуемый показателем Гельдера, Т, Q - введенные в предыдущем разделе множества "экспериментальных" времен и частот, v(dt), д(г/со) - некоторые неотрицательные весовые меры Радона, которые конечны и локально строго положительны на множествах Г, Q, и обладают скачками в точках t = 0, со = (т.е. удовлетворяют неравенствам v({0}) > 0, (!({«>})> 0). Наконец, М0 есть некое выпуклое слабо замкнутое подмножество совокупности M(R+) всех функций распределения. Рассмотрим следующие экстремальные задачи:
|И[Р](0-*(01М<*0 — min
т
при Р е М0
J|fi[P](oo) - G(CÙ)IV(Î/CO) — mi«
т
при Р е М0
(П.6)
cx'J|fi'[P](co) - С(а>)| V(<fa>) + ex" J|S"[P](co) - G"(co)r"|ï'(i/co) — min
(П.7)
при Р € М0, где а', а" > 0, а' + а" = 1 - весовые коэффициенты, А'[Р](со) = Ре ¿(со, т)ЛР(т), й"[Р](со) =
= |^1тЬ(со, хУР(т) - операторы, и вместо одного
комплекта (^, р) из показателя Гельдера и меры фигурируют (</', д'), (д\ ц"). Задача (П.7) обобщает обсуждавшийся в работах [3, 6] взвешенный функционал ошибок и включает как крайние частные случаи восстановление по одной функции С(со) или С'(со). Формулируемые ниже результаты равно справедливы для экстремальных задач (П.5), (П.6), (П.7), если сделать предположение о локальной строгой положительности на О. меры р.(^со) = р'(^со) + р"(^со) и о наличии у р" скачка в со = °о, т.е. ц"({оо})>0.
Для ОИП имеет место следующее свойство единственности и устойчивости
Предложение 3. Для всякой функции g из
(П 5) ^ч существует и единственно распределе-
ние Р минимизирующее (П.5), т.е. ОИП. Если последовательность gn слабо сходится в (Т, V) к g, то соответствующие ОИП слабо сходятся, т.е. Р„ => Р.
Минимум в задачах (П.5), (П.6) ищется среди монотонных функций класса M(RJ, принадлежащих подмножеству М0. Решение, коль скоро оно существует, называется оптимальным интегральным представлением g, G (ОИП) по ядрам a(t, х),
R С
¿(со, х) в банаховых пространствах Lq (Г, v), Lq (Q, р.) соответственно. Впрочем, для (П.6), ввиду того, что образ оператора В[Р] является комплекс-нозначной функцией, можно сформировать и гибридную задачу минимизации взвешенной суммы уклонений вещественной и мнимой части
Предложение 4. Для всякой функции G из Lq (Q, |i) существует и единственно распределение Р минимизирующее (П.6), (П.7), т.е. ОИП. Если последовательность G„ слабо сходится в
Lq (Q, ц) к G, то соответствующие ОИП слабо сходятся, т.е. Рп => Р.
Под слабой сходимостью gn(t) —- g(t), G„(со) —► G(co) понимается, что для всяких функций y(f) € LqKq_l}{T, v), Ч'(со) е 4%.u(Q,p) выполнено соответствующее соотношение сходимости интегралов
limJV(/)&,(f)v(A) = jv(0*(f)v(<fc)
г г
limJvF(co)G„(a>)p.(iito) = jV(co)G(co)n(Ao)
(П.8)
Сходимость по нормам пространств Lq (Т, v),
с
Lq (£2, |i) (т.е. сильная сходимость) влечет слабую, но не наоборот [21].
Смысл Предложений 3 и 4 состоит в том, что ОИП однозначно определены и чрезвычайно устойчивы. А именно, они устойчивы к интегрально малым погрешностям и даже к белому шуму конечной амплитуды, наложенному на истинные g(t), G(cо). Последнее обстоятельство связано с тем, что слабая сходимость (П.8) является адекватной абстракцией последовательности измерений по сгущающейся сети времен или частот, содержащих случайную шумовую составляющую, которая обладает нулевым математическим ожиданием на любом временном или частотном интервале.
При постановке экстремальных задач (П.6), (П.7) можно теоретически использовать любые конечные положительные весовые меры v(dt), ¡¿(г/со) (со скачками в точках t = О, О) = а также гельдеров показатель q, который регулирует чувствительность минимизируемых функционалов ошибок к "выбросам" в экспериментальных g(t), G(со). Эта чувствительность монотонно зависит от q. Множество М0, включенное в постановку задач (П.6), (П.7), позволяет вносить многие типы априорной информации об искомом PC. Отметим, что всегда можно решать эти задачи и для максимально широкого М0 = M(R+). Но задавая, например,
М0 = {Р: Р(х) = 0 при х< а, Р(Х) = const При ß < X < ОО },
(П.9)
мы тем самым потребуем, чтобы искомый PC был сосредоточен на сегменте [а, ß].
Можно указать, что искомый PC интегрально непрерывен при х = 0 с заданной априорной оценкой темпа роста, положив
М0 = [Р: Р(х) < о(х)Р(°°), 0 < х < оо) (П. 10)
для некой неотрицательной оценочной функции с(х) с условием limc(x) = 0. Мотивированное за-
тХО
дание М0 повышает точность восстановления вследствие сужения запаса распределений Р -кандидатов на аппроксимацию экспериментальных данных. Приведенные варианты (П.9) при а > 0 и (П. 10) задания М0 позволяют решать задачу восстановления в классе непрерывных при х = 0 распределений. Это может быть ценно при обработке динамических данных веществ, для которых об отсутствии компоненты мгновенной упругости заранее известно. Например, именно так обстоит дело для каучуков.
Характерным примером недопустимого способа задания ограничения непрерывности PC при х = 0 служит М0 = {Р: Р(0+) = 0}. Дело в том, что хотя М0 состоит из всех непрерывных в 0 распределений, оно не слабо замкнуто, поэтому задачи (П.5), (П.6), вообще говоря, неразрешимы!
ЧИСЛЕННАЯ РЕАЛИЗАЦИЯ
На практике имеем конечное множество времен и частот измерения Т= {i,}, £2 = {со,},/ = 1,..., п функций g(t), G((0). Произведем дискретизацию непрерывных задач, введя фиксированную сеть узлов {х,}, i = 1,..., m. Искомая функция распределения аппроксимируется ступенчатой комбинацией Р(х) = j р¡h (х - Xj), где р = {р,}, i = l,...,
m есть набор неизвестных неотрицательных коэффициентов.
Теперь вместо континуальных задач (П.5), (П.6) получим их дискретные приближенные аналоги
X,)p,.-g(f;) при р е М0
i = 1
min
(П. 11)
j = 1
X,)P1 - G(co;) при P G M0,
1 = 1
min
(П.12)
где {V,}, {ц,} , / = 1, ..., п есть фиксированные наборы фиксированных положительных весовых коэффициентов. Ограничение р е М0 означает здесь принадлежность вектора р некоторому вы-
пуклому замкнутому многогранному множеству,
г»/И
лежащему в конусе К+ неотрицательных т-мер-ных векторов р, заданному, например, конечным числом нестрогих линейных неравенств М0 = {р:
ХГ= 1с«р* - Р< - 1 = 1.....^'= ^ т} для
некоторой матрицы С = {си}\'™= 1 и вектора с1 =
Это задачи минимизации конечномерной строго выпуклой функции на выпуклом многограннике. При их решения эффективен метод логарифмических штрафов [22] в сочетании со спуском второго порядка для последовательности промежуточных безусловных задач минимизации, который приводит к быстро сходящимся итерациям (ориентировочно затрачивается 30 итераций ньютоновского типа для достижения минимума с машинной точностью по целевой функции). В наиболее обычном случае использования показателя ^г = 2 (квадратичный функционал ошибок) (П.11), (П. 12) относятся к квадратичному программированию, и для их решения применимы известные конечно-шаговые алгоритмы прямого, двойственного или комбинированного типа [23], в частности, прямой алгоритм NN1^ [13]. В предельном случае = 1 задачи (П.11), (П. 12) легко преобразовать к форме линейного программирования, поэтому для их решения здесь применим конечный алгоритм симплекс-метода.
Обсудим вкратце особенности (П.11), (П. 12). Прежде всего для достижения хорошей аппроксимации непрерывных задач (П.5), (П.6) дискретными аналогами необходимы достаточно густая сеть узлов т„ покрывающая область возможного сосредоточения релаксационного спектра, ориентировочно, десятки на декаду. Среди входящих в измерения соу желательно (необязательно) иметь £ = 0, со = а среди узлов дискретизации х = 0, с*>. Диапазоны сосредоточения положительных времен и частот со, должны быть в обе стороны на порядок шире диапазона введенных положительных узлов х,, чтобы избежать проявления вырождения столбцов матрицы [а(1р X,)}, {¿(со,, X,)}. В результате типичная в приложениях размерность задач (П.11), (П. 12) весьма высока (.т ~ 500, п ~ 300).
При практической численной минимизации (П.11), (П. 12) по любому из цитированных мето-
дов неизбежно приходится решать плохо обусловленные системы линейных уравнений (типичен разброс сингулярных чисел матрицы системы, превышающий 10 порядков). Причиной, вызывающей плохую обусловленность, является само требование хорошей аппроксимации непрерывной задачи соответствующей дискретной (использование густой сети х для достижения высокого разрешения), а также малый штрафной коэффициент на финальных этапах минимизации методом внутренних штрафов. Этот фактор может быть критическим при неудачной программной реализации из-за замедления или даже расходимости итераций при потере значащих цифр. В нашей реализации алгоритмов на языке Fortran-90 точность численного решения промежуточных линейных уравнений достигается сочетанием (^-преобразований [13] и аппарата арифметики расширенной (четверной) мантиссы при вычислении частных производных и на этапе осуществления обратной подстановки в гауссовом исключении. Использована библиотека программ расширенной арифметики Дэвида Бейли (D.H. Beiley, A fortran-90 double-double library, Available at http://www.nersc.gov/dhbailey/mpdist/mpdist.html). В результате итерации сходятся быстро, в типичном случае при m = 500 затрачивается 2-5 секунд/вариант на персональном компьютере класса AMD Athlon 800.
СПИСОК ЛИТЕРАТУРЫ
1. Бартенев Г.М., Бартенева А.Г. Релаксационные свойства полимеров. М.: Химия, 1992.
2. Виноградов Г.В., Малкин А.Я. Реология полимеров. М.: Химия, 1977.
3. Baumgaertel M., Shausberger A., Winter Н.Н. // Rheol. Acta. 1990. V. 29. P. 400.
4. Winter Н.Н. // J. Non-Newtoian Fluid Mech. 1997. V. 68. P. 225.
5. Malkin A.Ya., Masalova I. // Rheol. Acta. 2001. V. 40. P. 261.
6. Малкин A.Я. // Высокомолек. соед. Б. 2002. T. 44. №9. С. 1598.
7. Дубовицкий В. А., Милютина И А. Метод гистограмм для численного анализа многокомпонентной кинетики реассоциации ДНК. Препринт ОИХФ АН СССР. Черноголовка, 1985.
8. Дубовицкий В.А., Ермолаев К.В. // Журн. физ. химии. 1996. Т. 70. № 7. С. 1233.
9. Ахиезер Н.И. Классическая проблема моментов. М.: Физматгиз, 1961.
10. Фелпс Р. Лекции о теоремах Шоке. М.: Мир, 1968.
11. Гохман Э.Х. Интеграл Стилтьеса и его приложения. М.: Физматгиз, 1958.
12. Шилов Г.Е., Гуревич В.А. Интеграл, мера и производная. М.: Наука, 1967.
13. Лоусон Ч„ Хенсон Р. Численное решение задач метода наименьших квадратов. М.: Наука, 1986.
14. Иржак Т.Ф., КузубЛ.И., Иржак В.И. // Высокомо-лек. соед. А. 2002. Т. 44. № 7. С. 1101.
15. Stepanek Р. // Dynamic Light Scattering, The Method and Some Applications / Ed. by Brown W. Oxford: Clar-endron Press, 1993. P. 177.
16. Provencher S.W. // Comput. Phys. Commun. 1982. V. 27. № 3. P. 213.
17. Rots T., Maler D„ Friedrich Ch., Marth M., Honerkamp J. Ц Rheol. Acta. 2000. V. 39. № 2. P. 163.
18. Rots T., Marth M., Weese J., Honerkamp J. // Comput. Phys. Commun. 2001. V. 139. № 3. P. 279.
19. Weese J. // Comput. Phys. Commun. 1992. V. 69. № 1. P. 99.
20. Baumgaertel M., Winter H.H. // Rheol. Acta. 1989. V. 28. P. 511.
21. Колмогоров A.H., Фомин C.B. Элементы теории функций и функционального анализа. М.: Наука, 1972.
22. Фиакко А., Мак-Кормик Г. Нелинейное программирование: методы последовательной безусловной минимизации. М.: Мир, 1972.
23. Goldfarb D., Idnani А. // Mathematical Programming. 1983. V. 27. P. 1.
24. Иржак В.И. // Докл. РАН. 2002. Т. 385. № 4. С. 470.
On the Problem of Consistent Construction of Relaxation Spectra from the Data of Mechanical Relaxation of Polymers
V. A. Dubovitskii and V. I. Irzhak
Institute of Problems of Chemical Physics, Russian Academy of Sciences, pr. Akademika Semenova 1, Chernogolovka, Moscow oblast, 142432 Russia
Abstract—A new formalism for description of relaxation curves is advanced; this approach is based on the notion of the Stieltjes integral with respect to the relaxation time distribution function of exponential modes. This approach allows a consistent description of all possible types of relaxation spectra, including discrete and continuous spectra, the initial instantaneous decay of a kinetic curve, and its stationary plateau. The key aspect of this approach concerns the existence of one-to-one continuous correspondence between the monotonic relaxation time function and the time or frequency dependence of the relaxation modulus of elasticity. Owing to this correspondence, the inverse problem of relaxation-time distribution reconstruction becomes mathematically well-posed in its classical meaning and does not require application of a regularization technique. An efficient algorithm of numerical solution is described. The examples of model kinetic data treatment for the modified Rouse chain and the experimental data on the frequency dependence of elastic modulus of linear polymers are presented. Therefore, a method for relaxation spectrum construction on the basis of the data on mechanical relaxation of polymers is proposed. The correctness of this approach is limited only by the adequacy of the applied mathematical model and by a systematic error and completeness of the experimental data.
