Вычислительная неустойчивость алгоритма метода начальных функций Текст научной статьи по специальности «Математика»

ВЕСТНИК САНКТ-ПЕТЕРБУРГСКОГО УНИВЕРСИТЕТА

Сер. 10. 2010. Вып. 4

УДК 519.6+539.3 А. В. Матросов

ВЫЧИСЛИТЕЛЬНАЯ НЕУСТОЙЧИВОСТЬ АЛГОРИТМА МЕТОДА НАЧАЛЬНЫХ ФУНКЦИЙ

Введение. В работе [1] представлен алгоритм построения общего решения на основе суперпозиции двух решений, полученных методом начальных функций (МНФ), для анизотропного упругого прямоугольника, находящегося в условиях плоской задачи. Это решение применяется в алгоритме расчета тел сложной конфигурации, составленных из анизотропных упругих прямоугольных областей [2]. Известно, что реализация тригонометрическими рядами алгоритма МНФ обладает вычислительной неустойчивостью, которая проявляется при вычислениях с большими гармониками. В [3] разработаны две вычислительно-устойчивые модификации алгоритма МНФ, основанные на представлении рассчитываемой прямоугольной области в виде слоистой структуры. Однако использовать их в процедуре построения общего решения для прямоугольной области не удалось в связи с их численной ориентацией для традиционного алгоритма МНФ, при котором точно удовлетворяются граничные условия на двух противоположных гранях, тогда как на двух других граничные условия диктуются выбранными для представления начальных функций тригонометрическими функциями.

Вычислительная неустойчивость алгоритма МНФ проявляется тогда, когда расчеты выполняются на подмножестве вещественных чисел с ограниченной длиной мантиссы. В традиционных системах программирования (FORTRAN, С, СН—Ь, Java и др.) используются числа с обычной (8 десятичных цифр в мантиссе) или удвоенной (16 десятичных цифр в мантиссе) точностью. Для вычислений с мантиссой большей длины следует применять либо функции дополнительной библиотеки (FORTRAN, С), либо специальные встроенные пакеты с классами (Java), реализующими подобные вычисления. С появлением систем аналитических вычислений (Maple, Mathematica и др.) появилась возможность без специального перепрограммирования выполнять вычисления с любым количеством цифр в мантиссе, всего лишь меняя значение определенной системной переменной (в Maple, например, Digits).

Для использования МНФ в универсальном алгоритме построения общего решения для упругой анизотропной прямоугольной области следует быть уверенным, что все вычисления выполняются «точно», без погрешностей, обусловленных вычислительной неустойчивостью. Для этого, как и в исследовании по вычислительной неустойчивости метода начальных параметров [4], была проведена серия вычислительных экспериментов для определения при заданной длине мантиссы максимальной гармоники,

Матросов Александр Васильевич — кандидат технических наук, доцент кафедры информационных систем факультета прикладной математики—процессов управления Санкт-Петербургского государственного университета. Количество опубликованных работ: 51. Научные направления: численно-аналитические алгоритмы механики деформируемого твердого тела, вычислительная устойчивость алгоритмов. E-mail: [email protected].

© А. В. Матросов, 2010

при которой вычисления в соответствии с алгоритмом МНФ остаются устойчивыми при выборе тригонометрических функций в качестве начальных.

Математическая модель. Основное соотношение МНФ для линейно-упругого анизотропного тела, находящегося в условиях плоской задачи, в прямоугольной декартовой системе координат может быть получено в следующем виде [1]:

U = LU0. (1)

Здесь U = {u, v, ax, ay, тху} - вектор-столбец перемещений и напряжений в произвольной точке упругого тела, U0 = {u0,v0,aXX,т°у} - вектор-столбец начальных функций, определенных на грани x = 0, L = [Lij] (i = 1,..., 5, j = 1,..., 4) - матрица операторов МНФ вида

Ж

Lij = E k (Aps, dy) xk. (2)

k=0

Коэффициенты ij (Aps,dy) зависят от оператора дифференцирования ду по переменной y и упругих постоянных анизотропного тела Aps из обобщенного закона Гука

а = H £, (3)

где а = {ах,ау,тху} - вектор-столбец компонентов тензора напряжений, е = {ex,ey,

Г All Ai2 Ai6 "

Yxy} - вектор-столбец компонентов тензора деформаций, а H = A12 A22 A26

Ai6 A26 A66 матрица упругих постоянных материала [5].

Соотношение (1) получается из представления решения уравнений равновесия Ламе для плоской задачи

Wu = 0 (4)

через две функции F = {F1 (x, y), F2 (x, y)} следующим образом:

u = W F. (5)

