MATHEMATIQUES ET CAO

Volume 3

CALCUL SCIENTIFIQUE

par

Pierre Germain-Lacour Le chapitre 4 de cet ouvrage est du a P. L. George, F. Pistre

Collection dirigee par Yvon Gardan Professeur, Universite de Metz

Hermes Publishing

Volume 4

COURBES ET SURFACES

par

Pierre Bezier

Математика и САПР

В двух книгах

Книга —— 2

ВЫЧИСЛИТЕЛЬНЫЕ МЕТОДЫ

ГЕОМЕТРИЧЕСКИЕ МЕТОДЫ

Перевод с французского

канд. физ.-мат. наук В. В. Коваленко

и канд. физ.-мат. наук С. Д. Чигиря

под редакцией д-ра физ.-мат. наук Н. Г. Волкова

Москва «Мир» 1989

ББК 32.97

МЗЗ

УДК 681.3.082.5

Авторы: Жермен-Лакур П., Жорж П. Л., Пистр Ф., Безье П.

Математика и САПР: В 2-х кн. Кн. 2. Пер. с МЗЗ франц./Жермен-Лакур П., Жорж П. Л., Пистр Ф., Безье П. — М.: Мир, 1989. — 264с., ил. ISBN 5-03-000465-3

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

Для специалистов в области автоматизации процессов проектирования и кон- струирования на ЭВМ, а также для студентов вузов, изучающих САПР.

ББК 32.97

М 2М4060000-271111-89 041(01)-89

Редакция литературы по информатике и робототехнике

ISBN 5-03-000465-3 (русск.)

ISBN 5-03-000464-5 (русск.)

ISBN 2-88601-042-6 (франц.)

© Hermes Publishing, 1986, 1987

© перевод на русский язык, «Мир», 1989

ПРЕДИСЛОВИЕ РЕДАКТОРА ПЕРЕВОДА

Вторая книга монографии «Математика и САПР»0 состоит из двух ча- стей (части 3 и 4), которые посвящены общим методам вычислений в САПР и методам построения изображений.

В предисловии к первой книге я говорил о том, что программное обес- печение САПР призвано решить две задачи: 1) описание геометрических характеристик объектов, необходимое при их производстве; 2) описание функционирования объектов, необходимое при выборе их оптимальных ха- рактеристик. В соответствии с этим программное обеспечение САПР мож- но разбить на две части: программы геометрического моделирования и программы функционального моделирования. В первых большое место за- нимают алгоритмы и программы построения кривых и поверхностей, об- ладающих необходимыми свойствами гладкости. Ко вторым относятся программы численного решения уравнений, описывающих функционирова- ние систем.

Часть 3 посвящена численным методам, применяемым в САПР. Пред- ставлены численные методы решения систем обыкновенных уравнений, уравнений в частных производных. В отличие от большинства пособий на эту тему изложение носит сугубо прикладной характер.

Часть 4 посвящена разработанным автором методам аппроксимации кривых и поверхностей полиномами, зависящими от параметра. Эти мето- ды достаточно просты в реализации и дают неплохие результаты.

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

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

Книгу перевели В. В. Коваленко (часть 3) и С. Д. Чигирь (часть 4)..

Н. Г. Волков

Математика и САПР: в 2-х кн. Кн. 1. Пер. с франц. — М.: Мир, 1988.

Часть 3

ВЫЧИСЛИТЕЛЬНЫЕ МЕТОДЫ

П. Жермен-Лакур, П. Л. Жорж, Ф. Пистр

Глава 1

Введение

1.1. ОБЩИЕ ПОЛОЖЕНИЯ

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

Современные коммерческие САПР (CATIA, EUCLID, STRIM 100, SAPHIR, PAFEC, CADAM, CADDS4X, ANVIL, IDEAS) решают главным образом задачи геометрического представления изделий, а оценку их техни- ческих характеристик осуществляют с помощью вычислительных про- грамм, не входящих в состав указанных систем. Желательна более тесная связь между САПР и этими вычислениями, так как САПР предоставляет хорошие возможности пользователю: в частности, она позволяет описы- вать геометрические характеристики изделия и его частей, дает возмож- ность применять различные базы данных и средства графического пред- ставления результатов расчета. Эта связь, однако, затруднена рядом факторов:

1) если система геометрического описания объектов универсальна, то описание функционирования и физических свойств зависит от типа изделия и привлекает самые разные дисциплины (механику, термодинамику, элект- ричество и т. д.);

2) в ряде случаев время вычислений слишком велико и не дает возмож- ности использовать интерактивный режим работы;

3) если базы данных имеют узкую направленность, их трудно адаптиро- вать к новым задачам;

4) ресурсы графических средств ограниченны и часто с трудом справля- ются с такими задачами, как многократное воспроизведение линий тока, поверхностей уровня, переходных явлений и т. п.;

5) современные САПР велики по объему, что затрудняет их обслужива- ние и развитие. Так, например, система CADDS4X фирмы Computervision содержит около шести миллионов операторов исходной программы.

Введение i

Из сказанного ясно, что вычислительные средства, необходимые для описания физических характеристик объектов, должны развиваться как в рамках САПР, так и вне их.

В данной части книги излагаются методы вычислений, предназначенные для описания физических характеристик и функционирования объектов. Проблемы геометрического описания будут рассмотрены в следующей ча- сти. Метод изложения — от простого к сложному, при этом мы будем стараться избегать углубления в технические детали. Изложение будет ве- стись на инженерном, а не на математическом языке, хотя для каждой за- дачи будет приводиться ее математическое решение при необходимости с включением алгоритма или программы на языке Фортран. При этом мы будем отдавать предпочтение простым и надежным методам.

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

Итак, наша основная цель — ознакомить пользователей и специалистов с методами вычислений, используемыми в САПР.

1.2. ВЫЧИСЛЕНИЯ И САПР

Любое автоматизированное вычисление основано на построении моде- ли, которая описывает изучаемое явление и позволяет предсказать его раз- витие. С помощью алгебраических, дифференциальных и других уравнений модель представляет те свойства явления, которые необходимы для пред- сказания, и пренебрегает свойствами, не оказывающими влияния на иско- мые результаты. Модель содержит два типа физических величин: исходные данные и результаты. Для автоматизации вычислений каждую величину необходимо представить конечным множеством скалярных значений. Так, например, функция определяется набором ее значений для конечного мно- жества точек в области ее определения, а значения в других точках нахо- дятся интерполяцией.

Любое автоматизированное вычисление есть алгоритмическая обработ- ка конечного множества исходных данных с получением конечного множе- ства результатов. Вся эта информация (исходные данные и результаты) мо- жет располагаться — полностью или частично — в используемой САПР. Поэтому можно различать шесть уровней интеграции вычислений в зависи- мости от возможностей конкретной САПР. Они перечисляются ниже.

1.2.1. Внутренний алгоритм

В этом случае все данные содержатся в команде, запускающей вычисле- ния, и в базе данных САПР, а множество результатов помещается в ту же базу. Такой подход, дающий полную возможность использовать интер- фейсы с базой данных и функциональные средства самой системы (опцио- ны, графику, интерактивность и т. д.), типичен для разработчиков САПР. Фактически он сводится к увеличению системы на одну команду.

Глава 1

1.2.2. Макроинструкции

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

1.2.3. Внешний алгоритм с интерфейсом

В этом случае данные и результаты тоже размещаются в базе данных, но их обработку производит внешняя программа (в интерактивном или па- кетном режиме, это определяется спецификой вычислений), используя ин- терфейсы с базой данных, которые поставляются разработчиком САПР для этой цели. Примером может служить CATGEO для САПР CATIA. Та- кой подход часто практикуется специалистами по вычислительной технике предприятия, на котором применяется эта система.

1.2.4. Передача геометрических данных

Только часть данных и результатов размещается в базе данных САПР. С помощью САПР производится отбор данных, касающихся главным об- разом геометрии объекта, и формируется файл данных, который обраба- тывается обычным способом вместе с другими данными. Часть результа- тов может быть аналогичным образом возвращена в САПР. Этот метод удобен, если время обработки велико.

1.2.5. Перенос геометрических данных

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

1.2.6. Автономные вычисления

Обработка и управление данными и результатами производится вне САПР. Этим предоставляется большая свобода в организации вычисли- тельного процесса, но, с другой стороны, необходимость управлять всем — данными, результатами, визуальным отображением, документаци- ей и обучением пользователей — может оказаться слишком обремени- тельной.

Численные методы

Все сказанное объясняет интерес к расширению возможностей, которые дает стандартизация интерфейса, базы данных и языка пользователя. Сле- дующий этап — создание САПР, основанной на стандартной системе управления базой данных. К сожалению, многочисленные технические трудности тормозят это развитие, которое может быть столь же полезно, как и усилия, направленные на стандартизацию графических функций (GKS, PHIGS).

Глава 2

Численные методы

2.1. АЛГЕБРАИЧЕСКИЕ УРАВНЕНИЯ

Рассмотрим решение уравнений типа F(x) = 0. В большей части задач, приводящих к системе п линейных уравнений с п неизвестными вида Ах = Ь, матрица А является симметричной и ленточной. Ленточной назы- вается такая матрица, у которой все элементы вне полосы, расположенной вдоль главной диагонали, равны нулю. Полную информацию о ленточной матрице содержит таблица ее ненулевых элементов размерами п х /, где п — порядок матрицы, / — число отличных от нуля диагоналей.

При решении системы Ах = Ъ классическим методом Гаусса в эту таб- лицу вместо матрицы А помещается матрица U разложения А = LU. В данном разложении, которое является единственным, L—нижняя треу- гольная матрица с единичной диагональю, V—верхняя треугольная матри- ца. Реализующая метод Гаусса программа на Фортране (она приведена ни- же) вычисляет решение за время, пропорциональное п12, поэтому уравнения и неизвестные следует располагать так, чтобы лента была как можно бо- лее узкой.

SUBROUTINE GAUSSBD(A,В,N,LB,ND,EPS,KER) М М Е Е Е Е S * IMPLICIT REAL*8(A-H,0-2) DIMENSION ACND.LB) , B(ND) •

*Решение методом Гаусса системы А*Х = В с симметричной и ленточной матрицей А «вычисляется разложение А = LU, диагональные члены L равны 1

*А: на входе — строки А, начиная с диагонали, на выходе — строки U

* В: на входе — правая часть, на выходе — вектор решения »М: число неизвестных

* LB: ширина ленты, начиная с диагонального члена

*NB: размерность в вызывающей программе

*EPS: минимальное значение модуля ведущего элемента

*KER: код ошибки (0 при нормальном завершении)

10 Глава 2

К = -1 IF( N .LE. 0 .OR. N .GT. NO ) GOTO 9 IF( LB .LE. 0 .OR. LB .GT. N ) GOTO 9 KER = 0 NLB1 = N - LB + 1 JSUP = LB DO 4 К = 1 , N-1 PIV = A(K,1) IF( ABS(PIV) .LE. EPS ) GOTO 9 IMAX = К + LB - 1 IF( К .LE. NLB1 ) GOTO 1 IMAX = N JSUP = JSUP - 1 1 JMAX = JSUP II = 2 DO 3 I = K+1 , IMAX Т = A(K,II) / PIV JMAX = JMAX - 1 JJ = JSUP DO 2 J = JMAX , 1 , -1 A(I,J) = A(I,J) - Т * A(K,JJ) 2 JJ = JJ - I B(I) = B(I) - Т * B(K) 3 II = II + 1 4 CONTINUE К = N PIV = A(N,1) IF( ABS(PIV) .LE. EPS ) GOTO 9 B(N) = B(N) / PIV JMAX = 1 DO б I = N-1 , 1 , -1 S = B(I) JMAX = JMAX + 1 IP( I .LE. NLB1 ) JMAX = LB JJ = I + 1 DO 5 J = 2 , JMAX S = S - A(I,J) * B(JJ) 5 JJ = JJ + 1 6 В<1) = S / A(I,1) RETURN 9 KER = К RETURN END

