Search This Blog

Wednesday, November 13, 2013

Применение инновационных технологий при проектировании грузовых судов

Kuno van den Berg, менеджер проектов, Gijsbert Jacobse, инженер, Michiel Verdult, инженер, компания Vuyk Engineering Rotterdam b.v., Нидерланды

Специалисты голландской инжиниринговой компании Vuyk Engineering Rotterdam (VER) использовали компьютерное моделирование для разработки новых и модернизации существующих судов, их проверки на соответствие международным промышленным стандартам, а также устранения разнообразных проблем, возникающих при эксплуатации судна.

Специалисты Vuyk остановили свой выбор на ПО ANSYS после тщательного анализа конкурирующих программных продуктов — принимая во внимание гибкость программного кода, спектр используемых приложений, репутацию компании и использование программных продуктов ведущими мировыми компаниями.

С 2007 г. в компании VER используется ПО ANSYS AQWA для расчетов гидродинамических нагрузок на объекты морского строительства, включая анализ возникающих напряжений и усталостный расчет.

Внедрение ANSYS AQWA позволило инженерам проводить комплексный гидродинамический анализ и получить более полную картину, отражающую эксплуатационные характеристики судна. Кроме того, специалисты получили возможность проводить расчеты намного быстрее, изменяя параметры и сравнивая различные варианты проектов.
Увеличение грузоподъемности
В одном из проектов специалисты VER использовали программный комплекс ANSYS для увеличения грузоподъемности самоходного плавучего крана Matador 3, используемого в портах, на объектах морского строительства, при удалении обломков затонувших судов, а также при строительстве мостов на реках и каналах.

Специалисты использовали ПО ANSYS Mechanical для оптимизации грузоподъемности крана — подъёмная рама моделировалась с помощью балочных элементов, рассчитывались силы противодействия в точках изгиба и напряжения во всех деталях конструкции. Проведя параметризацию модели, команда инженеров смогла быстро ввести различные угловые переменные для получения кривых подъема и грузоподъемности.

Кроме того, специалисты VER выполнили ряд расчетов для проектов в открытых водоемах, включая анализ движения крана Matador при различной волновой нагрузке при подъеме груза, а также транспортировке и установке высоковольтной станции ветровых турбин в Северном море. Комплекс ANSYS AQWA использовался для определения максимальной высоты волн, допустимой для различных периодов волны (время между вершинами волны).
Распределение напряжения во всей подъемной конструкции Matador 3. Технологии - AQWA и ASAS. Журнал ANSYS Advantage. Русская редакция. 15 2011 
Распределение напряжения во всей подъемной конструкции Matador 3

Распределение напряжения в отдельных деталях - в опорных плитах. Технологии - AQWA и ASAS. Журнал ANSYS Advantage. Русская редакция. 15 2011
Распределение напряжения в отдельных деталях - в частности, в опорных плитах

Расчеты на ранних стадиях проектирования
Специалисты компании использовали программные комплексы ANSYS Mechanical и ANSYS AQWA в одностороннем связанном расчете, в котором гидродинамические нагрузки от давления на внешнюю сторону корпуса судна, полученные в ANSYS AQWA, прямо передавались в ANSYS Mechanical для определения прочности конструкции дноуглубительного судна. В частности, необходимо было проверить продольный изгиб в средней части судна, рассчитать напряжения при расчёте эквивалентного бруса, а также провести анализ усталостной долговечности на основе полученных напряжений.

Чтобы получить распределение напряжений по всему периметру корпуса при воздействии волн на судно, в ANSYS AQWA проводился трехмерный дифракционный анализ. За основу была взята конечноэлементная модель корпуса — обеспечивающая полную совместимость между конечноэлементным и дифракционным расчетами.