Функции Fi удовлетворяют одному и тому же дифференциальному уравнению четвертого порядка:

(det W) Fi = 0, i = 1, 2. (6)

В (4) u = {u (x, y), v (x, y)} - вектор-столбец перемещений вдоль координатных осей, AiidX + 2Ai6dxdy + A66ду Ai6дХ + (Ai2 + A66) дхду + A26d^ Ai6дХ + (Ai2 + A66) дxдy + A26дХ A 66 дХ + 2A26дxдy + A22 дy рица операторов Ламе, в которой символы дx и дy обозначают операторы диф цирования соответственно по переменным x и y, а константы Aij являются упругими константами тела в случае плоской деформации либо выражены через них в случае плоского напряженного состояния; det W в (6) обозначает определитель матрицы W, элементы которой Wj являются алгебраическими дополнениями элементов Wji матрицы W.

W

мат-ферен-

Компоненты тензора напряжений а представляются через введенные две функции Е^(х, у) с использованием обобщенного закона Гука (3), соотношения Коши е = Си, " дх 0

и представления (5) таким образом:

С

0 ду

ду дх

а = ЯР, (7)

здесь Я = Н С\¥ - матрица дифференциальных операторов.

Уравнение (6) решаем, используя символический способ составления решения дифференциального уравнения в частных производных [6]. Полагая в нем оператор дифференцирования ду символьной константой, его можно рассматривать как обыкновенное дифференциальное уравнение относительно функции, зависящей от одной независимой переменной х:

(Бд + В±д1 + В2дХ + Б1дх + Бу) ВД = 0, (8)

где коэффициенты Б^, ] = 0,...,4, являются дифференциальными операторами:

Бо = ЛпЛее - А216,

= 2(лЦл16 - ЛиЛ1б) ду, Б1 = [ЛцЛи - Л21 + 2 (Л16Л26 - ЛиЛее)] ду, Бз = 2(Л16Л1 у - ЛцЛу6) д3, Б4 = (ЛцЛ66 - Л26) дУу.

Функции ¥1 будем искать в виде степенных рядов

ОО к

^ = (9)

к=0

с неизвестными коэффициентами /к, которые определяются из рекуррентного соотношения четвертого порядка

ЭД+у + Б 1/1+з + Б2/+2 + Б1/1+1 + Б у/1 =0, г = 1, 2, к = 0,...ж,

получаемого после подстановки степенных рядов (9) в дифференциальное уравнение (8).

Все неизвестные коэффициенты / можно выразить через первые четыре неизвестных коэффициента /0, /{, /2 и /3, являющихся соответственно значениями самой функции ^ и ее первых трех производных на линии х = 0. Следовательно, и сами функции ^ могут быть представлены как линейная комбинация своих начальных значений [7]

= мо/о + м/ + м2 /2 + мз/3, (10)

где Мгч = т\к (Лр8,ду) хк (д = 0,..., 3) есть степенные ряды с коэффициентами

к=0

тгдк, зависящими от упругих констант материала и символа дифференцирования ду.

Начальные значения являются произвольными функциями переменной у, не имеющими никакого механического смысла для задач теории упругости. Однако их можно выразить через компоненты напряженно-деформированного состояния (НДС)

и0 (y) = и (0,y), v0 (y) = v (O, y), a<x (y) = ax (0,y) и т°у (y) = тХу (0,y), отредететш-ми на начальной линии x = 0. Эти компоненты называются начальными функциями, давшими название методу. Подстановка (10) в (5) и (7) и последующее вычисление с помощью полученных представлений компонентов НДС на линии x = 0 приводит к выражению начальных функций через начальные значения функций Fi в виде

u0 = Wn ¿ m\o f¡+Wu ¿ m\ofl,

q=0 q=0

v0 = W21 m\o f¡+W22 ¿ m/

q3=0 Г (11)

a0x = Su E m^o/l+Si^ m\of¡,

q=0 q=0

3 3

тХу = S3lE m1o/l+S32T. m1o/l-

y q=0 q q q=0 q q

Выражения (11) можно трактовать как систему линейных алгебраических уравнений относительно начальных значений функций Fi. Система не доопределена (восемь неизвестных /¿, /2., /3, /§, /2, /° и /3 при четырех уравнениях), а поэтому для ее разрешения можно любые четыре неизвестных выбрать произвольным образом. Анализ матрицы системы показывает, что следует положить равными нулю начальные значения /О = // = /0 = /2 =0, так как в этом случае определитель системы (11) не будет содержать символ оператора дифференцирования ду, а будет выражен только через упругие постоянные Aij. Это условие необходимо, так как для получения решения системы необходимо выполнять деление на определитель системы. Выразив неизвестные величины /21, /31, /22 и /32 из системы (11) через начальные функции, получим представление компонентов вектора НДС U через вектор начальных функций U0 в виде (1), снова используя соотношения (5) и (7).

Если известны все компоненты вектора начальных функций, то соотношение (1) позволяет вычислить все компоненты НДС в любой точке упругого тела. Однако обычно не все начальные функции известны, так как для любой плоской статической задачи на каждой из граней упругого тела известны обычно два компонента НДС. Таким образом, на начальной плоскости, если она совмещена с одной из граней тела, из условий задачи обычно известны два компонента. Остальные два можно определить из граничных условий на противоположной грани x = h. Например, если на начальной линии x = 0 заданы напряжения (у),т°у (y), а на грани x = h - напряжения аХХ (у),тХХу (y), то для определения недостающих начальных функций и0 (y) и v0 (y) можно получить систему двух дифференциальных уравнений бесконечного порядка, вычислив значения соответствующих компонентов напряжений на линии x = h на основе соотношения (1):

L31 L, и0 + L32 Lh v0 = аХ - L33 \x=h а° - L34 \x=h т°у,

L51 \x=h u0 + L52 \x=h v0 = тХХу - L53 \x=h аХ - L54 |

0

(12)

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

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

Для тел с произвольной степенью анизотропии операторы матрицы L в соотношении (1) замкнуты в случае плоской задачи относительно линейной комбинации тригонометрических функций sm _ sin (mny/a) и cm _ cos (mny/a). Это означает, что если все начальные функции будут представлены в виде суммы тригонометрических рядов по синусам и косинусам

тт0 _ Y^ Tj0>sm i T Т0,cos

mm

m=0

u0,sin _ fb1,sin „a ^2,sinca b3,sin „a b4?sinca 1

(13)

m X m sm, m cm, m sm, m cm J ,

u0,cos _ Ihlj cos ca ь2,cos„a b3, cos ca k4,cos„a \

m X m cm, m sm, m cm, m sm),