Теоретически задача Ax = b разрешима, если матрица А регулярна, и в общем случае не имеет решения, если А сингулярна. С вычислительной точки зрения эта задача более или менее проста в зависимости от того, насколько А далека от сингулярности. Выявить эту зависимость можно, исследуя связь Дх (вариаций решения х) с АЬ (вариациями и) и ДЛ (вариа- циями А).

Проведем преобразования:

b + Д&, - ДЛ-Л-^, Ах= b, (А + ДЛ )(х + Дл-) ' (А + ДЛ)Дх= Дй Ах= (I + A~\ •ДЛ)-1 •Л-1 -(Дй - ДЛ -А'^). (2.1) (2.2) (2.3)

Численные методы 11

Далее мы воспользуемся неравенством

11(/+ М) (2.4)

1 - ИМИ '

которое выводится следующим образом:

(1+М)-(1+ М)~1 = I, (I+ М)-1 = I - М-(1 + М)~1, 11(7 + М)-Ч| < 11/11 + 11М-(/+М)-Ч| $ 1 + IIMII X IKZ+Mr'll.

Отсюда следует неравенство (2.4).

Используя выражения (2.3) и (2.4), можно написать

11Дл:11 s: IK^-A-1-^)--1!! х 11Л-Ч1 х (11Д611 + 11ДЛ11 х IIA-'&II), НА-' цду|| ^——"^—"——(ЦА611 + 11ДЛ11 х 1Ы1), 1 - ИЛ-^ДЛИ 11ДЛ11\ 11Л11 /' ИДд-11 ^ МП х \\А -'11 / \\Ш 11ДЛ1 Ы ^ 1 - \\А~\\\ х 11ДЛ11 \11А11 х 1Ы1 + 11Л11 но 11Й11 < 11Л11 х Ил-11, поскольку 11Ду|1 НА II Х НА"1 II /ИДйИ 11ДЛ11\ Ы ^ 1 - IIA-'ll х 11ДЛ11\ 11611 IL411/

Величина С(А) = \\А\\ х ИЛ"1 II называется числом обусловленности матрицы А. Если II ДЛИ мало, то IIAyll 1Ы1 . /11ДЙ11 ^ 11ДА11\ \ Ш \\А\\)'

поэтому можно сказать, что С(А) есть коэффициент усиления относитель- ных вариаций х по сравнению с относительными вариациями b и А. Оче- видно, что С (А) больше единицы и С(Л)-> + °о, если матрица А становит- ся сингулярной.

Из сказанного следует, что численное решение системы Ax = b с плохо обусловленными матрицами затруднено тем, что малые возмущения ис- ходных данных ведут к большим возмущениям решения. Примером может служить классический матричный тест — матрица Хессенберга, которая очень плохо обусловлена, хотя и регулярна, симметрична и положительно определена:

H, = 1 1 1 Т п 1 1 1 Т т п + '1 1 1 п 2л - 1

12 Глава 2 Численные методы 13

При решении систем с такими матрицами метод Гаусса неудобен, так как он не позволяет повысить устойчивость результата. Кроме того, число неизвестных системы Ax = b может быть столь велико, что матрицу А нельзя разместить в оперативной памяти (хотя матрица обычно разрежена, т. е. содержит много нулей).

В таких случаях представляет интерес применение итерационных мето- дов: Гаусса — Зейделя, верхней релаксации, сопряженных градиентов. При использовании этих методов в памяти хранятся только ненулевые элемен- ты матрицы А, достаточные для вычисления линейного преобразования w = Av. На каждой итерации выполняется преобразование и качество про- межуточного результата контролируется нормой остатка г = Ax - b.

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

— задать А и b (ненулевые элементы);

— задать х° и е;

— вычислить г° = Ах° - b;

— решить систему Мг° = г°;

— положить d° = z°;

— для k = О, 1, ... , п повторять:

v" = Ad", die = <^, />/, X*41 =Xic- Qkd", гk+l=rk-ekvk,

решить систему Мг* +1 = г*+',

0^i= (r^1, г^Лг*, ^>, ^+i =^^ +Р^^

до выполнения условия <г<:+l, г*"1"1) < е2 х (г°, 1я).

где D — диагональ А, Е — ее нижний треугольник, u — действительное число в пределах 1 < ш < 2. В последнем случае нет необходимости стро- ить матрицу М, а все вычисления можно проводить с элементами А. Опти- мальное значение ш находится анализом конкретного класса задач. Еще од- но достоинство метода сопряженных градиентов состоит в том, что этот метод позволяет эффективно использовать преимущества современных век- торных и параллельных компьютеров.

Описанный метод теоретически гарантирует сходимость за конечное число итераций. Но не все итерационные методы обладают этим свой- ством, поэтому во многих случаях приходится последовательно вычислять приближения л-о, х\, ..., Xk, ... к неизвестному решению по формулам типа Xk+i = f(Xk), где каждое приближение зависит от предыдущего.

Итерационный метод сходится, если

Ve, 3N такое, что п > N=> \Хп - х\ < е.

Метод имеет сходимость порядка се > 0, если

Ve, 3./V и с | 0 < с < 1 такие, что Xn+l - Х\ < Е. п > N= \Хп - Х

В частности, если се = 1, сходимость называется линейной, если а = 2, — квадратичной.

Метод Ньютона — Рафсона для решения нелинейного уравнения

F(x)=o: W ^+1=^-^

имеет квадратичную сходимость, но можно получить и формулу со сходи- мостью порядка 3, разложив F(x) в степенной ряд в окрестности решения:

Теоретически при отсутствии ошибок округления число итераций не больше порядка матрицы А. На практике же, если хорошо выбрана матрица М, итерационный процесс заканчивается гораздо раньше. Этот метод основан на замене системы Ах = b системой L ~ \AL ~ '(L'x) = L ~~ lb, где матрица L связана с матрицей М соотношением М = LL' и М близка к А. Матрицу М следует выбирать так, чтобы, с одной стороны, сущест- венно улучшить обусловленность L'~\AL~t по сравнению с обусловлен- ностью А, ас другой стороны, сделать время решения системы Отсюда следует M7k'\r\ = rk+l как можно меньшим. Среди предложенных в литературе ва- риантов выбора М следует отметить два особенно простых и эффектив- ных: диагональная матрица из диагональных элементов А и матрица верхней релаксации

M = 1-ш F(x) =0= F(u) + SF'(u) + ^^"(и),

где и близко кхи6=х-и. Если пренебречь членом с б2, получим форму- лу Ньютона. Пренебрегая членами более высокого порядка, получим

S-x-u-- F(u) - l ^("^WY F'(u) 2 F'(u) V'(")/ '

F(Xk) _ 1 F" (Xk) ( F(Xk) \2 +' * F' (Xk) 2 F' (Xk) \F' (Xk)) '

Xk+l =

В общем случае формула Xk+i = f(Xk) обеспечивает сходимость к х та- кому, что х =f(x), если - 1 f(x) = РХ + 1 - p^Xk+l = pXk + 1 - р,

14 Глава 2 Численные методы 15

где приближения сходятся к х = 1, если - 1 < р < 1. Если/' (х) = 0, имеем квадратичную сходимость, классическим примером которой служит формула

^•-К^)'

которая сходится к х = Ve.

Приведем еще две формулы для решения уравнения F(x) мощью приближений x^+i = f(Xk):

FF' f(x) •• х - (FT •• 0 с по-

Пусть А — симметричная и положительно определенная матрица, В симметричная матрица, yi — произвольный вектор. Последовательно вычисляются:

1 Ш = (Y'lAyi) 2, vi = ai у i, ^ уг = BAvi - (v\ABAvi)Vi, си = (у'гАуг) 2, VT, = агуг, Уз = BAv-t - (v'zABAv-i^Vi- (y\ABAv-i)v\, аз = (У'зАу,)2, из = тзуз, У4 = BAV} - (v'sABAv3)V3 - (v'zABAv^vi, 04 = (У'4Ау4) 2, V4 = ОС4У4,

сходимость которой имеет третий порядок, но становится линейной, если корень F(x) кратный, и

^"-(F^^FF"-

сходимость которой всегда квадратична независимо от кратности корня.

Эти формулы можно более или менее просто обобщить на случай, ког- да неизвестное х имеет несколько составляющих. Для примера напомним метод Бэрстоу0 для нахождения корней полинома

р(х) = аах" + dix"'1 + ... + dn-ix + an.

Метод состоит в том, что р(х) делится на трехчлен

t(x) = х1 - sx+ р^р(х) = 1(х)с[(х) + (w + п),

где го и /•I зависят от s и р, и методом Ньютона отыскивается вектор (^,), обращающий в нуль вектор {^).

Все описанные итерационные методы обеспечивают сходимость в неко- торой окрестности решения, которая может быть очень малой. По этой причине в определенных случаях более устойчивыми могут оказаться прос- тые методы, как, например, метод дихотомии для нахождения корня F(x):

1) найти а < b такие, что F(a) х F(b) < 0;

2) вычислить с = ——— и F(c);

3) если F(a) х F(c) > 0, положить д = с, в противном случа положить Ь=с;

4) вернуться к 2), если (b - а) > е.

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

1) Другое название — метод Хичкока. — Прим. ред.