Инженеры VER сравнили нагрузки в стоячей воде с рабочими волновыми нагрузками, полученными в ANSYS AQWA. Затем эти данные использовались в ANSYS Mechanical для определения напряжений и деформации балок. В расчете выяснилось, что основная концентрация напряжений приходится на заднюю часть основной палубы. Добавление более плотных плит настила и балок позволило укрепить конструкцию. В последствии усталостный расчет показал, что измененная конструкция является более прочной и долговечной.
При дифракционном анализе определяется раскачивание и устойчивость поднимаемых грузов. Технологии - AQWA и ASAS. Журнал ANSYS Advantage. Русская редакция. 15 2011При дифракционном анализе определяется раскачивание и устойчивость поднимаемых грузов. Технологии - AQWA и ASAS. Журнал ANSYS Advantage. Русская редакция. 15 2011


При дифракционном анализе определяется раскачивание и устойчивость поднимаемых грузов

Оптимизация конструкции волновой энергоустановки

Bradford S. Lamb, президент, и Ken Rhinefrank, вице-президент, компания Columbia Power Technologies, LLC, Corvalis, США

Компания Columbia Power Technologies (COLUMBIA POWER) занимается разработкой волновых генераторов.


Преимущества использования энергии волн:
  • Удельная мощность: энергия волн намного плотнее по сравнению с другими источниками возобновляемой энергии. Подобные установки, размещенные на маленькой территории, вырабатывают большие объемы электроэнергии.
  • Прогнозируемость: количество энергии, вырабатываемой волнами, можно точно определить на несколько дней вперед.
  • Стабильность: в отличие от солнечной энергии, океан вырабатывает энергию 24 часа в сутки.
  • Близость к центрам потребления: энергия волн не требует существенных затрат на передачу электроэнергии, поскольку 37% населения Земли живет на расстоянии до 60 миль от побережья, а 70 % — на расстоянии 200 миль.
Однако волновая энергетика изучена гораздо меньше по сравнению с остальными способами генерации энергии. Компаниям-операторам необходимо значительно повысить эффективность и сократить затраты на проектирование, чтобы продемонстрировать потенциальным клиентам и инвесторам жизнеспособность и конкурентность продуктов. Компаниям, работающим в других отраслях энергетики для этого понадобились десятки лет — однако для волновой энергетики это слишком долго. Для достижения поставленных целей, необходимо быстро совершенствовать проекты при ограниченном бюджете на проектирование.

Основной проблемой стала оптимизация конструкции бакена — она должна быть оптимальной для выработки максимального количества энергии. Относительная ширина захвата — мера эффективности устройства при улавливании доступной энергии волны. Если эта величина равна 1, это означает, что бакен улавливает 100% доступной энергии волн.

Для оптимизации конструкции бакена, специалисты компании изучили возможности различных гидродинамических кодов и выбрали ANSYS AQWA из-за простоты использования и хорошего согласования результатов с экспериментальными данными. Инженеры COLUMBIA POWER оценили возможность получения в ANSYS AQWA решений в частотной области и нестационарной постановке.
изменение формы установки в зависимости от периода волны. Технологии - AQWA и ASAS. Журнал ANSYS Advantage. Русская редакция. 15 2011изменение формы установки в зависимости от периода волны. Технологии - AQWA и ASAS. Журнал ANSYS Advantage. Русская редакция. 15 2011
изменение формы установки в зависимости от периода волны. Технологии - AQWA и ASAS. Журнал ANSYS Advantage. Русская редакция. 15 2011изменение формы установки в зависимости от периода волны. Технологии - AQWA и ASAS. Журнал ANSYS Advantage. Русская редакция. 15 2011
Изменение формы установки в зависимости от периода волны

Этапы оптимизации конструкций в ANSYS

Оптимизации конструкций в ANSYS