то и все компоненты вектора НДС U также будут являться суммами тригонометрических рядов по синусам и косинусам

u _ £ umm+u

cos

mm 0

usin _ Í т 1,sin ^a T2jsiUca т 3,sin „a т 4,sin „a T5,siUc^ X

m \Tm sm, Tm cm, Tm sm, Tm sm, Tm cm J ,

ucos _ í т 1,cos c,a г 2,cos „a т3,cos r.a т 4,cosr.a т 5,cos „a \

m X^m cm, тm „m, тm cm, тm cm, тm „m J '

Коэффициенты Ь^™ и bl¡cos (i _ 1,..., 5) суть степенные ряды по переменной x, вычисляемые как результат воздействия соответствующих операторов начальных функций (2) на соответствующие тригонометрические функции „m и cm согласно соотношению (1) и представлению начальных функций (13) с последующим приведением подобных членов при тригонометрических функциях.

Для ортотропных, трансверсально-изотропных и изотропных тел операторы МНФ

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

00

и косинусам, но и только по синусам (U0 _ £ Umsin) или только по косинусам (U0 _

m=0

£ Umcos). При этом коэффициенты в представлении компонентов вектора НДС U

m=0

записываются намного проще, чем в случае общей анизотропии:

^ 4 ^ 4

тг,sin hj,sinjk,sinxk ji,cos hi,cosjk,cosxk (14)

^m / j / j m ij,m > ^m / j / j um ij,m ^ •

k=0 j=1 k=0 j=1

тг i i 7k,sin 7k, cos

Причем коэффициенты iijm и iijm вычисляются как результат воздействия операторов Ik (Aps,dy) из (2) на тригонометрические функции представления начальных функций по синусам или косинусам (13).

Уравнения типа (12) для вычисления неизвестных начальных функций в этом случае преобразуются в серию систем линейных алгебраических уравнений для нахождения неизвестных коэффициентов в их представлениях тригонометрическими рядами -для каждой гармоники своя система:

V^ jk,sin ь i b i,sin I jk,sin fc I ,2,sin _ ah

Z^ 131,mx \x=h bm + Z^ l32,mx \x=h bm _ axm

k=0 k=0

oo o .

чр , k,sin ь i b3,sin _ V^ jk,sin xk I

Z^ l33,mx \x=h bm Z^ l34,mx \x=h un

k=0 k=0

oo

V^ jk,sin xk I b1,sin i jk,sin b2,sin _ Th _

Z^ l51,mx \x^h bm ' Z^ L52,mx \x=h bm ~ 1 xym k=0 , k=0 , m