yi+i = BAvi - (v',ABAv,)vi - (v,-iABAvi)v,-i,

Vn = OInVn.

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

VAV = I ^A~\ = W,

V'ABAV= Т = V~\BAV,

где Т — трехдиагональная матрица.

Решением системы Ax = b является

Решив эту систему, мы почти сразу получаем собственные векторы и собственные значения матрицы ВА, так как эквивалентная ей матрица Т — трехдиагональна. Если В = I, одновременно получаем обращение мат- рицы А, ее собственные значения и векторы. Если А = К и В = М~1, а М — диагональная матрица, получаем обращение К и собственные значе- ния и векторы М~ \K. Все это достигается за п последовательных итера- ций, состоящих главным образом в выполнении преобразований BAv,, ко- торые производятся довольно быстро, если матрицы А и В разрежены, хотя их нулевые элементы могут быть рассредоточены.

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

16 Глава 2

2.2. ЧИСЛЕННОЕ ИНТЕГРИРОВАНИЕ

Задача состоит в вычислении значения определенного интеграла

1= [ f(x)dx J n в одно-, дву- и трехмерных пространствах. Рассмотрим сначала случай одномерного пространства: х /= [ f(x)dx.

Для начала выберем внутренние точки xi, x-i, ..., Xn-i такие, что Хо < х\ < ... < Xn-i < Хл, и выпишем интеграл в виде точной суммы

f(x)dx,

с тем чтобы на каждом из отрезков [х, -1, xi] сделать функцию f(x) более регулярной. На всех отрезках [Xi-i, х,] проведем линейную замену перемен- ной интегрирования по формуле b - а Xi - Xi-\ и введем функцию Гх, - х,-\, . Л • ————(и - в) +x,-i , V b- a J у, — JC. 1 I \Xi — Xi - l „./,,•> _ •"' "'-' —1— f ————— (и - а) + х,-1 8'vл' ~ b- а ш(и) •/ 1 6 - а \ где а)(м), а и b — произвольны. Тогда * Во многих частных случаях простая формула b t-p \u(u)g(u)du = / ^A,g(ai), " ,=i в которой значения а,е[а, b] и А, определяются конкретным методом, дает точный результат для всех полиномов степени ниже 2р. Численные методы 17 На практике для численного интегрирования применяют следующие методы: • Гаусса — Лежандра: ш{и) = 1, а = - 1, Ъ = + 1; • Гаусса — Чебышева: ш(и) = ————— , а = - 1, b = + 1; VT^t? • Гаусса — Лагерра: ш(и) = е"", а = 0, b = + °о; • Гаусса — Эрмита: ш(и) = е""2, а = - °о, b = + оо. Среди формул, использующих крайние точки отрезка интегрирования, " р \,w(u)g(u)du = Aig(a) + ^A,g(ai) + Ap+ig(b), J a 1=2 проще других формулы метода Гаусса — Якоби4 для различных значений р (к ним относятся формула трапеций и формула Симпсона). Самым удобным и чаще всего применяемым является метод Гаусса — Лежандра. Приведем узлы и веса этого метода для различных значений р: р = 1, порядок 2, а\ = О, А\ = 2; р = 2, порядок 4, - ui = Д2 = -,=, А\ = A-i = 1; V3 Р = 3, порядок 6, - oi = «з = ;|а2=0, А, =Аз=|, Л2=| и т. д. Для вычисления определенного интеграла в дву- или трехмерном про- странстве следует разделить область интегрирования П на простые элемен- ты — треугольники и четырехугольники (в двумерном) и четырехгранники и шестигранники (в трехмерном). Для интегрирования по квадрату или ку- бу достаточно использовать формулы Гаусса — Лежандра: +1+1 „ „ Г Г Х^ V f(u, v)dudu= ^ /.А,А/(а„ dj), v v | f(u, v, w)dudv dw AjAkf(di, в/, д*,), где А, и а, — те же, что и для одномерных интегралов, а и, v и w локальные переменные, соответствующие рассматриваемому элементу. В методе Гаусса 13> - 1. — Прим. ред. Якоби а = - 1, Ь = 1, ш(и) = (1 - и}" (1 + и}8, а > -1, 18 Глава 2 Подобные же формулы для каждого значения р существуют для тре- угольников и тетраэдров. Ниже приведены две такие формулы третьего порядка. 1) Треугольник (рис. 2.1): f I = \f(u, v)dudv = Л[А/(а) + Bf(b) + Cf(c) + Df(d)] Рис. 2.1. Плоские барицентрические координаты •Ь \ •а / • • \ / /^ d \ - а = b = 1' 3 '3 5 ' 1 1 3 '1 1315" 1 53 А = В =с=D= 27 4825 5 '3 5 '1 5"' 4825 с - "5' 1 48 25 d =^- "5' 5 48 Д — площадь треугольника. 2) Тетраэдр (рис. 2.2): / = \f(u, v, w)dudvdw = V\Af(a) + Bf(b) + Cf(c) + Df(d) + Ef(e)] ^ Д. ^^ n a в \ --^ • • d* Рис. 2.2. Численные методы 19 Объемные барицентрические координаты 1 1 1 1 4 Q — ,4 —/-^ — - ~, > 4 "т у 4 4 ' '4 5 b= 1 1 1 1 д= 9 ~2' 6 ' 6' 6 20 с = 1 1 1 1 с= 9 ~~ 6 ' 2 ' 6' ^ 20 d = 1 =6' 1 '6' 12 ' 1 6 ?) = 920 1 1 1 1 9 6 = "6' 6' 6' 2" J7? = '20' F — объем тетраэдра. 2.3. ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ С НАЧАЛЬНЫМИ УСЛОВИЯМИ Пусть требуется решить задачу, представленную в канонической форме: У' = f(x, У), У = (х = хо) = уо, (2.5) где у — неизвестная скалярная или векторная функция х. Следует отметить, что если функция / не зависит от у, можно во- спользоваться методами, описанными в разд. 2.2, так как х У(х) =Уо+ \ f(t)dt. х. (2.6) В общем случае используется численный метод, состоящий в замене производной у'(х) = dy/dx конечными разностями. Эту замену можно про- вести различными способами. Простейший из них Уп-1 - Уп Хп + I — Хп при равенстве всех шагов Ал- = h приводит к явной формуле Эйлера: Уп-и = Уп + hf(xn, Уп). (2.7) Легко можно получить и другие формулы: • Неявная формула Эйлера Уп +1 = Уп + hf(Xn +1, Уп +1), (2.8) которая имеет тот же порядок погрешности, что и формула (2.7), но менее 20 Глава 2 удобна, так как для вычисления уп +1 требуется знать Уп +1, а эта рекурсия замедляет вычисления. •Формула трапеций h. Уп+1 = Уп + .[f(Xn, Уп) + f(Xn+l, Уп+l)}. • Формула (2.9) Уп+1 = Уп-1 + 2hf(Xn, Уп), (2.10) которая имеет два преимущества — явную форму и погрешность метода третьего порядка. Погрешность метода e(h) определяется разностью между точным зна- чением левой части и правой частью приведенных формул0. Из формулы Тэйлора следует /г2 А3 Уп+ i = Ул + hy'n + —Уп' + —Уп'+ - , 2 6 поэтому для явной формулы Эйлера e(h) = уп+1 - уп - hyn = ^уп = 0(h\ а для формулы (2.10) E(h) =уп+1- уп-i - 2уп ^ hyS= 0(h3). К сожалению, формула (2.10) не является A-устойчивой. Метод называ- ется А-устойчивым, если, применяя его к уравнению у' = — у, получаем Уп-^О при п-»оо независимо от величины h. На первый взгляд пригодны все методы порядка 2 и выше. Действи- тельно, погрешность О (А2) на один шаг приводит к погрешности О (h) на п шагов, и, поскольку х = nh, можно получить сколь угодно малую по- грешность для любого х, положив h достаточно малым. На самом деле этого недостаточно: нужна еще и устойчивость формулы, т. е. свойство, которое гарантирует ограниченность вариаций результата, вызываемых ограниченными вариациями данных. Все приведенные выше формулы устойчивы, но комбинируя две устойчивые формулы второго порядка Уп+1 = Уп-\ + 2hyn и уп = Уп-i + hy'n, получим формулу, шаг которой зависит не только от уп, но и от уп -1, Уп+i = -Уп + 2уп-1 + 3hf(Xn, Уп), (2.11) не обладающую свойством устойчивости. Применив эту формулу к реше- нию уравнения у' = 1 с уо = О и у\ = h + Е, получим 1) Это называется также погрешностью аппроксимации. — Прим. ред. Численные методы 21 y-i='ih- е, уз = 3h + Зе, у^= 4h - 5е, ys = 5/г + Не, Уп = nh е. Отсюда видно, что для данного х =.nh при /1->0итг-*оо^„не стремится к точному решению nh. Приведенные выше простые формулы почти не используются в числен- ном анализе из-за малого порядка погрешности. На практике применяют более сложные формулы, в том числе: • Формула Рунге — Кутта четвертого порядка ko = hf(xn, Уп), ki=hf(xn+^, у. ki = hf^xn + -, Уп+ — ks = hf(Xn + h, yn+ ki), 4), •)• Уп+i = Уп + ^(ko + 2ki + 2k2 + А-з), 6 (2.12) которая имеет погрешность РЫТОГО гюрядка и Л-устойчива, если h < 2,75. Применение этого метода Рунге — ^утта иллюстрируется приведенной ни- же программой на Фортране / «программа решения уравнения Y /-:д*у приА=-50 * методом Рунге - Кутта четвертого порядка *• метод А-устойчив, если Н <2."/5/ (—А) IMPLICIT REAL*8(A DATA A/-50.0DO/ P(X,Y)=A*Y O-Z) X=O.ODO Y=l.ODO ^тй'^20^--^2^" д0 ^'^'/(11X,"VALEUR DE H = ",F10.4)') H VsDEXP/"'1"'10 Е^АУ^*10 ^"'^(V-YVV „/ГЕ(9,'(4С20.6)') X,Y,V,E r/'»H*P(X,Y) /.,1=H"F(X+H/2 . ODO, Y+GO/2 . ODO) G2=H*P(X+H/2.0DO,Y+G1/2.0DO) G3=H"F (2.20) Уп-и = Уп + h[6yn'-n + (1 - в)уЛ которые при оптимальных значениях параметров се = - и в = - обеспечива- 6 2 ют погрешность пятого порядка. Ниже приведена программа, иллюстри- рующая реализацию этих формул в три итерации (число итераций выбрано произвольно). * Программа YQTDIF2: * Решение дифференциального уравнения *второго порядка: Y" = G (х, у, у ) » Рассматривается на примере A: Y° =i *- [OMEGA**2] *Y, YO, YPO «Используются формулы: *Y(N+1) = Y(N)+H*Y'(N)+0.5*H*H* [A*Y"(N+1)+ *(1-A)* Y" (N)1 »Y'(N+1) =Y'(N)+H*[T*Y"(N+ 1)+(1-T)*Y"(N)] * где Y"(N) обозначает G(X(N), Y(N), Y'(N)) * Известно, что метод оптимален при А=1/6 иТ=1/2 » IMPLICIT REAL*8(A-H,0-Z) G(X,Y,YP)=-OMEGA2"Y 24 Глава 2 OMEGA=2.0DO*3.14159265358979DO OMEGA2=OMEGA*OMEGA YO=1.0DO YPO^O.ODO A=1.0DO/6.0DO T=1.0DO/2.0DO H=0.01DO WRITE(9,'(IX)') WRITE(9,'(11X, "VALEUR DE OMEGA = ",F10.1)') OMEGA WRITE(9,'(11X, "VALEUR DE Y0= ",F10.4)'> YO HRITE(9,'(11X, "VALEUR DE YPO = ",F10.4)') YPO WRITE( 9, •( 11X, "VALEUR DE A= ",F10.4)') A WRITE(9,'(11X, "VALEUR DE Т = ",F10.4)') Т WRITE( 9,'(11X, "VALEUR DE H= ",F10.4)') H WRITE(9,'(IX)') X=O.ODO Y=YO "P=YPO YS=G(X,Y,YP) DO 10 1=1,1.0 V=DCOS(OMEGA*X) E=DABS(V-Y)/V WRITE(9,'(1G20.6)') X,Y,V,E YS1=YS Y1=Y+H*YP+0.5DO*H*H*(A*YS1+(1.ODO-A)*YS) YP1=YP+H*(T*YS1+(1.0DO-T)*YS) YS1=G(X+H,Y1,YP1) Y1=Y+H*YP+0.5DO*H*H*(A*YS1+(1.ODO-A)*YS) YP1=YP+H*(T*YS1+(1.0DO-T)*Y5) YS1=G(X+H,Y1,YP1) Y1=Y+H*YP+0.5DO*H*H*(A*YS1+(I.ODO-A)*YS) YP1=YP+H*(T*YS1+(1.0DO-T)*YS^ X=X+H Y=Y1 YP=YP1 YS=G(X,Y,YP) 10 CONTINUE WRITE(9,'(IX)') END Проводя сравнение различных методов решения дифференциальных уравнений с начальными условиями, следует руководствоваться следующи- ми критериями: • Каков порядок погрешности метода? • Требуются ли для реализации метода не предусмотренные в нем вы- числения? • Дает ли метод критерий оценки величины шага? • Позволяет ли метод легко адаптировать шаг интегрирования А? • Сколько раз надо вычислять функцию на каждом шаге (этим опреде- ляется длительность вычисления)? Численные методы 25 2.4. ПРОЕКЦИОННЫЕ МЕТОДЫ Рассмотрим задачу, определенную внутри и на границе области П в од- но-, двух- и трехмерном пространстве: (L(u)=f(x) Ул-еП, \S{u) = g(x) VA-es = эп. Задачи такого рода возникают в механике, электростатике, термодина- мике, акустике, теории диффузии и в других физических дисциплинах. В качестве меры близости вычисленного решения к точному будем ис- пользовать величину е(х) = L(u) - f(x), называемую невязкой. Удобно ис- кать решение и(х) в виде линейной комбинации независимых функций, каж- дая из которых удовлетворяет граничным условиям п и(х) = ^о^(д-). (2.22) *=i Все описанные в этом разделе методы пригодны для решения различ- ных дифференциальных уравнений в одно-, дву- и трехмерном простран- стве, но для простоты мы будем их рассматривать на примере простого частного случая (рис. 2.3): CL(u) = - и" + и = f(x) = 1 +(е - 1)х, [ и(0) = и(1) = 0, (2•23) имеющего точное решение и(х) = 1 + (е - 1)х - е". (2.24) Простейшие базовые функции для этого случая имеют вид ^(х)=х(\ -x)xk-\ но можно их выбрать и по-другому. Введем векторное пространство V, порождаемое линейными комбина- циями ipic, и скалярное произведение в этом пространстве 1 Г (v, w) = \v(x) x w(x)dx. V ' fw ° \uix) е •—————,, / О I •с ~~о'^ ~~~~~^*^•г' Рис. 2.3. 26 Глава 2 2.4.1. Метод коллокации: <е, 6(лс,)> = О Невязка Е(х), которая при точном решении должна быть равна нулю везде, приравнивается нулю в заданных точках, число которых равно числу слагаемых в формуле (2.22). Если это число достаточно велико и точки плотно расположены в 0, приближенное решение будет близко к точному. Ограничившись в выражении (2.22) одним членом, получим приближе- ние к решению задачи (2.23): и(х) = х(1 - x)oti. Отсюда е(х) = - и" + и - f(x) = (2ai - 1) + (ai - e + 1).»: - aiX2. Положив e ( . ) =0, получим 4 , . e - 1 -9°1 =\+-Г• Отсюда ai = 0,826285. Сравним полученное приближенное решение и с точным и^. х и и пр т 0,25 0,155 0,146 0,50 0,207 0,210 0,75 0,155 0,172 Если взять в выражении (2.22) два члена с коэффициентами ai и иг, то и(х) = х(\ - х)(оц + агх). Отсюда -и"+и=(1+х- х2)»! +(-2+6х+х2'- х^)аг и е(х,) = = 0 ^ (2 + х, - xf)ai +(-2+6х,+х}- xj)oit = 1 + (e - l)x„ Положив Xi = ( ., . ), получим 2,22222ai + 0,07407а2 = 1,57276 Cm = 0,698537, 2,22222щ + 2,14814a2 = 2,14552 => [02 = 0,276152. Сравним полученное приближенное решение и с и^. х и и пр т 0,25 0,144 0,146 0,50 0,209 0,210 0,75 0,170 0,172 Численные методы 27 Если выбрать х, равными, например, - и -, получим систему линейных 4 4 уравнений для ai и сег: 2,1875ai - 0,453125п2 = 1,42957 (оц = 0,711042, 2,1875ai - 2,640625 сч = 2,28871 => [tt2 = 0,277702. Сравним теперь полученное приближенное решение и с точным и^. х "пр "т 0,25 0,146 0,146 0,50 0,212 0,210 0,75 0,172 0,172 И наконец, если точек х, больше, чем коэффициентов а*, то и уравнений больше, чем неизвестных. В этом случае достаточно минимизировать сум- му квадратов разностей между правыми и левыми частями этих уравнений. Так, если сохранить в выражении (2.22) два члена и выбрать /I 1 3\ х, = I., ,, .1, получим 2,1875ai - 0,453125 П2 = 1,429570 m.o.m 2,25 .,+1,125 .2 = 1,859140 - [а1 = °' 42' 2,1875 а, + 2,640625 осг = 2.288711 la2 = 0'277702- Сравнение и и и^ дает "пр "т 0,25 0,14480 0,14555 0,50 0,21043 0,21042 0,75 0,17084 0,17171 a 2.4.2. Метод наименьших квадратов: -—<с, е> = О иск* Ограничиваясь двумя членами в выражении (2.22) E(x} = ai(2 + х - х1) + са(- 2 + 6л- + х2 - л-3) - (1 + (e - 1)х) и прирав- 28 Глава 2 i нивая нулю производные Q= e\x)dx no ai и а2, получим О 1 л -эе- = 0-» \е(х)(2 + х- xi)dx = О, доч J о i -дQ-=0^ Е(Х)(- 2 + 6х + х2' - xs)dx = О, да-г. о 564ai + 282а?. = 130е + 130^ (ai = 0,716224, 987ni + 1796o;2 = 441e + 14 => [»2 = 0,281653, x "пр "т 0,25 0,147 0,146 0,50 0,214 0,210 0,75 0,174 0,172 2-М Ш^А ^Ж^Ш: t'e, ^-1) =й Сохранив в выражении (2.22), как и ранее, члены с коэффициентами ai и «2, потребуем ортогональности невязки к -Xй и х1: i f \ e(x)dx = 0, J о i \ e(x)xdx = О, fJ о что дает 26tti + 13a2 = 6e + 6 С oci = 0,717222, 65ai + 63 02 = 20е + 10s' [сц = 0,281686. Приближенное и точное значение соответственно равны Численные методы 29 х "пр "т 0,25 0,148 0,146 0,50 0,215 0,210 0,75 0,174 0,172 2.4.4. Метод Гвлеркина: <е, ф,) = О Потребуем ортогональности невязки Е(х) и первых п функций ф,(х). Для 1анного примера: i Г е(х)х(1 - x)dx = О, J о i С 1 e(x)xl(l - x)dx = 0, V О 1то при п = 2 дает 22ai + 161 П2 = 5е + 35 Гщ = 0,769976, 77 ai + бОод = 21е + 14 => [02 = 0,196596. Отсюда Мд и и^ равны ^ "пр "т 0,25 0,154 0,146 0,50 0,217 0,210 0,75 0,172 0,172 Следует отметить, что метод Галеркина часто оказывается точнее дру- 'их, его сходимость лучше изучена и он служит основой метода конечных •лементов. ^.S. Метод Рэлея — Ритца: минимизация J(v) Часто можно найти такой функционал J(v), т. е. отображение введенно- о ранее пространства Ув пространство действительных чисел R, что иско- 30 Глава 2 мое решение минимизирует J(v). В нашем примере такой функционал существует0: 1 J(v) = [ [-^[(У)2 + ^] - [1 + (е - l)x]v}dx. о При п = 2 и м = х(1 - x)(cii + с<2-<) ^(У) есть симметричная положи- тельно определенная квадратичная форма от ei и аг: J(u) = 1^, а2)(Л)М - (bi b^"1] , 2 (tt2J (a2J минимум которой можно найти, решив уравнение w-0^ [:]-&]• В нашем примере значения cri и осг совпадают с решением, даваемым методом Галеркина, вследствие линейности оператора L(u) в формуле (2.23). Приведем другой пример применения метода Ритца — расчет деформа- ции балки на двух опорах с сосредоточенной нагрузкой в центре: F _Л_______У В Д N Д оо оо I -^————————————————»• Выберем для деформации v(x) выражение, удовлетворяющее граничным условиям . . . тгх . Зтгх v(x) = ai sin— + azsm——. Полная (внутренняя + потенциальная) энергия балки запишется в виде -if ($)—<*"• о Искомое решение v(x) дает минимум функционала T(v). Этот функцио- нал есть симметричная положительно определенная квадратичная форма 1) Уравнение (2.23) является уравнением Эйлера для функционала J(v) (см. разд. 2.8). — Прим. ред. Численные методы 31 от oii и а2, значения которых определяются решением двух линейных уравнений ?=0•^=0• dai do'2 Выполнив необходимые вычисления, получим 2F I3 2F I3 /1\4 "^^•^-^У- Точное решение в точке М представляет собой ряд с бесконечным чис- лом коэффициентов а* ^l^H^-)-^. Замечание. При использовании методов Галеркина и Ритца не всегда можно аналитически вычислить коэффициенты систем линейных уравне- ний, представленные в интегральной форме. В таких случаях следует ис- пользовать методы численного интегрирования (разд. 2.2). 2.5. МЕТОД КОНЕЧНЫХ РАЗНОСТЕЙ Этот метод состоит в замене производных соответствующим образом выбранными разностями. Изучив два приведенных здесь примера, чита- тель легко сможет решать и другие задачи. 2.5.1. Одномерная задача Рассмотрим задачу, обобщающую пример из разд. 2.4: (-и" +c(x)u=f(x), p25) [ и(0) = a, u(L) = р, где с(х) ^ -у >0. Разложим и(х + h) в ряд по степеням h: и(х + h) = и(х) + hu'(x) + ^и" (х) + ^и'" (х) + -^^(х) + ... . Отсюда следует ч"(х) = ^(х + h) - 2и(х) + и(х - h)] - -^t/4^) h 12 "Ри $6 [х - h, x+ h]. 32 Глава 2 Обозначим У| = u(i7i), где /i = ——,. Тогда n + 1 -'(-у,-^ 2м, - u,+ i) + сцч = fi. (2.26) ft Уравнения (2.26) образуют систему линейных уравнений относительно неизвестных и\, иг, ... , Un, при этом йо = cr и Un+i = 0. Эта система имеет ленточную симметричную трехдиагональную матрицу, что позволяет про- водить вычисления быстро и точно. В нашем примере легко показать сходимость приближенного решения к истинному при м-»оо, т. е. h-^0. Действительно, вычитая выражение (2.26) из (2.25), получим ,д(- 5,-i + 25, - 5i+i) + с,5, = - е,, h2 где 5, = и(х,) - Ui, е, = -.-г/4^,). Отсюда следует (— + с, | 5, = --(5,-1 + 5,+1) + е,, \h~ / /г (— + 7) 15,1 ^ 1 (— + Ci\S,\ < ^(Si.i\ + 15,+il) + le,l. Пусть Hell = Maxle,! и 11511 = Maxl5,l, тогда j j (-^ + т) 11611 ^ -^ (11811 + 11511) + Hell, 11511 ^llell =0(h2)•^0, Ч что и доказывает сходимость. Чтобы повысить порядок погрешности метода, разложим и(х + 2/г), и(х + h), и(х — h), и(х - 2/г) в ряд Тейлора в окрестности х: \[uk-i - 4uic-i + f>u„ -4ut+i + м*+2) = u^ + ^-uP + 0(/г4). я 6 Если учесть, что —(- «t-i + 2t——<>——II——1>——II г »>—II—•I—1>—1> Л= .^, 4 i " .А. \^. V Неизвестными являются значения и во внутренних узлах сетки, которые ^жни пронумеровать от 1 до л вдоль горизонтальных линий слева напра- в0- Для каждого <-го узла требуется знать номера соседних узлов и знать, 34 Глава 2 лежат ли они на границе или внутри области. В общем случае окружен: внутреннего ;-го узла выглядит следующим образом: k-1, ,k ,k+l i-1, ,i .1+1 J-1. •J •J+1 Для замены (Дм); можно использовать две аппроксимации - (Ди), = —(- itj - ui-i + 4у, - u,+i - Uk), h - (Дм), = —-(- iij-i - Uj+i + 4y, - Uk-i - Uk+i), 1h1 а также их линейную комбинацию с весами - и . : |~\| \~\Т| -1 -4 -1 -- - 1 4-1 +-—— 4 =—— -4 20 -4 2/г2 _i 32Й2 _i _i 6A2 _i _4 -i Погрешность аппроксимации последней формулы составляет -^^-й»——' Главная часть этого выражения представляет собой ДДм,, равное - Д/,. Поэтому можно оценить погрешность порядка О (и4), используя формулу —д-(- iij-i - Auj - Uj+i - 4u,-i + 20м, - 4м,+1 - w - Uk-i + Auk - Uk+i) = fi + — ( - fj - fi-1 + 4/, - /,+1 - A), которая ненамного сложнее реализуется, но обеспечивает меньшую по- грешность. Итак, метод конечных разностей прост в реализации, легко адаптирует- ся к различным типам обыкновенных дифференциальных уравнений и урав- нений в частных производных, включая нелинейные, и позволяет использо- вать недекартовы координаты (цилиндрические, сферические и любые кри- волинейные), но его использование ограничивается трудностью наложения регулярной квадратной сетки на область интегрирования. Эту трудность можно устранить, применив преобразование координат для приведения об- ласти к простой форме. Численные методы 35 Пусть дано преобразование: °\ Г' !/' Г^1 -_-~^/ п) ^r-J L^~—" ——————————————*> и 1———————————————»> х (да да \ х = а(и, v) (dx) _ Эу ^у \ Cdu) у = 0(и, v) ^ [dy\ = W Э^ [dv\ss ди ди / i дЦ _ да\ Cdu} 1 ду ду\ (dx} _ , Ы=^^[ ^ ^}[dyY rw D(u> v) = \ ди ди I _ да д^_да д^ ди ди ди ди Исходная краевая задача fZ(z)=/(J«r, у) в П, iz(x, у) = z(x, у) на Г преобразуется в (L"(w) = g(u, и) в П*, (^w(u, v) = ~w(u, v) на Г*. При этом очевидно, что w(u, v) = z(a(u,v), P(u, v)), g(u, v)=f(a(u, v), Р(и, и)), w(u, v) = z(a(u, v), t3(u, v)). Остается преобразовать производные г: , dz, дг. Qw 1 [др. да., 1 dz = —dx + —dy = — — —dx - —ду\ + дх ду ди D [ду ди ] дw I [ ЭР . да „ 1 + ^-— - —dx + —ду . ди D [^ ди ди -J 36 Глава 2 Отсюда следует 8w3f3 _ Qw8l3 az = G(u v) =JUJV JVJU = A("' v) дх ' ' ' дадР _ЭаЭ0 D(u, v)' ди ду ду ди Эй'дос 9w Эп ^ ^ „. . ^ 'ду~ди~~ди~ди ^ В (и, и) ду ' ' ' да_д^ _ ЭаЭ/З D(u, v) ' ди ду Эг ди Для преобразования вторых производных вычислим ю = _L i^-D - A aD^ дн=-L (^D - В31^ ди D1 \ди ди} ' ди ~ д2 \ди ди) ' ac_J_/an _ ао\ дн_ i /дв дв\ Qv~ С^^дУ дv}'дv~Dl\дv Qv) ' fz _ 1 /дОдр д0др\ дх2' D(u, У)\ди ду ~ди~ди} ' ^т. _ 1 /ЭО ди _ Ю да\ 1 /ЭЯЭр _ дН д13\ дхду D(u, и) \ди ди ди ду/ D(u, v) \ди ди ди ди/ Э2г = 1 (дНди _ дНди\ ду1 D(u, v) \ду ди ди ди/ Эти общие формулы в частных случах упрощаются. Так, при использс вании полярных координат имеем a^j^^a2^ i э2^ дх2 ду2 Эг2 г2 дв2 ' Заменой координат можно пользоваться при любом методе решенй но особенно она эффективна в сочетании с методом конечных разносте!' где необходима простота области интегрирования. 2.6. ДИСКРЕТНЫЕ СИСТЕМЫ Пусть требуется автоматизировать определение равновесного состоянг системы, состоящей из простых элементов (сопротивлений, импедансо! пружин, балок, труб и т. д.). Численные методы 37 рассмотрим сначала частный случай: с ——4ID 4 Р——— ——- г, ^ \ \\" U- ' /t. W/WA /• a В этой схеме имеются пять узлов и пять резисторов и наложены два условия по току (/2 = b и h = с) и два по потенциалу (fi = 0 и Vs = a). Для любого сопротивления R разность потенциалов и ток связаны урав- нениями i i " - --1 Г ^t Г ^^ CuJ-u,=Rh^ ^ R R ivic-VJ^Rik-j _1 1 у i R R] L J L ^fc •7 Отметим, что эта элементарная матрица симметрична и сингулярна. Определение равновесного состояния проводится за два шага. На пер- вом шаге записываются уравнения для токов (для каждого узла втекающий извне ток равен сумме внутренних токов, вытекающих из этого узла). Это приводит к системе линейных уравнений: Т7П-г-г^ 7 '•1 Г4 Г] Г4 -JL JL+-L --L „ /, ''i /•i /•2 гг .-_L -L+-L --L .3 = /з Г2 Г-, ГЗ ГЗ ^.27) --L -_L -L+-L+1 -^ ^ /4 Г4 Гз /•3 Г4 Г-, Г5 -1- ^ v, Л __ Г5 Г5 |__ 38 Глава 2 На втором шаге записываются условия для потенциалов (vi = 0 и vs = а), что дает J-+-1- | -± || U |71 [7 ri гг гг -\- -1- + -!- -\- У} = о + о . (2.28; /•2 Гг Гз Гз 1111 а -—— —— + —— + —— V4 С —— Гз Гз Г4 FS Г5 h= (}-+\\\ vi-^-vz-^-v^, (2.29) \Г1 /4/ Г1 7-4 1 ! Is = - —— V4 + —— Vs . rs rs Из регулярной системы (2.28) можно определить потенциалы u-i, УЗ и у», а из уравнений (2.29) токи /i и /5, втекающие в узлы, в которых наложе- ны условия по потенциалу. Из анализа этого примера можно сделать следующие выводы: 1) для двух соединенных параллельно резисторов складываются их про- водимости; 2) для последовательно соединенных резисторов складываются их со- противления. В нашем примере г-г. и гз примыкают к узлу 3, в который не втекает внешний ток и на который не подано напряжение. Найдем экви- валентное сопротивление для соединения гг и гз: - ' —1- 0 1 Г иг 1 Г h ~ Г Гг 1111 , -— — + — -— УЗ = 1з • гг гг Гз Гз о - -1- -!- J L уд J L h ГЗ Гз ^lt 2 Так как /з = 0, то УЗ = ——————. Подставляя Уз в уравнения для h и Гг + Гз /4 , получим j Численные методы • 39 Г i _ 1 i гг + гз гг + гз М ^ [li} 1 + 1 Ы Ы ' гг + Гз гг + Гз Итак, соединение г, г, I _' J *2——^——'^ можно заменить его эквивалентом Гг+г, *2——^ Такая операция называется сжатием внутреннего узла 3 по отношению к внешним узлам 2 и 4. Ниже мы рассмотрим другие применения этой операции. 3) Ясно, что нельзя для одного и того же узла задать внешние условия на напряжение и ток одновременно. 4) Если ни на один из узлов не подано напряжение, то сумма всех внеш- них токов должна равняться нулю. Физически это очевидно, математически следует из того, что матрица системы (2.27) сингулярна, как и любая эле- ментарная матрица. 5) Каждому узлу, на который подано напряжение, соответствует одно неизвестное и одно уравнение в системе (2.27), но ему же соответствует и одно дополнительное неизвестное—ток, который должен втекать в узел извне, чтобы поддерживать заданное напряжение, 6) Ширина ленты матрицы в системе (2.27) зависит от нумерации узлов, а нумерация сопротивлений не влияет на эту матрицу. Итак, определение равновесного состояния системы выполняется в пять этапов: 1) Пронумеровать сопротивления и узлы. 2) Для каждого сопротивления построить характеризующую его элемен- тарную матрицу и аддитивно распределить ее элементы в глобальной мат- рице, при этом места элементов определяются номерами узлов, которые соединяет данное сопротивление. 3) Записать граничные условия для узлов, на которые поданы на- пряжения. 4) Решить систему уравнений и получить значения напряжений в узлах. 5) Вычислить дополнительные неизвестные—токи, втекающие в узлы, на которые поданы напряжения. Эти действия можно распространить на системы, состоящие из импе- Дансов, пружин, балок и даже нелинейных величин, как, например, в ги- дравлике, где разность входного и выходного давлений гидравлического Узла пропорциональна квадрату расхода. Эту зависимость можно записать в форме Л - Рг = kd1 = (kd)d и искать положение равновесия итеративно: 40 Глава 2 приняв начальные значения для расходов, вычислить сопротивления kd и, решая систему, получить новые значения расходов и давлений, т. е. использовать итерационные приближения типа х = f(x). Можно также по- следовательно нагружать систему и переходить от одного состояния к дру- гому, используя малые приращения, как в методе Ньютона — Рафсона. Эти замечания сохраняют силу и при решении непрерывных задач мето- дом конечных элементов. 2.7. МЕТОД КОНЕЧНЫХ ЭЛЕМЕНТОВ ДЛЯ ОДНОЙ ПЕРЕМЕННОЙ Рассмотрим задачу (см. разд. 2.4 и 2.5): С-и" +u=f(x)= 1 +(е- 1)х, [м(0) = «(1) =0. Эта формулировка задачи, называемая локальной, не очень удобна для применения метода конечных элементов. Мы введем две другие формули- ровки, эквивалентные первой, но в большей мере соответствующие этому методу. Пусть v(x) — произвольная дифференцируемая функция, принимающая, как и и(х), нулевые значения на границах. Тогда 1 1 Г Г 1 (- и" + u)vdx = \ f(x)v(x)dx, V v 6 V, J «J о о и после интегрирования по частям первого слагаемого в левой части получим задачу (вторую формулировку) 1 1 Г Г \(и'и' + uv)dx= \fvdx, W^V, (P2) J J о о где V—пространство функций с непрерывными производными на отрезке [О, 1] и с нулевыми значениями на его границах. Введем функционал 1 Г Г J(v) = ^ \((v')2 + ^dx - \fvdx, J J о о который определяет третью формулировку задачи: Найти такую и(х), что J(v) ^ J(u), УгбИ (РЗ) 1 Численные методы 41 1 i Г Г Действительно, пусть а(и, v) = \(и' v' + и v)dx и , J .] " о о можно записать другую формулировку задачи J(u) ^ J(v) Vye V. (РЗ) Все три формулировки (PI), (P2) и (РЗ) эквивалентны. Следует отметить, что в формулировке (РЗ) J(u) интерпретируется как полная энергия— внутрен- "яя энергия плюс потенциальная энергия внешних сил, а формулировка (' ^—прямое выражение принципа возможных перемещений. Введем теперь функции, соответствующие разбиению [О, L] на п отрезков ^•(х) = О, \^,(x)=-(x-x,^)2^h-2(x-x,)], Щх) = (х - х,^)2^ + 2(х -х,)], h w,(x) =0, если х < (i - l)h, если (i-\)h ^ х < ih, если ih $ х < (i + l)h, если х ^ (i + l)h '38 Глава 4 Рис. 4.98. Поверхности Ci, Сг. ... , разделенные кривыми L и Л. Рис. 4.99. Обобщение вспомогательных поверхностей. К наиболее простому решению привело бы разбиение клеток кривыми L(u; 0,5) и Л(0,5; v). В этом случае полная поверхность состояла бы из множества клеток, описанных в разд. 4.5.2 (вспомогательные поверхности), и от исходных клеток остались бы только их центральные точки Р(и = 0,5; v = 0,5). Заключение В настоящее время системы автоматизированного проектирования уже не являются чисто техническим средством. Они оказывают большое эконо- мическое и социальное воздействие на развитие промышленности. Благода- ря им появилась возможность провести исследования в области математи- ки, механики, информатики и др. (рис. 1—10). Можно, по-видимому, уже сейчас сказать, что САПР являются реальной силой технического прогрес- са. Поэтому целесообразно сделать обзор и подвести некоторые итоги раз- вития САПР. 1. Скульптура, созданная с помощью ЭВМ. (Коллекция автора, 1968 г.) 2. Специальный обрабатывающий станок, используемый как построитель изобра- жений в трехмерном пространстве (SOFERMO}. 3. Макет легкового автомобиля (масштаб 1:8) и кожух сцепления (масштаб 1:1) в процессе обработки (SOFERMO). Система UNISURF 4. Детали кузова легкового автомобиля (SOFERMO). Глава 4 5. Макет транспортного средства будущего (RNUR). 6. Деталь паровой турбины (RNUR). Система UMSURF 9. Поверхность, называемая «Верне». Система UNISURF РОЛЬ САПР В ЭКОНОМИКЕ 247 10. Цветок. Еще недавно САПР считались не более чем модным увлечением, однако в настоящее время они применяются всюду и их экономическая эффектив- ность полностью доказана. Эта эффективность определяется не только до- стигнутой экономией, но и повышением точности изготовления изделий и ускорением освоения их производства. Повышение точности изготовления приводит к уменьшению общего времени освоения новой продукции, что позволяет более гибко реагировать на запросы рынка и конъюнктуру. Ко- личественно оценить эффективность от применения САПР довольно труд- но. Во всяком случае, это нужно делать не в книге, посвященной собствен- но САПР, а в каком-либо экономическом обзоре или непосредственно на производстве. В такой развитой отрасли, как автомобилестроение, вот уже 20 лет ши- роко используются станки с числовым программным управлением для про- изводства опытных образцов и серийных изделий. Они произвели настоя- щую революцию в производстве. Высказывавшиеся ранее опасения относи- тельно их низкой эффективности оказались беспочвенными. То же относится и к САПР при условии, если выбранные средства соответствуют поставленной цели. ОБУЧЕНИЕ ПОЛЬЗОВАНИЮ САПР Обработка в САПР в основном делается автоматически, но тем не ме- нее всегда существует возможность неудачного решения, полученного алго- ритмическим путем. В этом случае необходимо вмешательство оператора, что требует уделять больше внимания подготовке квалифицированного персонала для обслуживания САПР. Содержание этой подготовки известно уже давно, однако до сих пор в учебные программы большинства учебных заведений включены только краткие разделы, касающиеся САПР. Они содержат лишь геометрию в описательном изложении и основы анализа. Это считается большим пре- пятствием для широкого использования САПР в промышленности. В действительности эти опасения лишены оснований: для использования САПР необязательно знать теоретические основы ее работы. Практика по- казала, что сбоев при эксплуатации САПР, обусловленных низкой квали- фикацией персонала, почти не бывает. САПР вполне доступны таким ка- тегориям пользователей, как конструкторы, лаборанты, операторы ЭВМ. КАЧЕСТВО ПРОЕКТОВ, СОЗДАВАЕМЫХ С ПОМОЩЬЮ САПР В последние годы появилось большое количество публикаций, посвя- щенных обработке кривых и поверхностей. Часть из них имеет чисто мате- матическую направленность и непосредственно не связана с запросами 248 Глава 4 пользователей САПР. Это, казалось бы, лишает их права на существова- ние, однако опыт показывает, что многие из таких разработок, часто явля- ющихся маленькими открытиями, влекут за собой далеко идущие по- следствия. Что же касается чисто практических разработок, то их можно разбить на две категории. В разработках первой категории обрабатываемые объек- ты представляются традиционным образом с помощью отрезков прямых и дуг окружностей, комбинациями простейших геометрических тел. Не- смотря на простоту исходных положений, применение таких методов ста- новится весьма сложным для случаев определения зон сшивки между по- верхностями. В разработках второй категории решение основывается на па- раметрическом представлении кривых и поверхностей, что приводит к неполным решениям, обладающим невысокой точностью. Теоретические и практические исследования долгое время развивались независимо друг от друга, однако в настоящее время возникает потреб- ность в их сближении. Этот процесс уже начался, но еще потребуется вне- сти глубокие изменения в традиционные подходы к проблеме всех, кто уча- ствует в ее решении: математиков, конструкторов, специалистов по инфор- матике и другим областям науки и техники. Все идет к тому, что скоро станет необходимым получение с помощью САПР трехмерных объектов с помощью голографических средств или станочной обработки как для эстетической оценки объектов, так и для проведения их испытаний. Обра- батывающие станки станут необходимым дополнением к дисплейным экра- нам и графопостроителям, ДВА ПОДХОДА К ИСПОЛЬЗОВАНИЮ САПР Существуют два различных подхода к использованию САПР. При пер- вом сохраняются традиционные средства изготовления опытного образца. Если конструируемый объект имеет техническое назначение, сначала изго- тавливаются его чертежи, затем специалисты высокой квалификации изго- тавливают его с помощью шаблонов, выполняя плавные переходы между плоскими поверхностями. После этого проводят измерение координат большого числа точек, расположенных на поверхности объекта, и, наконец, специалисты по САПР переводят эту информацию в цифровую форму, до- ступную ЭВМ. Если же объект имеет только эстетическое назначение, чаще всего он предварительно изготовляется самим автором из воска или глины, затем модельщиками — из твердых материалов и только после этого он может быть представлен в виде цифровых данных. При втором подходе конструкторы и дизайнеры учатся использовать средства САПР в процессе разработки. Таким образом, замысел реализует- ся самими авторами, и информация остается в такой форме, в какой она возникла в процессе разработки. При этом исключается процесс перевода в цифровую форму, который иногда приводит к появлению ошибок. Литература 249 Необходимо отметить, что в настоящее время первый подход применя- ется более часто. Его особенность заключается в том, что конструктору не требуется менять привычный стереотип действий, сложившийся на про- тяжении долгого времени, и приобретать новые знания, которые, по его мнению, лежат в стороне от его основной деятельности. В то же время непосредственное использование САПР при работе позволяет ускорить процесс конструирования, избежать затрат на перевод данных в цифровую форму и связанных с этим ошибок. Эти преимущества указывают на безус- ловную предпочтительность второго подхода. В технике и технологии время от времени происходят перемены, и со- здание САПР является этому хорошим примером. Важно, однако, пони- мать, что этот процесс перемен происходит непрерывно и постоянно. ЛИТЕРАТУРА Часть 3 Глава 3 1. Conor J. J., Brebia C. A. Finite element techniques for fluid flow. Newnes- Butterworths, 1976.. 2. Wilkinson J. H. The algebraic eigenvalue problem, Clarendon Press, 1965. [Имеется перевод: Уилкинсон Дж. X. Алгебраическая проблема собственных значе- ний. — М.: Наука, 1970.) 3. Bathe К. J., El Wilson. Numerical methods in finite element analysis. Prentice Hall, 1976. 4. Zienkiewicz О. С. The finite element method in engineering science, MacGraw Hill, 1977. [Имеется перевод: Зенкевич О. Метод конечных элементов в технике. — М.: Мир, 1975.] 5. Gallagher R. H. Finite element analysis: fundamentals. Prentice Hall, 1975. [Имеет- ся перевод: Галлагер Р. Метод конечных элементов: Основы. — М.: Мир, 1984.] 6. Dhatt G., Touzot G. Une presentation de la methode des elements finis, Maloine S. A., 1981. 7. Imbert J. F. Analyse des structures par elements finis, Cepadues Editions, 1979. 8. Gardan Y. Mathematiques et CAD, Hermes, 1985. [Имеется перевод: Математика и САПР: В 2-х кн. 1 — М.: Мир, 1988.] 9. Yeysseyre H. Analyse Numerique, Cours ECP, 1986—1987. 10. Laurent-Gengoux P. Analyse de quelques problemes mathematiques en mecanique. Cours ECP, 1986—1987. 11. Mammeri В. Exercices d'analyse numerique, Cours ECP, 1986—1987. 12. Bezier P. Mathematiques et CAD — Courbes et surfaces, Hermes, 1986. [Имеется перевод: Математика и САПР. Перевод книги содержится в данной книге в ч. 4.] 17-1172 250 Литература Глава 4 1. MODULEF: "Une bibliotheque modulaire d'elements finis", INRIA, 1985. 2. MODULEF: "Generation automatique de maillages", P. L. George, INRIA, 1986. 3. Diaporama: Bibliotheque Modulef: "Quelques realisations concretes de maillages", P. L. George, E. Saltel, INRIA. 4. "SIMAIL: "Logiciel de maillages" documentation technique et descriptif des commandes version 3.1., juillet 86, SIMULOG. Часть 4 1. American Machinist, "How Cross uses special computer in designing special ma- chines", 12 nov. 1966. 2. Andrighetto, Verge. Determination des contours apparents et des zones en contre- depouille. Memoire de 1'Ecole National Superieure des Arts et Metiers, 1974. 3. Ansaldi S., de Floriani L., Falcidieno B. Geometric Modelling of Solid Objects by Using a Face-adjacency Graph Representation, Computer Graphics, July, 1985. 4. Araki H., Amano Y. Replacement of physical models by highly accurate numerical models for manufacturing stamping dies in the integrated CAD/CAM system, Pro- lamat 1985. 5. Atkins D. A., Tapia R. A. Mathematical ship lofting, Society of Naval Architects and Marine Engineers, Santa Barbara Congress, 13 oct. 1962. 6. Auffray, Pichouron. Definition des surfaces developpables s'appuyant sur une courbe. Application a la determination des surfaces annexes des outils d'emboutissage, Memoire du Conservatoire National des Arts et Metiers, juin 1974. 7. Automotive Engineer, "Advanced techniques for body engineering (GMC)", march 1969. 8. Ball A. A., Consurf, Introduction of the conic lofting tile. Computer Aided Design, 6—4, oct. 1974. 9. Ball A. A., Storry D. J. T. Recursively generated B-spline surfaces, CAD, 84. 10. Barnhill R. E., Smooth interpolation over triangles. Computer Aided Geometric Design, Academic Press, 1974. 11. Barnhill R. E. Smooth interpolation over triangles, CAGD 1974. 12. Barnhill R. E. Surfaces in computer-aided geometric design a survey with new results, CAGD, sept. 1985. 13. Beeker E. Smoothing of shapes designed with free-form surfaces, CAD, may 1986. 14. Bernard F. Surfaces complexes lissees selon un critere d'energie minimale, Avions Marcel Dassault. 15. Bezier P. Numerical control, Mathematics and application, J. Wiley, 1972. 16. Bezier P., Sioussiou S. Semi-automatic system for defining free-form curves and surfaces, CAD, march 1983. 17. Bezier P. UNISURF, from styling to tool shop. Proceedings of C. A. P. E. 83 meeting. North Holland 1983. 18. Bezier P. Emploi de la commande numerique pour 1'etude des carrosseries et des outillages, J. de la Societe des I'Automobile (S. I. A.), Avril 1968. 19. Bezier P. Essai de definition numerique des courbes et des surfaces experimentales, These, Universite de Paris VI, 1977. 20. Board of Trade Journal, Electronic control of machine tools for contour milling (supplement, 1963). Литература 251 21. Boehm W. Generating the Bezier points of triangular splines, SGAGD, 1983. 22. Boehm W. Generating the Bezier-points of triangular splines. Surfaces in CAGD, North Holland, 1983. 23. Boehm W. The de Boor algorithm for triangular splines, SCAGD, 1983. 24. Boehm W. On the efficiency of knot insertion algorithm, CAGD, mai 1985. 25. Boehm W. Triangular spline algorithms, CAGD, mai 1985. 26. Boehm W. Multivariate spline method in CAGD, CAD, mars 1986. 27. Boehm W. Curvature continuous curves and surfaces, CAD, mars 1986. 28. Braid I. C., Hillyard R. C., Stroud I. A. Stepwise construction of polyhedra in geometric modelling, MMCGD. 29. Braid I. C. From geometric to product modelling, PROLAMAT 1985. 30. Brodlie К. W. A review of methods for curves and function drawing MMCGD, Academic Press, 1980. 31. Brunei P. Increasing the smothness of bi-cubic spline surfaces, CAGD, mai 1986. 32. CAD, Computer Aided Design, Butterworth. 33. CAD. 84, Actes du congres CAD de Brighton, 1984, Butteworth. 34. CAGD, Computer Aided Geometric Design, Elsevier-North-Holland. 35. CAGD. 74, Actes du colloque CAGD-Salt Lake City 1974, Academic Press. 36. Caja E., Garetti R. NC grinder makes press-tools, Metalworking Production, may 1965. 37. Casteljau P. (de), Mathematiques et CAD, Volume 2: formes a poles. Edition Hermes 1985. [Имеется перевод: Математика и САПР: В 2-х кн. Кн. 1 — М.: Мир, 1988.1 38. Catmull E. A class of local interpolating splines. CAGD 1974. 39. Chalvet M. Machine a dessiner et a sculpter, precede Unisurf Renault/CII, Ingenieurs et Techniciens, juin 1969. 40. Chang G. Z. Matrix formulation of Bezier techniques, CAD 14, № 6, nov. 1982. 41. Chang G. Z. An elementary proof of the convergence for the generalized Bernstein- Bezier polynomials, University of Science and technology of China. Hefei, Anhui (People's Republic of China). 42. Chang G. Z. Families of adjoint patches for a Bezier triangular patch, CAGD, mai 1985. 43. Chang G. Z. The application of Bezier's functions to a proof of the Bernstein's theorem, University of Utah Salt Lake City. feb. 1981. 44. Chang S. J. A class of matrix methods for surface representation. Academy of China. 45. Chemla B. Systeme conversationnel de determination mathematique de courbes et surfaces. These de doctoral, Universite de Paris VII, 1971. 46. Chen J., Wu J. H. Non-modal interpolation, an improved technique for the uni- forme B-spline method, CAD, 1984. 47. Chung W. L. A new method of view synthesis for solid modelling, CAD, 1984. 48. Coons S. A. Surfaces for computer-aided design of space forms, Technical report MAC-TR-41, M. I. T. Cambridge, Mass, juin 1967. 49. Coons S. A. Surface patches and B-splines curves, CAGD, 1974. 50. Dahmen W., Micchelli C. A. Line average algorithm: a method for the computer generation of smooth surfaces, CAGD, mai 1985. 51. Daulte J. J. Informations Schnittstellen zwischen Beischreibung und Fertigung gek- rummter Flachen, Ecole Federale Polytechnique de Zurich, These de doctoral, 1974. 52. Dokkan T. A design system for products with sculptured surfaces, CAD, 1984. 53. F-Mesh, Computer program manual, Fujistu Ltd. 54. Farin'G. Bezier polynomials over triangles and the construction of C' polynomials, Brunei University. Uxbridge, Great Britain (1980). 252 Литература 55. Farin G. Algorithm for rational Bezier-curves, CAD, march 1983. 56. Farin G. Smooth interpolation to scattered 3D data, SCAGD, 1983. 57. Farin G. A modified Cloud-Tocher interpolant, CAGD, sept. 1985. 58. Farin G. Piecewise triangular C' surface strips, CAD, fev. 1986. 59. Fayard G. Les fabrications mecaniques a partir des equations mathematiques, La Technique Moderne, avril 1961. 60. Ferguson J. Multivariable Curve interpolation, J. of the Association for Computing Machinery, avril 1964. 61. Filip D. J. Adaptive subdivision algorithm for a set of Bezier triangles, CAD, mars 1986. 62. Forrest A. R. Computational Geometry, achievements and problems, CAGD, 1974. 63. Forrest A. R. Recent work on geometric algorithms, MMCGD, 1980. 64. Forrest R. Notes on Chaikin's algorithm, Computational geometry Memorandum (CGM 74-1) University of East Anglia, G. B. 65. Forrest R. A. Interactive interpolation and approximation by Bezier's polynomials, Computer Journal, 151, 1972. 66. Furukawa S., Kobatake S., Ishizuka G. Succinct representation of 3-D objects and its application to hidden-line problem, PROLAMAT, 1983. 67. Gardan Y., Mathematiques et CAO — Volume 1. Methodes de Base. Editions Hermes. [Имеется перевод: Математика и САПР: В 2-х кн. Кн. 1. — М.: Мир, 1988.] 68. Germain Lacour P. et. al. Mathematiques et CAO, Volume 3: calcul scientifique. Editions Hermes. [Русский перевод содержится в данной книге в ч. 3.) 69. Goldman R. N. Using degenerate Bezier triangles and tetrahedra to subdivide Bezier curves, CAD, nov. 1982. 70. Gordon W., Riesenfeld R. Bernstein-Bezier methods for the computer-aided design of free-form curves and surfaces. J. of the Association for Computing Machinery, april 1974. 71. Gordon W. L. Blending-function Methods of Bi-variate and Multivariate Interpo- lation and Approximation. S. I. A. M. — J. numer. Anal., № 8, 1971. 72. Cordon W. L., Riesenfeld R. B-splines and surfaces, CAGD, 1974. 73. Gregory J. A., Smooth interpolation without twist constraints, CAGD, 1974. 74. Gregory J. A. Shape preserving spline interpolation, CAD, fev. 1986. 75. Guenn, Lagain, Representation automatique d'une forme tridimensionnelle, Effa- cement des parties cachees, Memoire de 1'Ecole Nationale Superieure des Arts et Metiers, 1974. 76. Haapaniemi A., Nagase H., Fujimoto M. et al. Development of a real time nume- rical controller for machining of sculptured surfaces, PROLAMAT, 1985. 77. Haschek J. Dual Bezier curves and surfaces, SCAGD, 1983. 78. Henderson M. R. Extraction and organization of form features, PROLAMAT, 1985. 79. Higashi M., Kohzeni J., Nagasaka J. An interactive CAD system for construction of shapes with high quality surfaces, PROLAMAT, 1983. 80. Ho-bin, Inputing constructive solid geometry representations directly from 2D orthographic engineering trawings, CAD, april 1986. 81. Hoschek J. Dual Bezier-curves and surfaces. Surfaces in CAGD, North Holland, 1983. 82. Hoschek J., Smothing of curves and surfaces, CAGD, mai 1985. 83. Jared G. E. M., Varady T. Synthesis of volume modelling and sculptured modelling, CAD, 1984. 84. Johnson W. L. Interactive computer graphics as an engineering tool, Document SAE 69 0192, Detroit congress, Jan. 1969. Литература 253 85. Kahmann J., Continuity of curvature between adjacent Bezier-patches, Surfaces in CAGD, North Holland, 1983. 86. Kahmann J. Continuity of curvature between adjacent Bezier patches, SCAGD, 1983. 87. Lasser D., Bernstein-Bezier representation of volumes, CAGD, mai 1985. 88. Lasser D. Intersection of parametric surfaces on the Bernstein-Bezier representation, CAD, mai 1986. 89. Le Tac G. Definition des courbes et surfaces an langage APT et ADAPT, Mecanique et Electricite, aout 1968. 90. Lon R. B-spline convex surfaces, Beijing University, 1980. People's Republic of China. 91. Machinery and Production Engineering, 15 sept. 1965, Large Cramic tape controlled drafting machine. 92. Machinery and Production Engineering, 4 June 1969, "ICL off-line system for computer produced drawings". 93. Machinery, April 1970, "New computer numerical control system". 94. Machinery, May 1969, "Chevrolet is using computers, NC, EDM, abrasive ma- chining and advanced casting techniques to reduce stamping die cost". 95. Marchesseau C. La programmation point par point — Bulletin du GAMI, № 3, 1967. 96. McDermott R. J. Graphical representation over triangles and rectangles. Computer Aided Geometric Design, Academic Press 1974. 97. Mehlum E. Curve and surface fitting based on variational criteriae for smoothness, Central Institute for Industrial Research, dec. 1969. 98. Mehlum E. Non-linear splines, CAGD, 1974. 99. Merris D. K. New approaches to computer-aided design, Product Engineering, 3 Feb. 1964. 100. Middletich A., Sears K. H. Blended surfaces for set theoretic modelling systems, Computer graphics, July 1985. 101. Minagawa M., Okino N.. Kakazu Y. Development of full automatic dimensioning system based on 3D solid geometry, PROLAMAT 85. 102. MMCGD: Mathematical Methods in Computer Graphics Drawing, Academic Press, I960. 103. Moreau R. Machine mathematique a usiner les formes compliquees, Science et Techniques, mai 1966. 104. Nielsen G. M. Some piecewise polynomial alternatives: splines under tension, CAGD, 1974. 105. Nielsen G. M., Franke R. Surface construction based upon triangulation, SGAGD, 1983. 106. Nilson G. Computation of v-splines, Arizona State University, juin 1974. 107. Parizot S., Bouffard-Vercelli R. Application de la commande numerique a 1'etude et a la realisation d'une carrosserie. Encyclopedic des Metiers, la Carrosserie, Edition des Compagnons du Tour de France. 108. Parkinson A. The use of solid models in BUILD as a database for NC Machining, PROLAMAT 1985. 109. Peters G. J. Interactive computer graphics application of the parametric bi-cubic surface to engineering design problems, CAGD, 1974. 110. Petersen D. P. Boundary to constructive solid geometry mapping — a focus on 2D issues, CAD, fevr. 1986. 111. Piegl L. Representation of quadric primitves by rational polynomials, CAGD, mai 1985. 254 Литература 112. Pilcher D. Smooth parametric surfaces, CAGD, 1974. 113. Ponte J. Le systeme "Profile Data", Mecanique et Electricite, aout 1968. 114. Prautzsch H. Generalised subdivision and convergence, CAGD, mai 1985. 115. PROLAMAT 69. Actes du Congres Prolamat de 1969 a Rome, North Holland. 116. PROLAMAT 83. Actes du Congres Prolamat de 1983 a Amsterdam, North Holland. 117. PROLAMAT 85. Actes du Congres Prolamat de 1985 a Paris, North Holland. 118. Prudhomme M. R. Applications industrielles de la CN. Societe des Ingenieurs de 1'Automobile (S. I. A.), mai 1963, mars 1964. 119. Richards Т. H., Onwubolu G. C. Automatic interpretation of engineering drawings for 3D surfaces representation in CAD, CAD, april 1986. 120. Riesenfeld R. Application of B-spline approximation to geometric problems of computer-aided design, These, mars 1974, Computer Science University of Utah, Salt Lake City, USA. 121. Ris G. Raccordement a 1'ordre N entre carreaux de surfaces definies par des polynomes parametriques a coefficients vectoriels. Application aux surfaces UNISURF, These de doctorat de 1'Universite de Nancy, April 1975. 122. Sabin M. Parametric splines in tension, British Aircraft Corporation, juillet 1970. 123. Sabin M., Spline surface, British Aircraft Corporation, Juin 1969. 124. Sabin M. Trinomial basis functions for interpolation in triangular regions (Bezier- triangles). Notes British Aircraft Corporation, July 1971. 125. Sabin M. A. Some negative results in л-sided patches, CAD, fev. 1986. 126. Sablonniere P. Interpolation d'Hermite pour des surfaces de classe C' quadratiques par morceaux, Universite d'Ascq-les-Lille, 1981. 127. Sablonniere P. Bases de Bernstein et approximants splines. These de doctorat, Uni- versite d'Ascq-les-Lille, 1982. 128. Sablonniere P. Bernstein-Bezier methods for the construction of bi-variate spline approximants, CAGD, mai 1985. 129. Sawai S. et al. Development of SMDS: SONY Mold Design System, PROLAMAT, 1985. 130. SCAGD: Surfaces in Computer Aided Geometric Design. Actes de la conference d'Oberwolfach an 1982, North Holland. 131. Schoenberg I. J. On variation diminishing approximation methods. On numerical approximation, University of Wisconsin Press, 1959. 132. Schoenberg I. J. Smoothing operators and their generating functions, Bulletin of the american mathematics Society, 1953. 133. Sederberg Т. W. Piecewise algebraic surface patches, CAGD, mai 1985. 134. Shi F. S., Hand K. On the derivation of the Bezier's function, Beijing Institute of Aeronautics and Astronautics, People's Republic of China, 1979. 135. Shu H., Hori S., Mann W. R., Little R. N. The synthesis of sculptured surfaces, PROLAMAT 1969. 136. Sioussiou S. Precede de description de formes par courbes polynomiales, These de doctorat de 1'Universite de Paris, 1982. 137. Skappel E. The Saab sculptured surfaces, system "Formela". Saab-Scania A. B. 138. Steel, 9 november 1964, NC diesinking is here. Are you ready?