Новая высокоэффективная технология многокритериальной оптимизации IOSO (Indirect Optimization on the base of Self­Organization; автор — д.т.н., профессор И. Егоров) от компании «Сигма Технология» базируется на собственных оригинальных методах и алгоритмах оптимизации, дающих возможность выполнить сложные задачи поиска оптимума, которые ранее не решались ввиду отсутствия эффективного метода. Стратегия решения задач оптимизации принципиально отличается от известных подходов нелинейного программирования, обладает более высокой эффективностью и обеспечивает существенно более широкие возможности. В соответствии с логикой работы алгоритмов IOSO на каждой итерации осуществляется построение поверхностей отклика критериев оптимизации и ограничиваемых параметров. Далее выполняется оптимизация с использованием полученных поверхностей отклика и в найденной точке проводится прямое обращение к математической модели исследуемой системы. В процессе оптимизации осуществляется накопление информации об исследуемой системе в окрестности оптимального решения, что приводит к повышению качества поверхностей отклика.
В данной статье мы рассмотрим процесс нахождения оптимальной конструкции металлической башни при помощи программного обеспечения ANSYS и алгоритмов оптимизации IOSO. Целью оптимизации является подбор наиболее экономически эффективных параметров конструкции антенной опоры башенного типа. Этот алгоритм также может быть применен для поиска оптимальной конструкции башенных кранов и различных вантовых сооружений.
С растущей в геометрической прогрессии зоной покрытия сигналом операторов сотовой связи требуется всё больше и больше точек для размещения антенного оборудования. Не всегда на такой территории есть дымовые трубы или подобные высокие объекты, подходящие для его размещения, поэтому часто приходится возводить новые конструкции, которые бы отвечали всем нормативным требованиям, но при этом их производство должно быть экономически обоснованно и эффективно, а это подчас два взаимоисключающих фактора.
Моделью является пространственная ферменная конструкция высотой 20 м, состоящая из десяти одинаковых секций. Секция имеет в плане квадратное сечение. Все элементы конструкции — равнополочные уголки. Направление раскосов и диагональных распорок по высоте чередуется для придания конструкции дополнительной жесткости (рис. 1). В качестве материала конструкции была принята сталь марки С245 (одна из самых часто применяемых при изготовлении проката).
 

Рис. 1. Конструкция башни
Для данного расчета был подготовлен входной файл, содержащий входные параметры, набор команд на языке APDL (Ansys Parametric Design Language) и позволяющий в автоматическом режиме произвести построение геометрической и конечно­элементной модели, вычисление и приложение внешних нагрузок, выполнить конечно­элементный анализ и обработку полученных результатов, а также записать выходной файл с результирующими параметрами, дающими возможность произвести проверку на соответствие условиям прочности и устойчивости несущих элементов. В качестве входных параметров использовались: ширина секции, количество ее делений (поясов) по высоте и номера профилей из сортамента в сечениях элементов.
 

Рис. 2. Ветровая нагрузка во фронтальном направлении на конструкцию башни
Нагрузками для башни являются: собственный вес, вес оборудования и ветровое давление на конструкцию и оборудование. Вес оборудования принят равным 100 кг, площадь наветренной поверхности антенн — 3 м2. Коэффициент надежности по нагрузке от собственного веса был принят 1,05, согласно СНиП 2.01.07­85* «Нагрузки и воздействия». Ветровая нагрузка также вычислялась по СНиП 2.01.07­85*. Для расчета был принят первый ветровой район с нормативным значением ветрового давления 23 кг/м2 (г.Москва), тип местности B — городские территории, лесные массивы и другие местности, равномерно покрытые препятствиями высотой более 10 м. Статическая и пульсационная составляющие ветровой нагрузки вычислялись в зависимости от изменения параметров конструкции и приводились к узловым силам, приложенным к каждой секции (рис. 2). Значение пульсационной составляющей вычислялось по формуле (9) в пункте 6.7 б) СНиП 2.01.07­85*, поэтому для каждого варианта предварительно выполнялся модальный анализ с учетом влияния массы установленного оборудования. Для нахождения максимальных усилий в элементах конструкции для каждого из вариантов производилось два расчета — при фронтальной и диагональной ветровой нагрузке.
Проверка несущих элементов осуществлялась согласно СНиП II­23­81* «Стальные конструкции». Раскосов и распорок — как центрально сжатых элементов, а стоек — как элементов, подверженных действию осевой силы с изгибом. Также проводилась проверка конструкции по перемещениям от ветровой нагрузки по п. 16.8, согласно которому они не должны превышать 1/100 от высоты всей конструкции, и проверка предельной гибкости элементов в соответствии с табл. 19.
 