oo

j k,sinxk I b3,sin V^ j k,sin xk I b4,sin Z^ l53,mx \x=h bm Z^ l54,mx \x=h bm ■

k=0 k=0

Здесь через ahm и т1^Ут обозначены коэффициенты при m-ных гармониках в разложении в ряды Фурье заданных функций напряжений на нижней грани прямоугольной области.

С увеличением номера гармоники оказывается невозможным получить достоверные значения компонентов НДС в области, отдаленной от начальной плоскости. Это связано как с накоплением ошибок при расчете числовых коэффициентов в рядах (2), так и с характером сходимости самих рядов, приводящим к катастрофической потере верных цифр в получаемой сумме этих рядов [8], используемых в качестве коэффициентов системы (15), которая становится плохо обусловленной с увеличением номера гармоники. Анализ и решение указанной проблемы показаны на примере модельной задачи деформирования изотропного тела.

Модельная задача исследования. Рассмотрим задачу расчета НДС прямоугольного в сечении с размерами h х a (h/a _ 1) изотропного тела (E _ E0, v _ 1/3), находящегося в условиях плоской деформации. По грани x _ 0 (начальная плоскость) к телу приложены усилия a° _ q0 sin (mny/a), r°y _ T0 cos (mny/a) (т0 _ 0), а по грани x _ h - усилия ah _ qh sin (mny/a), r^y _ Th cos (mny/a) (qh _ 0, Th _ q0). Неизвестные перемещения начальной плоскости x _ 0 ищутся в виде u0 _ U0 sin (mny/a), v0 _ V0 cos (mny/a) с неопределенными коэффициентами, которые определяются из решения системы уравнений (15), полагая в ней bmsin = U0, bmsin = V0, bmisin = q0, 6msin = t0, ahm _ qh и т^Ут _ Th. На двух других гранях тела граничные условия продиктованы выбранными тригонометрическими функциями представления компонентов вектора начальных функций.

Вычислительные эксперименты. Расчеты НДС с удвоенной точностью представления вещественных чисел (16 цифр в мантиссе) показывают наличие вычислительной неустойчивости при m ^ 10 (графики безразмерных перемещений uE0/q0h (рис. 1, I, а) в сечениях y _ a/2m и касательных напряжений тхy/q0 (рис. 1, I, б) в сечениях y _ a/m).

Простое увеличение в расчетах мантиссы вещественных чисел до 18 цифр позволяет получить достоверные результаты при m ^ 11 (рис. 1, II). Но при больших значениях гармоник вычисления продолжают оставаться неустойчивыми. Это связано с катастрофической потерей верных цифр в суммах рядов (14), являющихся коэффициентами линейной системы (15), число обусловленности которой растет с увеличением гармоники. Накопление ошибок при расчете сумм степенных рядов связано с характером сходимости последних - их частичные суммы сначала возрастают до чисел больших порядков, а потом убывают до не очень больших значений сумм рядов. Максимальные значения частичных сумм достигают порядка 107 для перемещений и 1014 для напряжений уже при m _ 10, что и служит причиной того, что суммы соответствующих степенных рядов вычисляются с большой погрешностью (не хватает цифр в мантиссе

Рис. 1. Безразмерные перемещения иЕо/доН (а) и касательные напряжения Тху/до (б) для различных гармоник при вычислениях с мантиссой длиной 16 (I) и 18 цифр (II)

для их расчета с необходимой точностью). Более того, суммируемые ряды знакопере-менны, а это приводит к вычислению разности больших близких величин, определенных с большой погрешностью. На рис. 2, а представлены в логарифмической шкале абсолютные значения частичных сумм степенных рядов (14) для вычисления перемещения и в горизонтальном сечении х = к при разных гармониках. Для сравнения на рис. 2, б приведены те же графики, вычисленные точно с использованием мантиссы длиной 90 цифр.

Рис. 2. Графики частичных сумм рядов (14) при вычислениях с мантиссой длиной 16 (а) и 90 цифр (б) для различных гармоник

Из этих же графиков частичных сумм степенных рядов (14) видно, что для вычисления правильного значения их сумм следует учитывать значительное число их членов, особенно при больших значениях гармоник. Так, выполняя вычисления с длиной мантиссы в 90 цифр, но учитывая только 89 членов в указанных рядах, получить достоверные результаты удается только при т = 9, что видно из графиков безразмерных перемещений иЕ0/к/д0 в сечениях у = а/2т и напряжений тху/до в сечениях у = а/2т, представленных на рис. 3. Значения перемещений и напряжений в нижней части плиты при т = 11 становятся столь большими, что пришлось ограничить отображаемый на рисунке их диапазон.