139. Sutcliffe D. C., Contouring over rectangular and skewed rectangular grids, MMCGD, 1980.

140. Toth D. L. On Ray Tracing Parametric surfaces. Computer Graphics, July 1985.

141. Varady Т. Operations to integrate free-form surfaces into the BUILD volumetric modeller, PROLAMAT 1985.

Литература 255

142. Veron M. Contribution a 1'etude des surfaces numeriques UNISURF, Conditions de raccordement. These de doctorat de 1'Universite de Nancy, 1973.

143. Wielinga R. F. Constrained interpolation using Bezier curves as a new tool for computer-aided geometric design, CAGD, 1974.

144. Wu J. H. Mathematical foundations of Bezier's curves, Beijing Institute of Aero- nautics and Astronautics, 1979. People's Republic of China.

145. Xiong Z. X. The spline functions of (2n + I)'11 degree, Beijing Institute of Aero- nautics and Astronautics, People's Republic of China.

146. Yajima A., Jonishi H., Tsuda J., Osada N„ MDM. 1: a computer-aided Mold Design and Manufacturing system. Graphics Interface' 82, National Computer Graphics Association of Canada.

147. Farin G. Triangular Bernstein-Bezier patches-CAD, aout 1986.

ПРЕДМЕТНЫЙ УКАЗАТЕЛЬ