Рис. 3. Различные варианты конструкции одной секции башни
При выполнении входного командного файла ANSYS производит два статических расчета в линейной постановке и модальный анализ (поиск собственных частот и форм колебаний), при этом для моделирования конструкции используется двухузловой балочный конечный элемент BEAM188. Всего в процессе подготовки модели к расчету по четырем независимым входным переменным автоматически вычисляется 84 зависимых параметра с использованием табличных данных из справочной и нормативной литературы, которые заложены во входной файл. По результатам расчетов перед завершением работы пакета ANSYS автоматически вычисляется семь выходных параметров: масса всей конструкции, три параметра — показатели прочности, отвечающие за проверку элементов по максимально возникающим в них напряжениям, еще два параметра — показатели гибкости, отвечающие за проверку элементов по критерию устойчивости и максимальное отклонение конструкции.
Программный пакет IOSO взаимодействует с решателем ANSYS напрямую, без использования интерфейса. Для этого программе IOSO нужно указать пути к входному файлу команд для пакета ANSYS, к выходному файлу и файлу запуска пакета ANSYS. Выходной файл содержит результирующие параметры в текстовом виде и автоматически создается по окончании работы расчетной модели. Файл запуска содержит командную строку запуска пакета ANSYS с соответствующими входными параметрами.
 

Рис. 4. Множество Парето
Постановка задачи оптимизации для IOSO состоит в следующем: независимые параметры: первый — ширина секции — варьируется в диапазоне 0,25...1,5 м с шагом 0,01 м; второй — количество поясов — меняется дискретно от 1 до 6; еще два независимых параметра — номера профилей из сортамента — варьируются от 1 до 116 для стоек и раскосов и распорок соответственно. Итого получается более 10 млн возможных вариантов исполнения данной конструкции в обозначенных диапазонах. Целевая функция представлена двумя критериями и пятью ограничениями: минимизация массы с ограничением до 4 т, минимизация отклонения мачты, все остальные выходные параметры ограничиваются значениями до 1, что соответствует удовлетворению условиям прочности и устойчивости. В качестве критериев остановки оптимизации можно использовать общее время ее выполнения, общее количество выполненных итераций (обращений к расчетной модели) либо оба критерия одновременно.
Получение результатов расчета для одного варианта конструкции занимает около 3­5 с. Если для расчета использовать многоядерный процессор, компьютер с несколькими процессорами или кластер, то можно применить параллельную версию программы — IOSO PM, которая позволяет одновременно запускать на расчет несколько расчетных вариантов. При этом для расчета каждого варианта конструкции также можно использовать несколько процессоров или ядер, указав соответствующий параметр в файле запуска пакета ANSYS. Это позволяет в разы увеличить скорость решения задачи.
 

Рис. 5. Оптимальный вариант конструкции башни
В процессе оптимизации происходит не просто перебор всех возможных комбинаций, а интеллектуальный отбор вариантов по алгоритмам IOSO, максимально приближающийся к целевой функции — снижению массы и отклонения с ограничением в виде проверок элементов на прочность и устойчивость. Таким образом, программа сама находит области поиска параметров внутри заданных пользователем диапазонов. Например, если все проверки удовлетворяются, то программа «понимает», что дальше увеличивать сечения элементов не нужно. В данном случае было выполнено около 500 итераций и получено множество Парето из 50 оптимальных вариантов исполнения конструкции (данное значение можно изменять), удовлетворяющих всем условиям, из которых можно выбрать наиболее подходящее (рис. 4). График множества Парето определяет ту самую золотую середину между надежностью конструкции и экономией денежных средств, но поскольку величина отклонения башни во всех случаях не превышает допустимого значения, равного 1, то выбрать можно любой вариант. В качестве оптимального варианта была выбрана башня с шириной основания 1,15 м, уголком 125Ѕ8 для стоек, 45Ѕ3 для раскосов и распорок и с секцией с одним делением по высоте. Масса конструкции составила ~1800 кг. Уровень напряжений в стойке не превышает 95% от максимально допустимых.
Таким образом, данная работа — только один из примеров применения программы оптимизации при проектировании изделий в строительной отрасли, но область ее приложения практически ничем не ограничена. Данный подход будет полезен как при разработке типовых изделий, что даст ощутимую экономию денежных средств, необходимых для массового производства изделия, так и при проектировании уникальных конструкций и сооружений. В настоящее время связка ANSYS и IOSO успешно применяется не только для промышленного и гражданского строительства, но и для решения производственных задач в сфере авиадвигателестроения и автомобилестроения и многих других. Оптимизация — это будущее проектирования и производства.
Задать интересующие вас вопросы вы можете на сайте Delcam­ural.ru в разделе «Вопрос­ответ» или на сайте Ansys.ru либо Ansysclub.ru  нашим специалистам. В следующей статье будет сделан обзор возможностей аэродинамического проектирования шахтного вентилятора в программной среде ANSYS Workbench.  