Для выявления предельных значений гармоник, при которых возможен устойчивый вычислительный процесс, была проведена серия вычислительных экспериментов, результаты которой сведены в таблицу. В ней приведены предельные «устойчивые» значения гармоник в зависимости от геометрических размеров пластины (а/к) и количества удерживаемых членов в степенных рядах (14). Вычисления во всех экспериментах выполнялись с мантиссой длиной 120 цифр.

Следует отметить, что вычислительные эксперименты проводились для «толстых» пластин, у которых высота больше ширины. Это связано с тем, что, во-первых, для «тонких» пластин вплоть до практически значимых значений гармоник вычислительная устойчивость или недостаточное количество удерживаемых членов не проявляется (например, для пластины с а/к = 2 при 101 удерживаемом члене в степенных

б

Рис. 3. Безразмерные перемещения иЕо/док (а) и касательные напряжения Тху/до (б) для различных гармоник при вычислениях с мантиссой длиной 90 цифр и удержанием 89 членов в степенных рядах (14) Обозначения см. на рис. 1

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

а/к Количество членов степенных рядов

101 201 301 401 501 601 701

1/1 9 21 33 45 57 69 81

1/2 5 11 17 21 27 33 39

1/3 3 7 11 15 19 23 27

1/4 5 7 11 13 17 19

1/5 3 5 9 11 13 15

1/6 3 5 7 9 11 13

1/7 3 3 5 7 9 11

1/8 1 3 5 7 7 9

1/9 1 3 5 5 7 9

1/10 1 3 3 5 5 7

рядах устойчивые вычисления возможны вплоть до гармоники 23), во-вторых, эксперименты проводились, имея в виду применение МНФ для построения общего решения упругой прямоугольной пластины [1], а в этом случае как раз и важен результат для пластин, у которых высота больше ширины, так как общее решение строится как суперпозиция двух решений МНФ, в которых используются степенные ряды по двум координатам.

Заключение. В работе рассмотрен один из возможных способов преодоления вычислительной неустойчивости алгоритма МНФ - увеличение длины мантиссы при числовых расчетах. Все расчетные схемы метода реализованы в системе аналитических

вычислений Maple, позволяющей легко манипулировать длиной мантиссы. Выявлены причины возникновения вычислительной неустойчивости алгоритма МНФ при его реализации в классе тригонометрических полиномов, связанные с характером сходимости степенных рядов, получаемых после воздействия операторов МНФ на тригонометрические функции, через которые выражаются начальные функции. Определены минимальные значения длин мантисс в представлении вещественных чисел при различном числе удерживаемых членов в степенных рядах решений, при которых вычисления по МНФ дают достоверные результаты для изотропных упругих пластин с разными геометрическими характеристиками.

Литература

1. Матросов А. В. Численно-аналитическое решение граничной задачи деформирования линейно-упругого анизотропного прямоугольника // Вестн. С.-Петерб. ун-та. Сер. 10: Прикладная математика, информатика, процессы управления. 2007. Вып. 2. С. 55—65.

2. Матросов А. В. Численно-аналитический алгоритм решения задач плоской деформации линейно-упругих тел сложной конфигурации // Вестн. С.-Петерб. ун-та. Сер. 10: Прикладная математика, информатика, процессы управления. 2008. Вып. 3. С. 70—84.

3. Galileev S. M., Matrosov A. V. Method of initial functions: stable algorithms in the analysis of thick laminated composite structures // Composite Structures. 1997. Vol. 39, N 3—4. P. 255—262.

4. Матросов А. В. Численно-аналитический алгоритм метода начальных параметров // Вестн. С.-Петерб. ун-та. Сер. 10: Прикладная математика, информатика, процессы управления. 2009. Вып. 2. С. 72-81.

5. Лехницкий С. Г. Теория упругости анизотропного тела. Изд. 2-е. М.: Наука, Гл. ред. физ.-мат. лит-ры, 1977. 416 с.

6. Лурье А. И. Пространственные задачи теории упругости. М.: Гос. изд-во техн.-теор. лит-ры, 1955. 491 с.

7. Galileev S. M., Matrosov A. V. Method of initial functions in the computation of sandwich plates // Intern. Applied Mechanics. 1995. Vol. 31, N 6. P. 469-476.

8. Форсайт Дж., Малькольм М., Моулер К. Машинные методы математических вычислений / пер. с англ. Х. Д. Икрамова. М.: Мир, 1980. 279 с.

Статья рекомендована к печати проф. Ю. М. Далем. Статья принята к печати 10 июня 2010 г.