Автоматизация поиска минимума 62 Автоматизированное вычисление 7 Автономные вычисления 8 Алгоритм 106 Алгоритмы построения 157—160 — — характеристической сетки 186—195 «Аналитические» зависимости 103—111 — тела 103, 114 Аппроксимация 106, 131 А-устойчивость метода 20

Базовая кривая 136 Барицентрические координаты 18, 19 Булевы комбинации 111 Внешний алгоритм с интерфейсом 8 Внутренний алгоритм 7 Вырожденная клетка 131 Вырожденные поверхности 175 Годографы 153—154, 182 Графическая система 72 Диагональное разбиение 233 Интерактивность 73 Интерполяционная функция Кунса 126 Интерполяция 106, 114 Квадратичная сходимость метода 13 Клетка с точкой излома 232 Коммерческие САПР 6 Ленточная матрица 9 Линейная сходимость метода 13 Линейное преобразование 135 Линии видимого контура 185 Локальная система баз данных 72 Макроинструкции 8 Матрица Хессенберга 11 Метод Бэрстоу 14 — Галеркина 29 — Гаусса—Лежандра 17 — дихотомии 14 — Инаба 121 — Кастельжо 131 — коллокации 26 — конечных разностей 31—36 — — элементов 40—45, 74 — Кунса 125 — моментов 28 — наименьших квадратов 27 — неопределенных множителей 62 — Ньютона 62 — Ньютона—Рафсона 13 — релаксации 62 — Ризенфольда 133 — Рэлея—Ритца 29 — сопряженных градиентов 12 — Файяра 121 — Фергюсона 123 — штрафных функций 62 Методы построения разбиения 74, 77—86 — спуска 62 Невязка 25, 26 Неявная формула Милна 22 — — Эйлера 19 Объемное моделирование 110 Ограничительные линии 112 Операционная система 72 Оптимизация 61 Отображения 201 Параметрические преобразования 197 Передача геометрических данных 8 Плоские кривые 107 Поверхности сшивки 113 Погрешность метода 20 Полиномиальные функции 152, 175 Полиномы Бернштейна 147 Последовательные годографы 185 Принцип возможных перемещений 53 Проекционные методы 25 Пространственные кривые 107 — тела 196 Разбиение 72 Разреженная матрица 12 Предметный указатель 257 Свойства функций Бернштейна 147—150 Система UNISURF 184 Смешанные годографы 183 Сплайны 114 Структура данных 74 Схема Кранка—Никольсона 59 Сшивка 210 Тестирование метода 73 Трехсторонние многоугольники 137 Формула Грина 45 — трапеций 20 — Хемминга 22 Функция Дирака 49 Характеристические сетки 172 Характеристический многоугольник 136, 141 Число обусловленности матрицы 11 Чувствительность решения 61