Wednesday, November 6, 2013

МОДЕЛИРОВАНИЕ КОНСТРУКЦИЙ В ANSYS










Линейный плоский треугольный элемент

В данной лекции будут рассмотрены особенности формирования системы линейных алгебраических уравнений метода конечных элементов (СЛАУ МКЭ) на примере трехузлового треугольного элемента с линейной интерполяцией перемещений, применяемого для решения плоской задачи теории упругости. Для краткости будем называть такой элемент линейным треугольным конечным элементом.
Этот элемент имеет ряд отличительных особенностей:
Он принадлежит к семейству так называемых изопараметрических элементов, о чем будет говориться в следующих лекциях;
Он позволяет получить выражения элементных матриц жесткости и элементных векторов сил в замкнутой форме, что означает отсутствие необходимости в численном интегрировании при вычислении элементных матриц жесткости и элементных векторов сил;
Точность решения, обеспечиваемая данным элементом, не может быть повышена путем добавления внутренних степеней свободы.
В дополнение хотелось бы отметить, что линейный треугольный конечный элемент имеет определенное историческое значение. Он был одним из двух первых конечных элементов, представленных в статье Мартина, Тернера, Клоха и Топпа в 1956 году. Эта публикация общепризнанно считается началом современного метода конечных элементов.
Хотя линейный треугольный конечный элемент в настоящее время реже используется при расчетах конструкции ввиду его низкой точности, тем не менее, он широко используется в тех случаях, когда нет необходимости в высокоточных расчетах, например, концентрации напряжений в конструкции. Другая причина широкого применения треугольного элемента состоит в том, что он очень удобен при использовании в алгоритмах автоматической генерации сетки, например, в широко известном и популярном алгоритметриангулизации по Делоне.


Tuesday, November 5, 2013

ЗАДАЧИ ТРЕХМЕРНОГО НАПРЯЖЕННОГО И ДЕФОРМИРОВАННОГО СОСТОЯНИЯ







ТРЕХМЕРНОЕ НАПРЯЖЕННОЕ И ДЕФОРМИРОВАННОЕ СОСТОЯНИЕ НА ОСНОВЕ МКЭ - Часть 1

1. Напряжённое состояние в точке


Вспомним, как нами было введено понятие напряжения. Рассмотрим тело, находящееся под действием системы уравновешенных сил (рис.48).
 Рис.1
 Будем исследовать внутренние силы в малой области окружающей точку , для чего проведём через данную точку сечение, рассекая тело на две части, и отбросим одну из них. Действие отброшенной части заменим внутренними силами (рис.49).



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



 Величина  1 зависит от размеров площадки. Чтобы избавиться от влияния размеров площадки deltaA , перейдём к пределу и будем стягивать площадку к точке B
  
                                                                                  (1)