Эволюционные уравнения 56

Экспериментальные зависимости 103

Явная формула Адамса 22 — Милна 11 — Эйлера 19

Оглавление

Предисловие редактора перевода ................................................. 5

Часть 3. Вычислительные методы

77. Жермен-Лакур

Глава 1. Введение ....................................................................... 6

1.1. Общие положения ........................................................ 6

1.2. Вычисления и САПР .................................................... 7

Глава 2. Численные методы ....................................................... 9

2.1. Алгебраические уравнения ............................................. 9

2.2. Численное интегрирование ............................................ 16

2.3. Дифференциальные уравнения с начальными условиями . 19

2.4. Проекционные методы ................................................. 25

2.5. Метод конечных разностей ........................................... 31

2.6. Дискретные системы .................................................... 36

2.7. Метод конечных элементов для одной переменной ........ 40

2.8. Метод конечных элементов для двумерного оператора Лапласа ....................................................................... 45

2.9. Уравнения с контурными интегралами .......................... 48

2.10. Статическая линейная трехмерная задача сопротивления материалов .................................................................. 51

2.11. Особые случаи ............................................................. 55

2.12. Нелинейная механика .................................................... 56

2.13. Эволюционные уравнения ............................................. 56

2.14. Оптимизация ................................................................ 61

Оглавление 259

Глава 3. Реализация численных методов в САПР ........................ 65

3.1. Построение разбиения .................................................. 65

3.2. Визуализация результатов ............................................. 66

3.3. Программирование ....................................................... 72

3.4. Связи .......................................................................... 73

3.5. Тестирование ................................................................ 73

Глава 4. Разбиения для метода конечных элементов..................... 74

4.1. Постановка задачи ....................................................... 74

4.2. Разбиение и структура данных ..................................... 75

4.3. Методы построения разбиения и методология проектирования ............................................................ 77

4.4. Методы постороения разбиения .................................... 80

4.5. Визуализация. Контроль. Примеры ............................... 86

4.6. Замечания о полном решении задачи ............................ 88

Заключение ................................................................................. 88

Часть 4. Геометрические методы

77. Безье

Предисловие ......................................................................... 96

Глава 1. Введение ....................................................................... 100

1.1. Задачи САПР ............................................................... 100

1.2. Геометрия, статика и динамика .................................... 100

1.3. Производство: подготовка и организация ...................... 101

Глава 2. Кривые и поверхности .................................................. 102

2.1. Основные проблемы. Определяемые зависимости ........... 102

2.2. Схемы получения решений ............................................ 105

2.3. Условия, учитываемые при поиске решения ................... 107