Величину  будем называть полным напряжением в точке  по площадке с внешней нормалью .
Совершенно очевидно, что если мы выберем другую площадку, проходящую через точку B, но ориентированную другим образом, то в общем случае вектор полного напряжения окажется иным.
Совокупность всех векторов полного напряжения по всем площадкам, проходящим через данную точку, составляет напряженное состояние в данной точке.
Напряжённое состояние в данной точке известно, если известны напряжения по трём взаимно перпендикулярным площадкам, проходящим через данную точку.
Докажем это важнейшее положение.
Нас интересует напряжённое состояние в точке . Выделим в окрестности этой точки малый прямоугольный параллелепипед (рис.50). Размеры параллелепипеда настолько малы, что напряжённое состояние в пределах параллелепипеда можно считать однородным, и что грани параллелепипеда являются площадками, проходящими через точку B и имеющими внешними нормалями оси X,Y,Z




 Рис. 3

Полное напряжение можно разложить на три составляющие, направленные по координатным осям. Всего будем иметь 9 компонент напряжённого состояния: три нормальных и шесть касательных напряжений. Нормальные напряжения обозначим sizma и припишем индекс, указывающий внешнюю нормаль.
Например: sizma— нормальное напряжение по площадке с внешней нормалью X.
Касательные напряжения обозначаются Tau с двумя индексами. Первый индекс указывает площадку, второй – направление напряжения.
Например Tau— касательное напряжение по площадке с внешней нормалью X, параллельное оси Y.
Нормальное напряжение считается положительным, если оно направлено по внешней нормали, т.е. является растягивающим.
Касательные напряжения считаются положительными, если при положительной внешней нормали они направлены в сторону положительных координатных осей.
Совершенно очевидно, что по противоположным граням параллелепипеда действуют равные по величине и противоположные по направлению напряжения.
Заметим также, что хотя на рис.50 компоненты напряжённого состояния показаны в виде векторов, они являются величинами скалярными.


Докажем, что касательные напряжения по взаимно перпендикулярным площадкам равны. Т.к. тело, из которого вырезан элементарный параллелепипед, находится в равновесии, то условия равновесия применимы и к элементу объёма. Запишем условие, что сумма моментов всех сил приложенных к элементарному параллелепипеду относительно координатных осей равна нулю

Применение метода Аппроксимации в оптимизации конструкций

          В этом параграфе мы рассмотрим, как алгоритм вычисления аппроксимаций Паде реализован в таких популярных математических пакетах, как Maple, Mathematica и Matlab.
В Maple встроенная процедура pade(f, x = a, [n, m]) позволяет найти аппроксимацию Паде типа (n, m) в точке x = a функции f(x). Если при обращении к процедуре pade запросить доступную для пользователя информацию, указав, например, значение параметра infolevel[pade]:= 1, то будет получено сообщение об используемом для вычисления аппроксимации алгоритме. В случае, когда коэффициенты ряда Тейлора функции f(x) не содержат параметров и чисел с плавающей точкой, используется быстрый алгоритм Кабая и Чоя[2]. В остальных случаях нахождение аппроксимации проводится по алгоритму Геддеса [3],использующему символьные вычисления.
В случае неодномерного ядра матрицы (1) вектор, составленный из коэффициентов знаменателя, может быть найден не надлежащим образом, что может привести к появлению лишних нулей знаменателя и числителя. Покажем это на примерах.
Пример 1
Найдем аппроксимацию Паде π2,3(f1) типа (2, 3) функции f1(x) = x+1,0001 (x+1,999)(x−2,001) в точке x = 0. В этом случае ядро матрицы (1) неодномерно, однако, как видно из примера, найденные знаменатель и числитель не содержат лишних нулей.
restart;
with(numapprox);
padef1:= pade
x+1,0001
(x+1,999)(x−2,001), x = 0, [2, 3];
−0,2500250626−0,2500000625x1,000000000+0,0005000002006x−0,2500000625x2
factor(denom(padef1), complex); solve(numer(padef1) = 0);
0,2500000625(x + 1,999000000)(x − 2,001000000)−1,000100000.
 Пример 2