Глава 3. Системы ....................................................................... 110

3.1. Объемное моделирование .............................................. 110

3.2. Интерполяция .............................................................. 114

3.3. Аппроксимация ............................................................ 131

260 Оглавление

Глава 4. Система UNISURF ........................................................ 135

4.1. Кривые ........................................................................ 135

4.2. Поверхности ................................................................ 172

4.3. Пространственные тела ................................................

4.4. Параметрические преобразования .................................. 197

4.5. Сшивка ........................................................................ 210

Заключение ................................................................................. 238

Литература ................................................................................. 249 Предметный указатель ................................................................ 256

Научное издание

Пьер Жермен-Лакур, П. Л. Жорж, Ф. Пистр, Пьер Безье

МАТЕМАТИКА И САПР

В 2-х книгах

Книга 2

Заведующий редакцией д-р техн. наук А. Л. Щерс Зам. заведующего редакцией Э. Н. Вадиков Ст. научный редактор И. М. Андреева Мл. научный редактор Н. И. Сивилева Художник Л. М. Муратова Художественный редактор Н. И. Иванов Технические редакторы Л. А. Тихомирова, М. А. Анциферова Корректор Р. Л. Вибке ИБ № 6883

Подписано к печати 21.08.89. Формат 60 х 90Уч, Бумага офсетная № 1 Гарнитура тайме. Печать офсетная. Объем 8,25 бум. л. Усл. печ. л. 16,50. Усл. кр.-отт. 30,80. Уч.-изд. л. 17,19. Изд. № 6/6143. Тираж 23000 экз. Зак.1172. Цена I р. 40 к. Издательство «Мир»

В/О «Совэкспорткнига» Государственного комитета СССР по делам издательств, полиграфии и книжной торговли.

129820, ГСП, Москва, 1-й Рижский пер., 2.

Набрано в Межиздательском фотонаборном центре издательства «Мир»

Можайский полиграфкомбинат В/О «Совэкспорткнига» Государственного комитета СССР по делам издательств, полиграфии и книжной торговли. 143200, Можайск, ул. Мира, 93.

УВАЖАЕМЫЙ ЧИТАТЕЛЬ!

Ваши замечания о содержании книги, ее оформлении, качестве перевода и другие просим присылать по адресу: 129820, Москва, И-110, ГСП, 1-Рижский пер., 2, изд-во «Мир».

ИЗДАТЕЛЬСТВО «МИР» ВЫПУСТИЛО КНИГУ

Математика и САПР: В 2-х кн. Кн. 1. Пер. с франц./ Шенен П., Коснар М., Гардан И. и др. — 1988. — 204 с., 1р. Юк.

Книга французских специалистов посвящена математическим основам методов графического построения кривых и поверхностей, используемых в САПР. Излагаются методы интерполяции, аппроксимации, сглаживания, метод конечных элементов и метод конечных разностей.

Для специалистов в области САПР и студентов высших учебных за- ведений.

ИЗДАТЕЛЬСТВО «МИР» ВЫПУСТИЛО КНИГУ

Кулон Ж.-Л., Сабоннадьер Ж.-К., САПР в электротехнике. Пер. с франц., 1988. — 208 с. 1 р.

В книге французских специалистов подробно изложены методы теории конечных элементов для расчета основных параметров электротехнических устройств, контакторов, полупроводниковых преобразователей, управления электрическими сетями, предполагающие интерактивный режим работы проектировщика с ЭВМ. Значительное внимание уделено созданию аппарат- ных и программных средств САПР применительно к предприятиям электро- технической промышленности. Приведены перспективные направления раз- вития САПР.

Для инженеров-электриков, а также аспирантов и студентов соответ- ствующих специальностей вузов.