Найдем аппроксимацию Паде π4,5(f2) типа (4, 5) функции f2(x) = (x−3,001)(x+1,9999) (x2+1)(x+4,0001) вточке x = 0. В этом случае ядро матрицы (1) неодномерно, и найденные знаменатель и числитель аппроксимации Паде имеют лишние нули (они выделены жирным шрифтом).
padef2:= pade
(x−3,001)(x+1,9999)
(x2+1)(x+4,0001) , x = 0, [4, 5];
−1,500387465−0,3753104091x−0,4121913908x
2−0,08614086896x3+0,1068576988x41,000000000+0,3333333332x+1,448275862x2+0,4401910336x3+0,4482758628x4+0,1068577004x5
factor(denom(padef2), complex); factor(numer(padef2), complex);
0,1068577004(x + 4,000100004)(x + 0, 09748654036 + 1, 526433054I) (x + 0, 09748654036 − 1, 526433054I)(x + 1,000000001I)(x − 1,000000001I) 0,1068576988(x + 1,999900006)(x + 0, 09748653975 + 1, 526433060I) (x + 0, 09748653975 − 1, 526433060I)(x − 3,001000018).
Отметим, что в обоих приведенных примерах мы находимся в так называемом сингулярном блоке таблицы Паде для аппроксимируемой рациональной функции и, согласно теории аппроксимаций Паде, π2,3(f1), π4,5(f2) должны были оказаться равными рациональным дробям f1(x) и f2(x).
Найденная с помощью Maple аппроксимация π2,3(f1) оказалась равной аппроксимируемой функции, и в этом случае не появилось лишних нулей числителя и знаменателя. Однако аппроксимация π4,5(f2) не совпала с f2(x), у числителя и знаменателя аппроксимации появились близкие нули – дуплеты Фруассара. Дело в том, что из неодномерного ядра матрицы Tn+1 (1) в первом случае (случайно) был выбран правильный знаменатель. Во втором же примере взятый знаменатель не совпал со знаменателем рациональной дроби f2(x).
Похожая картина наблюдается и при попытке найти π2,3(f1), π4,5(f2) с помощью системы Mathematica. Как показывают приведенные ниже примеры, при этом опять появляются дуплеты Фруассара.
Пример 3
Вновь, как в примере 1, найдем аппроксимацию Паде π2,3(f1) типа (2, 3) функции f1(x) = x+1,0001(x+1,999)(x−2,001) в точке x = 0. Появляющиеся при этом дуплеты Фруассара выделены жирным шрифтом.
In[1]:= PadeApproximant[
x+1,0001
(x+1,999)(x−2,001), {x, 0, {2, 3}}]
Out[1]= −0,250025−0,222305x+0,0276927x 2
1,000000000000000−0,110271x−0,250055x2+0,0276927x3
In[2]:= Roots[Denominator[Out[1]] == 0, x]
Out[2]= x == −1, 999||x == 2, 001||x == 9, 02765
In[3]:= Roots[Numerator[Out[1]] == 0, x]
Out[3]= x == −1, 0001||x == 9, 02765
Пример 4
Как в примере 2, найдем аппроксимацию Паде π4,5(f2) типа (4, 5) функции f2(x) =(x−3,001)(x+1,9999)(x2+1)(x+4,0001) в точке x = 0.
In[4]:= PadeApproximant[
(x−3,001)(x+1,9999)((x2+1)(x+4,0001) , {x, 0, {4, 5}}]
Out[4]= −1,50039−25,4213x−6,51318x
2+3,7662x
3+0,42731x4
1,00000000000000+17,0263x+6,90326x2+17,4536x3+5,90326x4+0,42731x5
In[5]:= Roots[Denominator[Out[4]] == 0, x]
Out[5]:= x == −9, 75487||x == −4, 0001||x == −0, 0599743||x == −1, 0I||x == +1, 0I
In[6]:= Roots[Numerator[Out[4]] == 0, x]
Out[6]= x == −9, 75487||x == −1, 9999||x == −0, 0599743||x == 3, 001
Итак, числитель и знаменатель найденной по алгоритму, реализованному в Maple и Mathematica аппроксимации Паде, могут оказаться не взаимно простыми.