авторефераты диссертаций БЕСПЛАТНАЯ БИБЛИОТЕКА РОССИИ

КОНФЕРЕНЦИИ, КНИГИ, ПОСОБИЯ, НАУЧНЫЕ ИЗДАНИЯ

<< ГЛАВНАЯ
АГРОИНЖЕНЕРИЯ
АСТРОНОМИЯ
БЕЗОПАСНОСТЬ
БИОЛОГИЯ
ЗЕМЛЯ
ИНФОРМАТИКА
ИСКУССТВОВЕДЕНИЕ
ИСТОРИЯ
КУЛЬТУРОЛОГИЯ
МАШИНОСТРОЕНИЕ
МЕДИЦИНА
МЕТАЛЛУРГИЯ
МЕХАНИКА
ПЕДАГОГИКА
ПОЛИТИКА
ПРИБОРОСТРОЕНИЕ
ПРОДОВОЛЬСТВИЕ
ПСИХОЛОГИЯ
РАДИОТЕХНИКА
СЕЛЬСКОЕ ХОЗЯЙСТВО
СОЦИОЛОГИЯ
СТРОИТЕЛЬСТВО
ТЕХНИЧЕСКИЕ НАУКИ
ТРАНСПОРТ
ФАРМАЦЕВТИКА
ФИЗИКА
ФИЗИОЛОГИЯ
ФИЛОЛОГИЯ
ФИЛОСОФИЯ
ХИМИЯ
ЭКОНОМИКА
ЭЛЕКТРОТЕХНИКА
ЭНЕРГЕТИКА
ЮРИСПРУДЕНЦИЯ
ЯЗЫКОЗНАНИЕ
РАЗНОЕ
КОНТАКТЫ


Pages:   || 2 | 3 | 4 |
-- [ Страница 1 ] --

А. В. Анкилов, П. А. Вельмисов, А. С. Семёнов

АЛГ ОР ИТ МЫ МЕ Т О Д О В

ВЗВЕ Ш Е ННЫ Х НЕВЯЗОК

ДЛЯ РЕШЕНИЯ ЛИНЕЙНЫХ ЗАДАЧ

МАТЕМАТИЧЕСКОЙ ФИЗИКИ

И

ИХ РЕАЛИЗАЦИЯ В СИСТЕМЕ MATHCAD

Учебное пособие

Ульяновск 2006

УДК 519.6 (075)

ББК 22.311 я7

A 67

Рецензенты: Кафедра прикладной математики Ульяновского государственного

университета (зав. кафедрой доктор физико-математических наук,

профессор А. А. Бутов);

Доктор физико-математических наук, проф. УлГУ В. Л. Леонтьев.

Утверждено редакционно-издательским советом университета в качестве учебного пособия.

Анкилов, А. В. Алгоритмы методов взвешенных невязок для решения линейных задач математической физики и их реализация в системе A 67 MathCAD: учебное пособие / А. В. Анкилов, П. А. Вельмисов, А. С. Семёнов. – Ульяновск : УлГТУ, 2006. – 168 с.

ISBN 5-89146-913- Пособие содержит изложение алгоритмов численного решения некоторых линейных краевых или начально-краевых задач математической физики.

Приведены постановки лабораторных работ, с использованием математически ориентированного пакета M athCAD. Даны примеры.

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

Работа выполнена на кафедре «Высшая математика» УлГТУ. Печатается в авторской редакции.

УДК 519.6 (075) ББК 22.311 я Учебное издание АНКИЛОВ Андрей Владимирович ВЕЛЬМ ИСОВ Петр Александрович СЕМ ЁНОВ Алексей Степанович Алгоритмы методов взвешенных невязок для решения линейных задач математической физики и их реализация в системе MathCAD Учебное пособие Подписано в печать 20.12.2006. Формат 60 84/16.

Бумага офсетная. Усл. печ. л. 10,00.

Тираж 200 экз. Заказ Ульяновский государственный технический университет 432027, Ульяновск, ул. Северный Венец, 32.

Типография УлГТУ, 432027, Ульяновск, ул. Северный Венец, 32.

© А. В. Анкилов, П. А. Вельмисов, А. С. Семёнов, © Оформление. УлГТУ, ISBN 5-89146-913- ОГЛАВЛЕНИЕ Предисловие ……………………………………………………………… Программы по лабораторным работам ….…………………..……….. 1. Математическое моделирование физических задач …………….. 1.1. Вывод уравнений одномерной теплопроводнос ти ………..….. 1.2. Постановка начально-краевой задачи одномерной стационарной теплопроводнос ти …………………………………… 1.3. Постановка начально-краевой задачи одномерной нестационарной теплопроводности ………………………………… 1.4. Постановка краевых задач двухмерной стационарной теплопроводности …………………………………………………… 1.5. Вывод уравнений поперечных колебаний струны …………… 1.6. Вывод уравнения продольных и крутильных колебаний cтержня ……………………………………………………………….

1.7. Постановка статических краевых задач для струны и стержня 1.8. Краевые задачи в теории колебаний струн и стержней ……… 2. Решение краевой задачи для линейного обыкновенного дифференциального уравнения второго порядка методом Галеркина ……………………………………………………. 2.1. Подс тановка задачи и алгоритм метода ………………………. 2.2. Построение систем пробных и поверочных функций ………. 2.3. Задание к лабораторной работе ………………………………. 2.4. Выполнение работы в компьютерном классе ……………….. 2.5. Порядок выполнения лабораторной работы …………………. 2.6. Тестирующий пример ………………………………………….. 2.7. Вопросы для самоконтроля ……………………………………. 3. Решение краевой задачи для линейного обыкновенного дифференциального уравнения второго порядка вариационным методом Ритца ……………………………………… 3.1. Подс тановка задачи и алгоритм метода ……………………… 3.2. Построение систем пробных и поверочных функций ………. 3.3. Задание к лабораторной работе ……………………………….. 3.4. Выполнение работы в компьютерном классе ……………….. 3.5. Порядок выполнения лабораторной работы …………………. 3.6. Тестирующий пример …………………………………………. 3.7. Вопросы для самоконтроля …………………………………… 4. Решение краевой задачи для линейного обыкновенного дифференциального уравнения второго порядка интегральным методом наименьших квадратов …………………… 4.1. Подс тановка задачи и алгоритм метода ……………………… 4.2. Задание к лабораторной работе ……………………………….. 4.3. Выполнение работы в компьютерном классе ………………… 4.4. Порядок выполнения лабораторной работы …………………. 4.5. Тестирующий пример ………………………………………….. 4.6. Вопросы для самоконтроля ……………………………………. 5. Решение начально-краевой задачи для одномерного параболического уравнения методом Галеркина..…………………. 5.1. Подс тановка задачи и алгоритм метода ………………………. 5.2. О построении функции u 0 ( x, t) ………………………………… 5.3. Задание к лабораторной работе ……………………………….. 5.4. Выполнение работы в компьютерном классе ……………….. 5.5. Порядок выполнения лабораторной работы …………………. 5.6. Тестирующий пример …………………………………………. 5.7. Вопросы для самоконтроля …………………………………… 6. Решение начально-краевой задачи для одномерного гиперболического уравнения методом Галеркина..……………….. 6.1. Подс тановка задачи и алгоритм метода ……………………… 6.2. Задание к лабораторной работе …………………………….…. 6.3. Выполнение работы в компьютерном классе ……………….. 6.4. Порядок выполнения лабораторной работы …………………. 6.5. Тестирующий пример …………………………………………. 6.6. Вопросы для самоконтроля ……………………………………. 7. Решение первой краевой задачи для двухмерного эллиптического уравнения методом Галеркина..………………….. 7.1. Подс тановка задачи и алгоритм метода ………………………. 7.2. Задание к лабораторной работе ……………………………….. 7.3. Выполнение работы в компьютерном классе ……………….. 7.4. Порядок выполнения лабораторной работы ………………….. 7.5. Тестирующий пример ………………………………………….. 7.6. Вопросы для самоконтроля ……………………………………. 8. Прикладной математический пакет «MathCAD»..……………… 8.1. О программе …………………………………….……………… 8.2. Основные понятия и функции ………………………………… 8.3. Операторы математического анализа ………………………… 8.4. Функции и операторы матриц ………………………………… 8.5. Создание декартовых графиков на плоскости и в пространстве ………………………………………………………. 8.6. Программные блоки …………………………………………… Приложение А ……………………………………………………………. Приложение Б ……………………………………………………………. Приложение В ……………………………………………………………. Приложение Г ……………………………………………………………. Заключение ………………………………………………………………. Библиографический список ………..…………………………………… Предисловие Целый ряд современных методов, предназначенных для решения самых разнообразных задач математической физики, базируется на идеях ученых Б. В. Галеркина и В. Ритца. К этим методам относятся, например, методы взвешенных невязок и вариационные методы [1,2].

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

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

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

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

Программы по лабораторным работам Для выполнения лабораторных работ разработаны в прикладной системе MathCAD 2000 professional файлы ODU.mcd, Parab.mcd, Giperb.mcd, Ellipt.mcd.

В начале каждого файла расписано задание на лабораторную работу, а в конце приведены варианты заданий, номера которых студенты получают от преподавателя. Тексты программ подробно описывают все дейс твия и операции. Для удобства изучения сами тексты с разобранными тестирующими примерами приведены в приложениях А, Б, В, Г данного пособия.

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

1. Математическое моделирование физических задач 1.1. Вывод уравнений одномерной теплопроводности Пусть дано материальное тело, расположенное между точками x = a и x = b оси Ox, продольный размер которого значительно превосходит размеры поперечного сечения, например, тонкий стержень, длинный трубопровод и т. д.

В дальнейшем будем называть это тело стержнем. Будем считать площадь S (x ) поперечного сечения (перпендикулярного оси Ox ) настолько малой, что всем точкам одного сечения в момент времени t можно приписать одну и ту же температуру u ( x, t). Будем считать, что стержень теплоизолирован вдоль боковой поверхности, а внутри стержня нет источников или стоков (поглотителей) тепла.

Рассмотрим элемент с тержня между его сечениями с абсциссами x и x + dx. Найдем количес тво тепла, которое накапливается в элементе за время dt. Согласно закону Фурье интенсивность q( x, t ) теплового потока в сечении x определяется выражением:

u ( x, t) q( x, t ) = K ( x ), x где K (x) – коэффициент теплопроводнос ти (K ( x) 0). Тогда разность dQ между количеством тепла, вошедшим в элемент через сечение x и вышедшим через сечение x + dx за время dt, будет равна:

u ( x, t) u ( x + dx, t) dQ = K (x )S ( x ) dt K ( x + dx )S (x + dx) dt.

x x Используя формулу Тейлора первого порядка с остаточным членом в форме Пеано для функций K ( x + dx ), S ( x + dx ), u x ( x + dx, t ), имеем u ( x, t ) dt + (K ( x) + K ( x)dx + o(dx ))(S ( x ) + S ( x )dx + o(dx )) dQ = K (x )S (x ) x (1.1) u( x, t ) u 2 ( x, t ) ( x, t) x + x 2 dx + o(dx) dt = x K ( x) S ( x) x dxdt + o(dxdt ).

Напомним, что символом o(x) обозначается величина бесконечно малая более высокого порядка, чем x.

С другой стороны, за счет притока тепла температура в элементе изменяется, и количество тепла dQ, поглощаемое элементом за время dt, равно x + dx C (v )S (v ) (v)(u (v,t + dt ) u (v,t ))dv, dQ = x где C (x ) – теплоемкость;

(x) – объемная плотность вещества стержня (C ( x) 0, (x ) 0 ).

Откуда, на основании теоремы о среднем для определяемого интеграла, получаем равенство dQ = c( ) ( )S ( )(u (, t + dt ) u (, t ))dx, x x + dx, которое при помощи теоремы Лагранжа о конечных приращениях преобразуется к виду u(, t + dt ) dQ = c( ) ( )S ( ) (0,1).

dtdx, (1.2) t Приравнивая, на основании закона сохранения энергии, выражения (1.1), (1.2) и осуществляя предельный переход при dt 0, получаем одномерное уравнение теплопроводности в виде u u KS = C S.

x x t Предположим теперь, что внутри стержня происходит выделение или поглощение тепла (это имеет мес то, например, при прохождении по телу электрического тока или вследствие происходящих в нем химических реакций).

Тогда количество тепла, накопленное в элементе с тержня за время dt за счет внутренних источников, будет равно dQ o = SF (x, t )dxdt, где F ( x, t ) – плотнос ть тепловых источников внутри стержня. Уравнение теплопроводности с учетом внутренних источников тепла принимает вид dQ = dQ + dQ o или u u C S = KS + SF.

t x x Предположим далее, что на боковой поверхности стержня происходит теплообмен с окружающей средой. Тогда тепловой поток, проходящий за время dt через боковую поверхность элемента, согласно закону Ньютона пропорционален разнос ти температур поверхности тела и окружающей среды и определяется выражением dQ* = * ( x)[u (x, t ) T ( x, t )]dxdt, где T (x, t ) – температура внешней среды;

* ( x ) – коэффициент теплообмена, зависящий от свойств материала стержня и внешней среды, режима взаимодействия (условий контакта) стержня с внешней средой, а также от геометрических характеристик поперечного сечения. Уравнение теплопроводности с учетом внутренних источников тепла и теплообмена на боковой поверхности имеет вид dQ = dQ + dQ* + dQ o или u u * (u T ) + SF (x, t ).

C S = KS (1.4) t x x Заметим, что если тепло распространяется в жидкости, которая движется со скоростью V ( x, t ) параллельно оси x, то уравнение теплопроводности запишется следующим образом:

u u u C S + V = KS * (u T ) + SF (x, t ).

t x x x Если тело однородно, т. е. С,, K – постоянные, и площадь сечения S постоянна, то уравнение (1.4) можно записать в виде (u T ) F (x, t ) u 2 u =a * +, CS C t x где a 2 = K / C – коэффициент температуропроводнос ти.

Аналогично уравнению (1.4) выводится уравнение, описывающее процесс распространения тепла в трехмерных телах u = div[K ( x, y, z ) grad u] + F ( x, y, z, t ) C t или в развернутой форме u u u u C = K + K + K + F ( x, y, z, t ). (1.5) t x x y y y y Для однородных тел уравнение (1.5) удобно представить в виде 2 u 2 u 2 u F ( x, y, z, t) u K =a 2 + 2 + 2 +, a2 =.

x z C C t y Для двухмерных тепловых полей в плас тинах, тонких плитах уравнение (1.5) примет вид u u u C = K + K + F ( x, y, t ) (1.6) t x x y y или для однородных плас тин 2 u 2 u F (x, y, t ) u K =a 2 + 2 +, a2 =.

x y C C t 1.2. Постановка краевой задачи одномерной стационарной теплопроводности Согласно (1.4) стационарное (установившееся во времени) распределение теплового поля в стержне постоянного поперечного сечения (S ( x) = const ) описывается уравнением L[y ] = [K ( x) y ] (x ) y = g (x ). (1.7) В (1.7) введены обозначения y (x) = u (x), ( x) = * (x ) / S, g ( x) = F ( x) ( x)T ( x).

Перечислим основные типы граничных условий (на примере левого конца стержня при x = a ).

а) Известная температура при x = a : y (a) = Ta.

б) Задана интенсивность теплового потока через торцевое сечение x = a :

K a y (a ) = q a, K a = K (a ). В частности, если стержень теплоизолирован при x = a, то y (a ) = 0.

в) На конце x = a имеет место теплообмен с окружающей средой известной температуры Ta :

K a y (a) = a [y(a) Ta ], K a = K(a).

Здесь a – коэффициент теплообмена на конце x = a. Последнее условие (условие Ньютона) означает, что тепловой поток, передаваемый в единицу времени с единицы площади поверхнос ти в окружающую среду, пропорционален разности температур поверхности тела и окружающей среды.

Аналогичные краевые условия могут быть заданы и на правом конце стержня при x = b. Например, условие теплообмена при x = b имеет вид K b y (b) = b ( y (b) Tb ).

В таблице 1.1 приведены возможные варианты краевых условий для определения стационарного распределения температуры в стержне согласно уравнению (1.7).

Напомним еще раз используемые в таблице 1.1 обозначения:

K a = K (a), K b = K (b) – коэффициенты теплопроводности;

а, b – коэффициенты теплообмена на левом и правом концах стержня соответс твенно;

Tа, Tb – температуры, которые поддерживаются на концах с тержня при x = a и при x = b ;

qа, qb – интенсивности тепловых потоков при x = a и при x = b.

Очевидно, что все приведенные в таблице 1.1 варианты краевых условий можно записать в виде a0 y( a) + a1 y (a) = a (1.8) b0 y (b) + b1 y (b ) = b2, при соответствующем выборе значений коэффициентов ai, bi.

Например, для первого варианта условий из таблицы 1.1 имеем a0 = 1, a1 = 0, a 2 = Ta, b0 = 1, b2 = 0, b2 = Tb ;

а для девятого – a0 = a, a1 = K a, a2 = a Ta, b0 = a, b1 = K b, b2 = bTb ;

Таким образом, математическая задача одномерной стационарной теплопроводности формулируется следующим образом: требуется найти функцию y (x), удовлетворяющую на отрезке [a, b] обыкновенному линейному дифференциальному уравнению (1.7), а на концах отрезка – граничным условиям (1.8).

Таблица 1. Варианты краевых условий для уравнения (1.5) x x=a x=b № y = Ta y = Tb y = Ta K b y = q b K b y = b ( y Tb ) y = Ta y = Tb K a y = qa K a y = qa K b y = q b K b y = b ( y Tb ) K a y = qa K a y = a ( y Ta ) y = Tb K b y = qb K a y = a ( y Ta ) K b y = b ( y Tb ) K a y = a ( y Ta ) 1.3. Постановка начально-краевой задачи одномерной нестационарной теплопроводности В пункте 1.2 рассмотрена краевая задача для одномерного с тационарного уравнения теплопроводности (1.7), которая представляет собой краевую задачу для обыкновенного дифференциального уравнения второго порядка. В случае нестационарной теплопроводности к краевым (граничным) условиям (1.8) добавляется начальное условие в некоторый начальный момент времени t = t (обычно t = 0 ) u ( x, t0 ) = ( x ), (1.9) и говорят, что задана начально-краевая задача для уравнения параболического типа (1.4).

1.4. Постановка краевых задач двухмерной стационарной теплопроводности Согласно (1.6) стационарное (установившееся во времени) распределение теплового поля в пластине описывается уравнением u u K + K = F ( x, y ). (1.10) x x y y При решении краевых задач для уравнения эллиптического типа (1.10) наиболее часто используются три типа краевых условий.

а) Краевая задача с граничными условиями первого рода (первая краевая задача).

Требуется найти решение уравнения (1.10) в некоторой области D, принимающее на границе этой области заданные значения. Т. е. нужно найти (применительно к рассматриваемой задаче) стачионарное распределение температуры внутри области, если задана температура на границе этой области u = g (x, y ). (1.11) D Здесь D – граница области D, g (x, y ) – известная функция.

б) Краевая задача с граничными условиями второго рода (вторая краевая задача).

Требуется найти решение уравнения (1.10) в некоторой области, на u границе которой задана внешняя нормальная производная (т. е. на границе n задана интенсивность теплового потока).

u u K = q (x, y ) или + q ( x, y) = 0, (1.12) n D n D q где q =. Здесь D – граница области D, q( x, y ) – интенсивность теплового K потока. При этом, если q 0, то тепловой поток направлен наружу, а если q 0, то тепловой поток направлен внутрь облас ти. При q = 0 имеем условие u =0.

теплоизоляции n в) Краевая задача с граничными условиями третьего рода (третья краевая задача).

Требуется найти решение уравнения (1.9) в некоторой области, которое удовлетворяет на границе условию ( ) ( ) u u = u T или = u T, K (1.13) n D n D D D где =. Здесь D – граница области D, на которой задан теплообмен с K окружающей средой, температура которой равна T ;

– коэффициент теплообмена.

Если на различных частях границы D заданы условия различного рода, то такие условия и соответствующие им задачи называют смешанными.

1.5. Вывод уравнений поперечных колебаний струны Рассмотрим тонкую гибкую упругую нить (струну), которая в положении равновесия занимает отрезок [a, b] оси Ox и концы которой закреплены.

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

Рассмотрим элемент струны между точками x и x + dx (рис. 1.1) и обозначим смещение точек струны через u ( x, t ), а длину элемента струны через ds.

Тогда x +dx u ds = 1 + dx, x x откуда, предполагая смещение струны u ( x, t ) малыми настолько, что u 1, (1.14) x получаем ds dx, т. е. в пределах принятой точности удлинения участков струны в процессе колебаний не происходит. Следовательно, согласно закону Гука величина натяжения в каждой точке струны не меняется со временем и является функцией только x, т. е. T = T (x).

Рис. 1.1. Иллюстрация к выводу уравнения колебаний струны Запишем условия динамического равновесия элемента струны, на который действуют в плоскости Oxu силы натяжения T1 = T (x ), T2 = T (x + dx ), внешняя распределенная по длине дуги с линейной плотнос тью F ( x, t ) поперечная сила и сила инерции, направленная вдоль оси Ou.

Проектируя силы на ось Ox, получаем T (x + dx ) cos( 2 ) T ( x) cos( 1 ) = 0. (1.15) Так как, согласно тождествам тригонометрии и геометрического смысла производной, 1 cos( 2 ) = =, 1 + tg 2 ( 2 ) u (x + dx, t ) 1+ x (1.16) 1 cos(1 ) = =, 1 + tg ( 1 ) 2 u ( x, t) 1+ x то, учитывая условие (1.15), из (1.16) получим T (x + dx ) = T ( x). Откуда, в силу произвольности выбора точек x и x + dx, следует, что величина натяжения не зависит и от x, т. е. является постоянной, T (x ) = T0 = const.

Проектируя теперь все силы на ось Ou, получаем 2u ( z,t ) x +dx x+ dx (z ) dz = T0 sin( 2 ) T0 sin( 1 ) + F ( z, t )dz, (1.17) t x x где (x ) – линейная плотность струны.

Аналогично формулам (1.16) устанавливаем u( x + dx, t ) tg( 2 ) x sin( 2 ) = =, 1 + tg ( 2 ) 2 u ( x + dx, t ) 1+ x u ( x, t ) tg(1 ) x sin( 1 ) = =, 1 + tg (1 ) 2 u ( x, t ) 1+ x откуда, согласно условию (1.14), имеем u( x + dx, t ) u( x, t ) sin( 2 ) sin( 1 ),.

x x Теперь, применяя для входящих в формулу (1.17) интегралов теорему о среднем, а для u ( x + dx, t ) – формулу Тейлора первого порядка с остатком в x форме Пеано, получаем 2 u( x, t ) 2 u (1, t ) dx = T0 dx + o(dx ) + F ( 2, t )dx, (1 ) x 2 t 1 и 2 принадлежат отрезку [x, x + dx ]. Почленно деля последнее где равенство на dx и осуществляя предельный переход при dx 0, получаем уравнение колебания струны следующего вида:

2u (x, t ) 2u (x, t ) ( x) = T0 + F ( x, t). (1.18) t 2 x Если струна дополнительно по всей длине связана с вязкоупругим основанием, то для описания ее колебания можно получить уравнение u( x, t ) u (x, t ) 2 2 t ( x) ( x, t ) u( x, t ) Q (x, t, )u (x, )d = T t x 2 (1.19) u( x, t ) (x, t ) + F ( x, t ), t где ( x, t ), ( x, t ) – коэффициенты жесткос ти и демпфирования основания;

Q ( x, t, ) – ядро релаксации, учитывающее изменение с течением времени физико-механических свойств материала основания (т. е. его старение).

Заметим, что при выводе уравнения (1.19) предполагалось, что реакция основания пропорциональна его деформации (модель Винклера).

В статистических задачах профиль струны u = u (x) определяется согласно (1.12), решением уравнения (x ) F (x ) u u=. (1.20) T0 T 1.6. Вывод уравнения продольных и крутильных колебаний стержня Для вязкоупругого тела при одномерном растяжении (сжатии) связь между деформацией (относительным удлинением) ( x, t ) и напряжением ( x, t ) представляется формулой ( x, t ) t ( x, t ) = E ( x, t ) ( x, t ) R( x, t, ) (x, )d + (x, t ), (1.21) t где E – модуль упругости;

R – ядро релаксации, учитывающее старение материала тела;

– коэффициент внутреннего трения. Заметим, если R 0 и 0, то получаем закон Гука для упругого тела.

Рассмотрим элемент с тержня (рис. 1.2), заключенный между поперечным сечением с координатами x и x + dx.

Рис. 1.2 Иллюстрация к выводу уравнения продольных колебаний стержня В сечении « x » на элемент действует сила N ( x, t ) = (x, t )S (x ), где S (x ) – площадь сечения, в сечении « x + dx » – сила N ( x + dx, t ) = ( x + dx, t )S ( x + dx )).

Предполагая, что на стержень действует внешняя нагрузка, распределенная по длине стержня с объемной плоскостью F ( x, t ), аналогично выводу уравнения (1.18) получаем уравнение продольных колебаний стержня следующего вида:

2 u( x, t ) = (S ( x ) ( x, t )) + S ( x) F (x, t ), ( x) S ( x) (1.22) x t где (x) – объемная плотнос ть материала с тержня;

u ( x, t) – продольное смещение сечения стержня с координатой x в момент времени t от положения, которое занимало это сечение, когда стержень находился в ненапряженном состоянии.

Учитывая, что u (x + dx, t) u ( x, t) u ( x, t) ( x, t ) = lim =, x dx 0 dx и подставляя (1.21) в (1.22), имеем 2 u( x, t) u (x, ) u ( x, t) t = S ( x )E (x, t ) ( x )S ( x) R (x, t, ) d + x x x t 2 u (x, t ) + S (x )F ( x, t ).

+ ( x, t )S (x ) xt Если боковая поверхность стержня скреплена с вязкоупругим основанием (модель Винкера), то приходим к следующему уравнению:

u( x, t) u (x, ) u ( x, t) t ( x )S ( x) R (x, t, ) d + = S ( x )E (x, t ) x x x t 2 2 u (x, t ) t ( x, t ) u( x, t ) Q (x, t, )u (x, )d + ( x, t )S (x ) (1.23) xt u ( x, t ) ( x, t) + S ( x )F ( x, t), t где (x, t ), ( x, t) – коэффициенты жесткос ти и демпфирования основания;

Q ( x, t, ) – ядро релаксации основания. Заметим, что форма записи уравнения (1.23) не изменится, если считать S и зависящими от времени t.

Статис тические продольные смещения сечений стержня u (x) определяются, согласно (1.23), решением уравнения [S (x )E ( x)u] ( x )u = S ( x) F (x ). (1.24) Для вязкоупругого стержня, находящегося в состоянии кручения (рис. 1.3), связь между напряжением, вызванным сдвигом образующей на угол, и этим углом может быть предс тавлена формулой ( x, t) t (x, t ) = G (x, t ) (x, t ) R ( x, t, ) (x, )d + ( x, t ), (1.25) t где G – модуль сдвига;

R – ядро релаксации стержня;

– коэффициент внутреннего трения.

Заметим, если R 0, 0, то получаем известный закон сдвига для упругого тела.

Если обозначить через u ( x, t) угол поворота сечения с координатой x, то (см. рис. 1.3) из равенства rdu = dx, имеем u =r. (1.26) x Рис.1.3. Иллюстрация к выводу уравнения крутильных колебаний стержня.

Крутящий момент M ( x, t), действующий в сечении стержня, S соответс твующем координате x, определяется формулой M ( x, t) = r ds.

S Отсюда, используя выражения (1.25), (1.26), получаем u (x, ) u ( x, t) t 2 u ( x, t ) R (x, t, ) d + ( x, t )J 0 ( x) M ( x, t) = J 0 ( x)G (x, t ), (1.27) x x xt где J 0 = r 2 dS – полярный момент инерции сечения.

S Рассмотрим элемент стержня, заключенный между поперечными сечениями с координатами x и x + dx (рис. 1.3). В сечении « x » действует крутящий момент M ( x, t), в сечении « x + dx » – M ( x + dx, t ). Предполагая, что на с тержень действует крутящий момент внешних сил, распределенный по длине стержня с линейной плотностью F ( x, t ), из уравнения динамического равновесия получаем 2 u(1, t ) M ( x, t ) (1 ) J 0 (1 ) dx + o(dx) + F (, t )dx, = x t где – плотность с тержня;

1 и 2 – принадлежат [x, x + dx ]. Откуда аналогично уравнению (1.18) получаем уравнение крутильных колебаний стержня 2u (x, t ) M ( x, t ) ( x) J 0 ( x ) = + F ( x, t ), x t которое, с учетом (1.27), принимает вид u ( x, ) 2 u ( x, t ) u( x, t ) t ( x )J 0 ( x) R ( x, t, ) d + = ( J 0 ( x )G ( x, t ) x x x t 2 u (x, t ) + ( x, t) J 0 (x ) ) + F ( x, t ).

xt Если боковая поверхность стержня скреплена с вязкоупругим основанием (модель Винклера), то для описания крутильных колебаний приходим к уравнению u ( x, t ) u (x, ) u ( x, t) t = J 0 ( x)G (x, t ) ( x )J 0 ( x) R (x, t, ) d + x x x t 2 u (x, t ) t ( x, t ) u( x, t ) Q (x, t, )u (x, )d ( x, t) J 0 (x ) + (1.28) xt u ( x, t ) ( x, t) + F ( x, t ), t где,, Q – коэффициенты жес ткости, демпфирования и ядро релаксации основания.

Заметим, что форма записи уравнения (1.27) не изменится, если считать и J 0 функциями двух переменных x и t.

Статические углы поворота u (x) сечений стержня при кручении определяются, согласно (1.28), решением уравнения [J ( x)G (x )u ] ( x)u( x ) = F (x ) (1.29) 1.7. Постановка статических краевых задач для струны и стержня В с татическом варианте профиль струны, продольные и угловые перемещения сечений стержня, согласно (1.20), (1.24) и (1.29), определяются решением уравнения L( y) (K ( x) y ) ( x ) y = g ( x), (1.30) где y ( x) = u ( x);

a x b ;

K ( x) = T0, g (x ) = F ( x), если рассматривается задача (1.20);

K ( x) = S ( x )E (x ), g ( x) = S (x )F ( x ), если – задача (1.24);

K ( x) = J 0 (x )G( x), g ( x) = F ( x), если – задача (1.29).

Перечислим основные типы граничных условий при x = a для уравнений (1.20), (1.24), (1.29) в обозначениях уравнения (1.30).

а) y (a) = 0 ;

это условие соответствует жесткому закреплению левого конца струны и стержня.

б) K (a ) y (a ) = qa ;

это условие соответствует заданию на левом конце стержня продольной силы N (a ) = q a для задачи (1.24) и заданию крутящего момента M (a ) = qa в случае задачи (1.29). В частности, если левый конец свободен, то qa = 0.

в) K (a ) y (a ) = a y (a ) ;

это условие соответс твует упругому закреплению левого конца стержня, когда qa = a y (a) ( qa или равно N (a ), или – M (a ) ), где a – соответс твующий задаче (1.24) или (1.29) коэффициент закрепления.

Аналогичные краевые условия могут быть заданы и на правом конце струны или стержня при x = b. Очевидно, что все возможные варианты краевых условий для уравнения (1.30) можно получить из условий (1.8) при соответс твующем выборе значений коэффициентов ai, bi.

Таким образом, рассматриваемые статические краевые задачи для струны и стержня математически формулируется так же, как и задача стационарной теплопроводности из раздела 1.2.

1.8. Краевые задачи в теории колебаний струн и стержней Предположим, что геометрические и прочностные характерис тики упругих тел (струны, стержня) и оснований, на которые они опираются, зависят только от x, и запишем уравнения движения без учета демпфирования и старения.

Уравнение продольных колебаний струны 2 u( x, t ) 2 u( x, t ) ( x) ( x)u (x, t ) + F ( x, t).

= T0 (1.31) t 2 x Уравнение продольных колебаний стержня 2 u ( x, t ) u ( x, t ) ( x )S ( x) ( x )u ( x, t ) + S ( x) F ( x, t ). (1.32) = S (x )E (x ) x x t Уравнение крутильных колебаний стержня 2 u( x, t ) u ( x, t ) ( x )J 0 (x ) (x )u ( x, t ) + F ( x, t ). (1.33) = J 0 ( x )G ( x) x x t Уравнения (1.31)–(1.33) являются уравнениями гиперболического типа.

Рассмотрим гармонические колебания упругих тел. В этом случае решение уравнений (1.31)–(1.33) и приложенную внешнюю нагрузку F ( x, t ) предс тавим в виде:

u ( x, t) = u* ( x) sin(t + ), F ( x, t ) = F* (x ) sin(t + ), (1.34) где (частота колебаний) и – постоянные. Тогда для u * (x ) = y (x ) получим уравнение (1.30), в котором F (x ) следует заменить на F* (x ), а (x) – на * (x ), где * ( x) = ( x) (x ) 2 соответствует уравнению (1.31), * ( x) = ( x) (x )S ( x ) 2 – уравнению (1.32), * ( x) = ( x) (x ) J 0 (x ) 2 – уравнению (1.33).

Приведем основные типы граничных условий при x = a.

а) u ( x, t ) = a (t );

это условие соответствует движению левого конца струны или стержня по закону a (t ).

u(a, t ) = q a (t );

это условие соответствует заданию на левом конце б) K (a ) x стержня продольной силы N (a, t ) = q a (t ) для задачи (1.32) и заданию крутящего момента M (a, t ) = qa (t ) в случае задачи (1.33). В частнос ти, если левый конец свободен, то qa = 0.

u(a, t ) = a [u(a, t) a (t )];

это условие соответс твует упругому в) K (a ) x закреплению левого сечения стержня, движущегося (вращающегося) по закону a (t ).

Предполагая функции a (t ), a (t), q a (t ) периодическими во времени, аналогично (1.34) положим a (t ) = a sin(t + ), a (t ) = a sin(t + ), qa (t ) = q 0 sin(t + ), 0 a где a, a, qa – постоянные. Тогда для u * (x ) = y (x ) будем иметь граничные 0 0 условия следующего вида:

а) y (a) = a ;

б) K (a ) y (a ) = q a ;

[ ] в) K (a ) y (a ) = a y (a ) a.

Условия на правом конце x = b задаются аналогично.

Замечание. Аналогичные краевые задачи получим в случае, когда u ( x, t ) = u 0 ( x) + u n ( x ) sin( n t + n ), n = F ( x, t ) = F0 (x ) + Fn ( x ) sin( n t + n ), n = a (t ) = n ( x) sin( n t + n ), n= a (t ) = n ( x) sin( n t + n ), n = qa (t ) = q0 + qn ( x) sin(n t + n ), n = где n, n, n, n, q n, q0 – постоянные;

u 0 ( x ) – решение с тационарных краевых задач, описанных в (1.5);

а u n (x ) – решение краевых задач, рассмотренных выше в этом параграфе.

Следует иметь в виду, что частота колебаний n являются в общем случае неизвестными величинами, определяемыми в процессе решения задачи.

2. Решение краевой задачи для линейного обыкновенного дифференциального уравнения второго порядка методом Галеркина.

Постановка задачи и алгоритм метода 2.1.

Рассмотрим следующую краевую задачу: требуется на отрезке [a, b] найти решение Y (x ) дифференциального уравнения L[y ] = y + p( x) y + q( x) y = f ( x), (2.1) удовлетворяющее условиям a0 y (a) + a1 y (a) = a2, (2.2) b0 y (b) + b1 y (b) = b2, где p( x), q( x), f (x ) – заданные функции, непрерывные на [a, b];

a0, a1, a 2, b0, b1, b2 – заданные действительные числа, причем a0 + a12 0, b0 + b12 0.

2 Напомним, что в отличие от имеющей всегда единс твенное решение задачи Коши для уравнения (2.1), краевая задача (2.1), (2.2) может иметь или одно решение, или бесконечно много решений, или, наконец, может совсем не иметь решений.

Везде далее будем предполагать существование единс твенного решения Y (x ) поставленной краевой задачи, что часто вытекает из физического смысла того явления или процесса, математическое моделирование которого привело к задаче (2.1), (2.2).

В методе Галеркина для нахождения приближенного решения рассматриваемой задачи строится функциональная последовательнос ть {y n (x)} из пробных решений y n (x ) следующим образом.

Задаемся на отрезке [a, b] некоторой системой дважды непрерывно дифференцируемых функций u 0 ( x), u1 ( x),..., un (x ) таких, что u0( x) удовлетворяет краевым условиям (2.2), а функции u1 (x ), u 2 ( x),..., un ( x), называемые пробными функциями, линейно независимы на [a, b] и удовлетворяют однородным краевым условиям a0u (a ) + a1u (a) = 0, (2.3) b0u (b) + b1u (b) = 0.

Составляем функцию n y n (x ) = u0 ( x) + C iui ( x) (2.4) i = с неизвестными пока постоянными коэффициентами C1, C 2,..., C n. Подчеркнем, что в силу линейнос ти условий (2.2), функция (2.4) при любых значениях C1,..., C n удовлетворяет этим условиям. Подставляя функцию y n (x ) из (2.4) вместо y (x) в уравнение (2.1), получаем функцию R (C1, C 2,..., Cn, x ) = L[u0 ] f ( x ) + Ci L[ui ], n (2.5) i = которая называется невязкой. Как видно из (2.5), невязка линейно зависит от параметров C1, C 2,..., C n и является характеристикой уклонения функции (2.4) от точного решения Y (x ) задачи (2.1), (2.2). Во всяком случае, если при некоторых значениях параметров C1, C 2,..., C n невязка на [a, b ] тождественно равна нулю, то Y (x ) yn (x ) в силу единственности Y (x ).

Однако в общем случае невязка оказывается отличной от нуля. Поэтому подбираем значения параметров C1,..., C n так, чтобы невязка в каком-то смысле была бы наименьшей. В обобщенном методе Галеркина значения параметров C1,..., C n определяются из системы уравнений (R(C1,...,C n, x ),Wk (x )) = 0, k = 1, n, (2.6) где b ( ( x ), g ( x)) = ( x )g ( x )dx, (2.7) a а W1 ( x),...,Wn (x ) – заданные непрерывные и линейно независимые на [a, b] функции, часто называемые поверочными функциями. Заметим, что если в качестве поверочных функций взять пробные, то получится метод Галеркина в авторском варианте [1]. Заметим также, что если W1 ( x),...,Wn (x ) входят в полную сис тему функций, то при n равенства (2.6) свидетельс твуют об ортогональнос ти невязки всем элементам полной системы [3]. Значит, невязка сходится при n к нулю в среднем, и можно ожидать сходимости последовательности (2.4) к точному решению Y (x ) в среднем, т. е.

(Y ( x) y n ( x), Y (x ) yn (x )) = 0.

lim n Записав условие (2.6) в развернутом виде, для определения значений параметров получаем неоднородную систему линейных C1,..., C n алгебраических уравнений n-го порядка n akj C j = bk, k = 1.n, (2.8) j = где [] b akj = ( L u j, Wk ( x)) = (u j + puj + qu j )Wk dx, a (2.9) b bk = ( f ( x) L[u 0 ], Wk (x )) = ( f ( x) u pu qu0 )Wk dx.

0 a Решив систему (2.8) и подставив определяемые этим решением значения параметров C1,..., C n в (2.4), заканчиваем построение пробного решения y n (x ).

Опишем теперь возможный алгоритм приближенного решения задачи (2.1), (2.2) методом Галеркина, предполагая, что y n (x ) сходится к Y (x ) при n.

1. Подготовленный шаг алгоритма. На этом шаге выбираем функцию u 0 ( x ), пробные функции u1 (x ),..., u n ( x ) и поверочные функции W1 ( x),...,Wn (x ).

Находим функцию R 0 ( x ) = L[u0 ] f (x ), т. е. невязку от подстановки u 0 ( x ) в уравнение (2.1). Если x [a, b ] : R0 ( x) = 0, то u 0 ( x ) = Y ( x), и вычисления заканчиваем. Если же R 0 ( x) 0, то переходим к следующему шагу алгоритма.

2. Первый шаг алгоритма. Строим y1 (x ) = u0 ( x) + C1u1 ( x), определив значение C1 из решения системы (2.8) при n = 1. Находим невязку R ( C 1, x ) = L [u 0 ] f ( x ) + C 1 L [u 1 ] = R 0 ( x ) + C 1 L [u 1 ].

Если x [a, b ] : R (C1, x ) = 0, то Y (x ) = y1 ( x), и задача решена, если же R (C, x ) 0, то находим max y1 ( x) u 0 (x ) = [a, b] или max R (C1, x ) = 12.

[a, b] Если 11 1 или 12 2, где 1 и 2 заданные меры точнос ти приближенного решения, то полагаем Y (x ) y1 ( x) и вычисления заканчиваем, если же 11 или 12 2, то переходим к вычислениям на следующем шаге и т. д.

Таким образом, на m -м (m 1) шаге алгоритма строим функцию m y m ( x) = u 0 ( x) + C iu i ( x ), i = определив значения C1,..., C m из решения системы (2.8) при n = m, и определяем невязку m R (C1,..., Cm, x) = R0 ( x) + C i L[u i ].

i = Если x [a, b ] : R (С1,..., C m, x ) = 0, то Y (x ) = ym ( x ), и вычисления заканчиваем.

Если R (C1,..., C m, x) 0, то находим m1 = max ym ( x ) y m 1 (x ) или m 2 = max R (C1,..., C m, x.

[a, b] [a, b ] Если m1 1 или m 2 2, то Y (x ) ym ( x ), если же m1 1 или m 2 – переходим к (m + 1) -му шагу.

2.2.Построение систем пробных и поверочных функций Известно, что степенные функции 1, x, x 2,..., x n,... линейно независимы на всей числовой прямой R и, следовательно, на любом ее отрезке [a, b] R.

Покажем, что на любом отрезке [a, b] линейно независима любая система многочленов последовательных с тепеней. Рассмотрим произвольную систему многочленов:

P0 ( x) = A00 0;

P1 (x ) = A11 x + A10, A11 0;

P2 ( x) = A22 x 2 + A21x + A20, A22 0;

LLLLLLLLLLLLLLLLLL Pn ( x ) = Ann x n + Ann1 x n 1 +... + An 0, Ann и решим относительно неизвестных 0, 1,..., n определенное на R тождество 0 P0 ( x) + 1P1 ( x) +... + n Pn (x ) 0. (2.10) Из условий тождественного равенства нулю многочлена n -й степени (равенство нулю коэффициентов при всех степенях x ) последовательно получаем n Ann = 0 n = 0;

n Ann1 + n1 An 1 n 1 = 0 n 1 = 0;

................................

n An 0 + n 1 An 10 +... + 1 A10 + 0 A00 = 0 0 = 0.

Таким образом, условие (2.10) выполняется тогда и только тогда, когда 0 = 1 =... = n = 0, т. е. система многочленов P0 ( x ),..., Pn ( x) и любая подсистема из них линейно независима на R и, следовательно, на любом [a, b] R.

Для построения u 0 ( x ) и линейно независимой на [a, b], системы пробных функций u1 (x ),..., un (x ), являющихся многочленами, можно применить метод неопределенных коэффициентов.

Например, предположим u 0 = A = P0 ( x), из условий (2.2) получаем систему линейных алгебраических уравнений относительно A a0 A = a2, b0 A = b2.

В том случае, когда эта система совместна, коэффициент A определяется. Если система не совместна, то ищем аналогичным образом u 0 ( x ) в виде u 0 ( x ) = A + Bx = P1 (x ) и т. д., до тех пор, пока не будет найдена u 0 ( x) = Pr0 (x ), удовлетворяющая условиям (2.2).

Далее, используя условия (2.3), методом неопределенных коэффициентов определяем последовательно так же, как и u 0 ( x ), u1 (x ) = Pr1 ( x), r1 0;

u 2 ( x) = Pr2 ( x), r2 r1 ;

LLLLLLLLLL u n ( x ) = Prn ( x), rn rn 1.

Пример 1. Построить u 0 ( x ) и систему из пяти пробных функций для задачи с краевыми условиями y(0) + y (0) = 1, (2.11) y (1) + y (1) = 4.

Решение. Пусть u0 (x) = A, тогда u = 0 и условия (2.11) дают несовместную систему из уравнений A = 1 и A = 4.

Пусть u 0 = A + Bx, тогда u 0 = B и условия (2.11) дают A + B = 1, A + B = 0, A = 6, A + 2 B = 4, B = 5, B = 5.

Итак, u 0 = 6 5 x.

Определяем u1 (x ). Если u1 = A или u1 = A + Bx, то однородные условия, соответс твующие условиям (2.11), выполняются, если u1 0, что невозможно из-за требования линейной независимости пробных функций.

Ищем u1 (x ) = A + Bx + Cx 2 (C 0), тогда u1 = B + 2Cx, и из однородных условий, соответствующих (2.11), получаем систему A + B = 0, A + 2 B + 3C = 0.

Решая ее методом Гаусса, имеем A + B = 0, B + 3C = 0.

Видим, что система имеет множество решений G = {(A, B, C ): A = 3, B = 3, C =, R }.

Выбираем одно решение из G при =, тогда u1 ( x ) = 1 x + x 2.

Аналогично, используя формулу u k = A0 + A1x +... + Ak +1 x k + 1, находим 1 1 1 u 2 ( x) = 1 x + x 3, u3 ( x) = 1 x + x 4, u 4 ( x) = 1 x + x 5, u5 ( x) = 1 x + x 6.

4 5 6 Пример 2. Построить u 0 ( x ) и систему из трех пробных функций для задачи с краевыми условиями y (0) + y (0) = 1, (2.12) y (2) y (2) = 2.

Решение. Если u 0 ( x ), то условия (2.12) приводят к несовместной системе A = 1, A = 2.

Предположим, что u 0 = A + Bx, тогда u = B и условия (2.12) дают A + B = 1, A + B = 1, A + B = 1, A + 2 B B = 2, A + B = 2, 0 = 1, тоже несовместную систему.

Полагаем u 0 = A + Bx + Cx 2, тогда u = B + 2Cx и условия (2.12) дают A + B = 1, A + B = 1, A + B = 1, A + 2 B B = 2, A + B = 2, 0 = 1, которая несовместна.

Ищем u 0 ( x ) в виде u 0 = A + Bx + Cx 2 + Dx 3, тогда u = B + 2Cx + 3Dx 2, и из (2.12) имеем A + B = 1, A + B = 1, A + 2 B + 4C + 8D B 4C 12D = 2, A + B 4D = 2.

Решаем полученную систему методом Гаусса в матричной форме, чтобы найти все решения системы.

Прямой ход метода:

ABCD ABCD ADCB 1 1 0 0 1 1 1 0 0 1 1 0 0 1 1 1 0 4 2 ~ 0 0 0 4 1 ~ 0 4 0 0 1.

Видим, что система совместна, ибо ранг матрицы системы (rg ) равен рангу расширенной матрицы и равен 2. Так как число неизвестных системы 4 больше rg = 2, то система неопределена, и все множество решений G 0 системы получаем обратным ходом метода Гаусса, придавая двум неизвестным C и B произвольные значения. Получаем G0 = ( A, B, C, D ) : A = 1 1, B = 1, C = 2, D = ;

1, 2 R.

Выбираем решение из G 0 при 1 = 2 = 0. Тогда u0 ( x) = 1 x3.

Определяем теперь u1 (x ). Если u1 (x ) = A 0, то однородные условия, соответс твующие условиям (2.12), выполняются при A = 0, что недопустимо.

Пусть u1 (x ) = A + Bx, u1 (x ) = B, и из однородных условий, соответс твующих условиям (2.12), имеем A + B = 0, A+ B =0.

A + B = 0, Эта система неопределена, ее множество решений G1 = {( A, B ) : A =, B = ;

R}.

Выбираем одно ненулевое решение при = 1, тогда u1 (x ) =1 x.

Ищем u 2 (x ). Пусть u 2 (x ) = A + Bx + Cx 2, (C 0 ), тогда u 2 = D + 2Cx и однородные условия дают систему A + B = 0, A + B + 0C = 0.

Решая ее методом Гаусса, находим множество решений G 2 = {( A, B, C ) : A = 1, B = 1, C = 2 ;

1, 2 R}.

Выбирая одно ненулевое решение ( C 0 ), при 1 = 2 = 1, получаем u 2 (x ) = 1 x x 2.

u 3 ( x) = A + Bx + Cx 2 + Dx 3, ( D 0 ), Находим Если то u 3 ( x).

u 3 ( x) = B + 2Cx + 3Dx 2, и из однородных условий имеем систему A + B = 0, A + B = 0, A + B + 0C 4D = 0, D = 0, которая противоречит условию D 0.

u 3 ( x) = A + Bx + Cx 2 + Dx 3 + Ex 4, E 0, Пусть теперь тогда u 3 ( x) = B + 2Cx + 3Dx 2 + 4 Ex 3, и из однородных условий получаем систему A + B = 0, A + 2 B + 4C + 8D + 16 E B 4C 12 D 32 E = 0, A + B = 0, A + B + 0C 4D 16 E = 0.

Решая ее методом Гаусса, получаем множество решений G3 = {( A, B, C, D, E ) : A = 1, B = 1, C = 2, D = 4 3, E = 3 ;

1, 2, 3 R}.

Выбирая одно ненулевое решение ( E 0 ) при 1 = 1, 2 = 1, 3 = 1, имеем u 3 ( x) = 1 x + x 2 4 x 3 + x 4.

Применяя метод неопределенных коэффициентов, можно строить системы функций, используя другие системы линейно независимых на R функций, такие как e1 x, e 2 x,..., en x,... ;

ex, xex,..., x n ex,... ;

1, cos( x), sin( x ), cos(2 x ), sin(2 x),... и т. п.

Важным ис точником для пос троения ортогональных на [a, b] пробных функций является множество решений задачи, называемой задачей на собственные значения для дифференциального оператора L[ y ] = y.

Рассмотрим конкретный пример такой задачи.

Пример 3. Требуется найти действительные значения параметра, при которых существуют нетривиальные решения дифференциального уравнения y + y = 0, (2.13) удовлетворяющие однородным условиям y (0) + y (0) = 0, (2.14) y (1) + y (1) = 0.

Решение. Пус ть = 0, тогда общее решение уравнения (2.13) будет иметь вид y = C1 x + C 2. Пытаясь удовлетворить условиям (2.14), получаем C1 + C 2 = 0, C1 + C2 = 0, C1 = 0, 2C1 + C 2 = 0, C 2 = 0, C 2 = 0.

Таким образом, = 0 не является собственным значением, так как ему соответс твует единс твенное тривиальное ( y 0 ) решение задачи (2.13), (2.14).

Пусть 0, тогда y = c1e x + c 2 e x, = | |, и условия (2.14) приводят к системе уравнений C1 + C 2 C1 + C 2 = 0, (1 )C1 + (1 + )C 2 = 0, C1 e + C2 e C1 e + C 2 e = 0, (1 )C1e + (1 + )C2 e = 0, (1 )C1 + (1 + )C 2 = 0, C2 = 0, (1 + )(e + e )C 2 = 0, C1 = 0.

Следовательно, среди отрицательных действительных чисел собственных значений задачи (2.13), (2.14) нет.

Пусть теперь 0. Тогда y = C1 cos(x) + C 2 sin(x), =, и краевые условия (2.14) дают C1 + C 2 = 0, C1 cos + C 2 sin C1 sin + C 2 cos = 0, C1 = C 2, C1 (cos sin ) + C 2 (sin + cos ) = 0, C1 = C 2, C1 = C 2, C 2 ( (cos sin ) + sin cos ) = 0, C2 (1 + ) sin = 0.

Видим, что существуют нетривиальные решения задачи (2.13), (2.14), если sin = 0, т. е. = n, n = 1, 2,...

Таким образом, множество собственных значений определяется формулой n = (n ) 2, n = 1, 2,..., а множество собственных функций, соответс твующих собственному значению n, имеет базисную функцию y n = sin( n x) n cos( n x ) = sin( nx ) n cos( nx ).

Для того чтобы убедиться в ортогональнос ти на [0, 1] функций y n ( x), y m ( x ) (n m), достаточно проверить, что ( y n, y m ) = y n ( x) y m ( x )dx = 0.

При выборе систем поверочных функций полезно вспомнить и о других системах функций, ортогональных на некотором отрезке. Например, известно [3], что многочлены Лежандра, определяемые формулой 1 dn Pn (t ) = n (t 1) n, n = 0, 1, 2,... (2.15) n 2 n! dt ортогональны на [–1,1]. Так что, если в качестве поверочных функций Wk (x) решено взять, например, первые пять многочленов Лежандра, ортогональных на [a, b], то в первые пять выражений из (2.15):

W1 = P0 (t ) = 1, W2 = P1 (t ) = t, W3 = P2 (t ) = (3t 2 1), 1 W4 = P3 (t ) = (5t 3 3t), W5 = P4 (t) = (35t 4 30t 2 + 3).

2 следует подставить a +b t= x.

ba 2.3. Задание к лабораторной работе Методом Галеркина найти наиболее точное приближенное решение краевой задачи d +d x d y + 0 2 1 y + y = 0, x 1 1 x a0 y (a) + a1 y (a ) = a2, (2.16) b0 y (b ) + b1 y (b) = b2, построенное при помощи системы из n пробных функций – многочленов и двух систем поверочных функций, одна из которых составлена из пробных функций, а вторая – из многочленов Лежандра. За меру точности выбрать (по указанию преподавателя) или max y m (x ) y m 1 ( x ), [ a,b ] или max R (C1,..., C m, x ), [ a,b ] Варианты заданий, определяемые различными наборами заданий параметров d 0, d1, d 2, a 0, a1, a 2, b0, b1, b 2, a, b задачи (2.16) приведены в таблице 2.1.

Лабораторная работа выполняется с использованием прикладной системы MathCAD, в которой реализуется алгоритм построения пробных решений y m (x) методом Галеркина.

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

Числовые данные:

a, b – концы отрезка интегрирования;

n – максимальное число параметров C i в пробном решении.

Значение параметра n задает преподаватель.

Строчные данные:

Аналитические выражения для функций u 0 ( x ),..., u n ( x) и для функций W1 ( x),..., Wn ( x ) ;

Аналитические выражения набираются клавишами дисплея по определенным правилам, которые программа при необходимости может напомнить.

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

Анализируя данные этих таблиц, необходимо найти обоснованный ответ на поставленную в лабораторной работе задачу.

Таблица 2. Варианты заданий лабораторной работы № a0 a1 a2 b0 b1 b2 d1 d d a b вар.

1 0 0,8 1 0 – 0,5 1 0 0,5 2 0 2 0 0,8 1 0 0 1 0 0,1 2 0 3 0 0,8 0 1 0,5 1 0 0,2 2 0 4 0 0,8 0 1 0 1 0 0,1 2 0 5 0 0,8 1 0 – 0,4 0 1 – 0,2 2 0 6 0 0,8 1 0 – 0,5 1 0 0,5 0 2 7 0 0,8 1 0 0 1 0 0,1 0 2 8 0 0,8 0 1 0,5 1 0 0,2 0 2 9 0 0,8 0 1 0 0 1 0,1 0 2 10 0 0,8 1 0 – 0,4 0 1 – 0,2 0 2 11 0 0,6 1 0 0,2 1 0 0,8 2 0 12 0 0,6 0 1 0,15 0 1 0,2 2 0 13 0 0,6 0 1 – 0,1 1 0 0,4 2 0 Продолжение таблицы 2. № a0 a1 a2 b0 b1 b2 d1 d d a b вар.

14 0 0,6 1 0 – 0,2 1 0 – 0,8 2 0 15 0 0,6 1 1 0,5 0 1 – 1,0 0 2 16 0 0,6 1 1 0,4 1 0 1 0 2 17 0 0,6 1 0 0,4 1 1 0,2 0 2 18 0 0,4 1 0 0 1 0 0,9 0 2 19 0 0,4 1 0 1 1 0 – 0,13 0 2 20 0 0,4 1 0 0 1 0 – 1,1 0 2 21 0 0,4 1 0 0 1 0 – 0,35 0 2 22 0 0,4 0 1 1 0 1 –1 2 0 23 0 0,5 1 –1 0,2 0 1 1 2 0 24 0 0,5 –1 0 0,4 1 0 0,6 2 0 25 0 0,5 0 1 2 1 0 0,4 2 0 2.4. Выполнение работы в компьютерном классе 1. Прежде чем начать выполнение лабораторной работы на ЭВМ, внимательно ознакомьтесь с данной инструкцией.

2. При необходимости включите сами (или попросите лаборанта) питание компьютера. После того, как система загрузится, запускаем двойным щелчком левой кнопки мыши на рабочем столе программу Mathcad, если же ярлык отсутс твует, тогда открываем программу через кнопку «Пуск» (Программы Mathsoft Mathcad).


3. Узнайте у лаборанта расположение файла ODU.mcd и откройте его (File Open или, если программа русифицирована, Файл Открыть). При любой ошибке ввода программы нужно обратиться к лаборанту.

4. Прочитайте в начале файла задание на лабораторную работу и просмотрите пример выполнения работы, для которого исследование уже проведено. Программа файла ODU.mcd состоит из шести пунктов «1.

Постановка задачи», «2. Получение точного решения», «3. Получение приближенного решения методом Галеркина», «4. Получение приближенного решения вариационным методом Ритца», «5. Получение приближенного решения интегральным методом наименьших квадратов», «6. Выводы». Для выполнения лабораторной работы «Решение краевой задачи для линейного обыкновенного дифференциального уравнения второго порядка методом Галеркина» необходимо использовать пункты 1, 2 и 3. Цели и задачи каждого из пунктов описаны ниже.

5. Для набора функций нужно либо воспользоваться всплывающим меню инструментов «Calculator», либо ввести ее с клавиатуры, используя следующие символы арифметических действий и стандартных функций: сложение – ‘+’;

вычитание – ‘–‘;

умножение – ‘*’;

деление – ‘/’;

возведение в степень – ‘^’;

квадратный корень – ‘\’;

синус – sin(x);

косинус – cos(x);

экспонента – exp(x);

натуральный логарифм – ln(x). При вводе числовых данных, являющихся десятичными дробями, целую и дробную части нужно разделять точкой (например, 0.5, 1.5 и т. д.).

6. Порядок выполнения работы Вам укажет программа подсказками и заданиями, выделенными красным цветом.

2.5. Порядок выполнения лабораторной работы Рекомендуется следующий порядок выполнения лабораторной работы.

1. Изучить разделы 1.1–1.2, 2.1–2.4 и подготовить ответы на контрольные вопросы из раздела 2.6.

2. Пройти собеседование с преподавателем, получить допуск к выполнению работы в диалоге с ПЭВМ, номер варианта задания и значение параметра n.

3. В соответс твии с вариантом задания выполнить подготовительный шаг алгоритма метода Галеркина и подготовить, если u 0 ( x) не является точным решением задачи, все числовые и строчные исходные данные для расчетов на ПЭВМ.

4. Выполнить основную расчетную часть лабораторной работы в системе MathCAD. Следует переписать с экрана дисплея значения коэффициентов C i пробных решений, составить итоговые таблицы пробных решений и их невязок.

5. Оформить и защитить отчет по лабораторной работе, который должен содержать титульный лист, математическую постановку задачи, результаты выполнения подготовительных расчетов, основные результаты расчетов на ПЭВМ и итоговый отчет.

2.6. Тестирующий пример Методом Галеркина найти на [0, 1] приближенное решение краевой задачи y 3 y + 2 y = 2 x 2 6 x + 2, y (0) + y (0 ) = 1, (2.17) y (1) + y (1) = 4.

1. Запускаем программу MathCAD. Открываем файл ODU.mcd (текст программы приведен в приложении A). В пункте «Пос тановка задачи» вводим числовые параметры a0, a1, a 2, b0, b1, b2, a, b и функции p( x), q (x ), f (x ), входящие в задачу (2.17) p( x) := 3, q( x) := 2, f (x ) := 2 x 2 6 x + 2, a := 0, b := 1, a0 := 1, a1 := 1, a 2 := 1, b0 := 1, b1 := 1, b2 := 4.

Замечание. Для задачи (2.16) необходимо еще ввести числовые параметры d 0, d 1, d 2, входящие в функции p( x), q (x ).

2. В пункте «Получение точного решения в системе MathCAD» записываем дифференциальное уравнение (2.17) в виде нормальной системы дифференциальных уравнений второго порядка y 0 = y1, y1 = 3 y1 2 y 0 + 2 x 2 6 x + 2.

Далее с помощью функции bvalfit (см. п. 8.2) краевую задачу приводим к задаче Коши, получая начальные условия y (0) = 0,153224, y (0 ) = 1 y (0 ) = 0,846776.

Решая полученную задачу Коши для нормальной системы второго порядка с помощью функции rkfixed (см. п. 8.2), находим решение дифференциального уравнения, разбив отрезок [0, 1] на N = 100 частей, т. е. с шагом интегрирования h = 0,01 (в дальнейшем будем называть его точным решением).

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

Таблица 2. Таблица точного решения задачи.

x 0 0,1 0,2 0,3 0,4 0, y 0,846776 0,86522 0,886703 0,905491 0,914371 0, x 0,6 0,7 0,8 0,9 1, y 0,863987 0,779347 0,632877 0,402849 0, График точного решения имеет вид Рис.2.1. График точного решения.

Заметим, что поставленная задача имеет единс твенное точное решение вида 1 e2 + 7 x 7 + e 2x y= e e + x 2 = 1,5403e x 0, 6935e 2 x + x 2, (2.18) 2 (e e ) 3(e e ) 2 которое получено аналитическим методом, известным из теории линейных дифференциальных уравнений с постоянными коэффициентами.

Если найти значения этого решения y (x) в тех же промежуточных точках, то получим значения таблицы 2.2, т. е. компьютерное решение найдено с точностью 10 6. Если сравнить точное аналитическое решение с решением, полученным в системе MathCAD, видим, что локальная погрешность не превышает 1,77 10 8. Следовательно, решение найдено достаточно точно.

Замечание. Компьютерное решение будет тем точнее, чем больше число точек разбиения введено в функцию rkfixed. Например, при N = 1000, локальная погрешность не превышает 3,45 10 12.

3. В пункте «Получение приближенного решения методом Галеркина»

вводим порядок пробного решения n := 5. В качестве пробных функций u 0 (x), u1 ( x),..., u5 (x ) используем функции, построенные в примере 1 раздела 2.2, 1 u 0 ( x) = 6 5x, u1 ( x ) = 1 x + x 2, u 2 ( x) = 1 x + x 3, 3 1 1 u 3 ( x ) = 1 x + x 4, u 4 ( x ) = 1 x + x 5, u5 ( x ) = 1 x + x 6.

5 6 Тогда f (x ) L[u0 ] = 2x 2 + 4 x 25, ( ) L[u1 ] = 2x 2 12 x + 17, ( ) L[u2 ] = 2 x 3 9 x 2 2 x + 20, ( ) L[u3 ] = 2 x 4 12 x 3 + 12 2 10 x + 25, ( ) L[u4 ] = 2 x 5 15 x 4 + 203 + 12x + 30, ( ) L[u5 ] = 2 x 6 18 x 5 + 30 4 14 x + 35.

3.1. Воспользовавшись указаниями из раздела 2.2, в качестве поверочных функций возьмем пробные u1 (x ),..., u5 ( x) 1 1 W1 (x ) = 1 x + x 2, W2 ( x) = 1 x + x 3, W3 ( x ) = 1 x + x 4, 3 4 1 W4 ( x) = 1 x + x 5, W5 ( x ) = 1 x + x 6.

6 В результате расчета по программе при n = 5 получим вектор коэффициентов С = (1,132936 2, 499320 2, 647392 0,073920 1, 213380), следовательно, пробное решение имеет вид y 5 (x ) = u0 (x ) + 1,132936u1 (x ) 2, 499320u2 ( x ) 2,647392u 3 ( x ) + + 0, 073920u 4 ( x) 1, 213380u5 ( x).

Основные результаты расчета при n 5 (т. е. подставляя последовательно n = 0,1, 2,3, 4,5 ) представлены в таблицах 2.3 и 2.4. В приложении A приведен пример при n = 5.

Таблица 2. Таблица значений пробных решений n=0 n=1 n=2 n=3 n=4 n= x 0,0 6,0 0,546243 0,878844 0,843793 0,846932 0, 0,1 5,5 0,573439 0,905606 0,860505 0,865612 0, 0,2 5,0 0,564277 0,951373 0,87879 0,887337 0, 0,3 4,5 0,518757 1,000717 0,895436 0,906106 0, 0,4 4,0 0,436879 1,038207 0,904176 0,914814 0, 0,5 3,5 0,318642 1,048413 0,895683 0,904636 0, 0,6 3,0 0,164046 1,015907 0,857573 0,864404 0, 0,7 2,5 – 0,026908 0,925257 0,774403 0,779985 0, 0,8 2,0 – 0,254220 0,761035 0,627672 0,633662 0, 0,9 1,5 – 0,517890 0,507811 0,395822 0,403509 0, 1,0 1,0 – 0,817919 0,150155 0,054236 0,062772 0, Таблица 2. Таблица невязок пробных решений n=0 n=1 n=2 n=3 n=4 n= x 0,0 25,00 – 5,904624 2,837772 – 0,583105 0,090419 – 0, 0,1 24,58 – 4,179480 1,126703 – 0,049648 – 0,014480 0, 0,2 24,12 – 2,567052 – 0,123462 0,191588 – 0,031421 0, 0,3 23,62 – 1,067341 – 0,943584 0,225923 – 0,010841 – 0, 0,4 23,08 0,319653 – 1,364521 0,132562 0,013263 – 0, 0,5 22,50 1,593931 – 1,417133 – 0,015406 0,022085 0, 0,6 21,88 2,755491 – 1,132279 – 0,151009 0,010772 0, 0,7 21,22 3,804335 – 0,540820 – 0,213391 – 0,012818 0, 0,8 20,52 4,740462 0,326387 – 0,147810 – 0,029362 – 0, 0,9 19,78 5,563873 1,438480 0,094356 – 0,009312 – 0, 1,0 19,00 6,274566 2,764601 0,555617 0,085862 0, Анализируя таблицы, видим, что наилучшее приближение к точному решению дает пробное решение y 5 (x ), для которого max Y ( x) y 5 (x ) 0,000041, [a,b ] max R (C1,..., C 5, x ) 0,010887, [a,b ] max y5 ( x) y 4 ( x) 0,000806.

[a,b ] 3.2. Воспользовавшись указаниями из раздела 2.2, в качестве поверочных функций возьмем многочлены Лежандра ( ) 3(2x 1)2 1, W1 ( x ) = 1, W2 (x ) = 2 x 1, W3 (x ) = ( ) ( ) 5(2 x 1)3 3(2 x 1), W5 ( x) = 35(2 x 1)4 30(2 x 1)2 + 3.

1 W4 ( x ) = 2 В результате расчета по программе при n = 5 получим вектор коэффициентов С = (1,135995 2,510843 2, 638116 0,080296 1, 220556), следовательно, пробное решение имеет вид y 5 ( x ) = u0 ( x ) + 1,135995u1 ( x ) 2,510843 u2 ( x ) 2,638116 u3 ( x ) + + 0,080296 u 4 ( x ) 1,220556 u5 ( x ).

Основные результаты расчета при n 5 представлены в таблицах 2.5 и 2.6.

В приложении A приведен пример при n = 5.

Таблица 2. Таблица значений пробных решений n=0 n=1 n=2 n=3 n=4 n= x 0,0 6,0 0,257143 0,902890 0,845255 0,846800 0, 0,1 5,5 0,312286 0,928405 0,861827 0,865451 0, 0,2 5,0 0,329143 0,974890 0,880057 0,887081 0, 0,3 4,5 0,307714 1,026393 0,896811 0,905721 0, 0,4 4,0 0,248000 1,066960 0,905854 0,914300 0, 0,5 3,5 0,150000 1,080636 0,897855 0,904025 0, 0,6 3,0 0,013714 1,051468 0,860381 0,863751 0, 0,7 2,5 – 0,160857 0,963503 0,777897 0,779356 0, 0,8 2,0 – 0,373714 0,800786 0,631773 0,633114 0, 0,9 1,5 – 0,624857 0,547364 0,400276 0,403070 0, 1,0 1,0 – 0,914286 0,187283 0,058574 0,062412 0, Таблица 2. Таблица невязок пробных решений n=0 n=1 n=2 n=3 n=4 n= x 0,0 25,00 – 7,542857 3,206936 – 0,578802 0,080526 – 0, 0,1 24,58 – 5,704000 1,374150 – 0,034754 – 0,019452 0, 0,2 24,12 – 3,981714 0,021919 0,210140 – 0,031797 0, 0,3 23,62 – 2,376000 – 0,881665 0,242606 – 0,007511 – 0, 0,4 23,08 – 0,886857 – 1,368509 0,143171 0,018995 – 0, 0,5 22,50 0,485714 – 1,470520 – 0,013838 0,028644 0, 0,6 21,88 1,741714 – 1,219607 – 0,160293 0,016448 0, 0,7 21,22 2,881143 – 0,647676 – 0,234266 – 0,009752 0, 0,8 20,52 3,904000 0,213364 – 0,180028 – 0,030533 – 0, 0,9 19,78 4,810286 1,331607 0,051950 – 0,016143 – 0, 1,0 19,00 5,600000 2,675145 0,504997 0,072243 0, Анализируя решения при n = 0,1,...,5, видим, что наилучшее приближение к точному решению дает пробное решение y 5 (x ), для которого max Y ( x) y 5 (x ) 0,000023, [a,b ] max R (C1,..., C 5, x ) 0,00879, [a,b ] max y5 ( x) y 4 ( x) 0,000387.


[a,b ] 4. Анализируя полученные погрешности, видим, что для задачи (2.17) многочлены Лежандра, как поверочные функции, дают лучшее приближение решения этой краевой задачи.

2.7. Вопросы для самоконтроля 1. Найдите решение краевой задачи (2.17) аналитическим методом.

2. Каковы отличия краевой задачи от задачи Коши?

3. Каким условиям должны удовлетворять пробные функции в методе Галеркина?

4. Как находится функция, названная в методе Галеркина невязкой пробного решения?

5. Какими свойствами должны обладать поверочные функции в методе Галеркина?

6. Как в методе Галеркина строится система линейных алгебраических уравнений для определения коэффициентов пробного решения? Проверьте истиннос ть формул (2.8), (2.9).

7. В каком случае невязка пробного решения сходится к нулю в среднем при n ?

8. Опишите алгоритм приближенного решения краевой задачи для линейного дифференциального уравнения второго порядка методом Галеркина.

9. Приведите пример пос троения пробных функций методом неопределенных коэффициентов.

10.Напишите два многочлена Лежандра, ортогональные [2,4], и проверьте их ортогональнос ть.

11.Напишите уравнение и краевые условия задачи на собственные значения.

12.Напишите две собственные функции задачи (2.13), (2.14) и проверьте их ортогональнос ть.

13.Приведите физические интерпретации изучаемой краевой задачи.

3. Решение краевой задачи для линейного обыкновенного дифференциального уравнения второго порядка вариационным методом Ритца 3.1. Постановка задачи и алгоритм метода Рассмотрим краевую задачу в следующей постановке.

Требуется на отрезке [a, b] найти решение Y (x ) дифференциального уравнения L[ y ] (K ( x ) y ) ( x ) y = g (x ), (3.1) удовлетворяющее двум краевым (или граничным) условиям a 0 y (a) + a1 y (a) = a 2, (3.2) b0 y( a) + b1 y (a ) = b2, где K (x), K (x), (x ), g (x ) – заданные непрерывные на [a, b] функции ( K ( x) 0 );

a 0, a1, a 2, b0, b1, b2 – заданные действительные числа, причем a 2 + a1 0, b0 + b1 0.

2 2 (3.3) Заметим, что краевая задача (2.1), (2.2) может быть сведена к задаче (3.1), (3.2) после умножения уравнения (2.1) на положительный множитель x p ( t )dt K ( x) = e a, (3.4) и тогда (x ) = K ( x )q( x), g (x ) = K ( x) f ( x ).

Идея вариационного метода состоит в замене краевой задачи (3.1), (3.2) равносильной задачей об отыскании дважды непрерывно дифференцируемой на [a, b] функции Y (x ), доставляющей экстремум следующему функционалу ( ) ( ) b J ( y ) = K ( x ) y 2 + ( x) y 2 + 2 g ( x) y dx + b y 2 (b ) 2Tb y (b) + (3.5) a ( ) + a y 2 ( a) 2Ta y (a) + 2 qb y(b) 2qa y (a), причем значения параметров a, b, q a, qb, Ta, Tb в этом функционале определяются в зависимости от значений a0, a1, a 2, b0, b1, b2 по таблице 3.1.

Таблица 3.1.

Значения параметров функционала a b № a0 a1 b0 b1 Ta Tb qa qb a2 b 0 1 0 0 0 0 0 a0 b b a2 0 0 K (b) 2 0 0 0 0 0 a0 b Продолжение таблицы 3.1.

a b № a0 a1 b0 b1 Ta Tb qa qb b a2 b 0 0 0 K (b) 3 0 0 0 a0 b0 b a b2 0 0 K (a) 4 0 0 0 0 0 b0 a a b 2 K (a) 2 K (b) 0 5 0 0 0 0 0 a1 b b0 a b2 2 K (a) 0 0 0 K (b) 6 0 0 0 b0 a b a a2 b2 0 0 0 K (a ) 7 0 0 0 a0 b0 a a b a2 0 K (a ) 0 0 0 K (b) 8 0 0 0 a0 a1 b a b a2 b2 0 K (a ) 0 0 0 0 K (b) 9 0 a0 b0 a1 b В методе Ритца для нахождения приближенного решения краевой задачи (3.1), (3.2) строится функциональная последовательнос ть {y n (x )} из пробных решений y n (x ) следующим образом.

Как и в методе Галеркина, задаемся на [a, b] функцией u 0 ( x ) и пробными функциями u1 (x ),..., u n ( x), такими, что u 0 ( x ) удовлетворяет условиям (3.2), а u1 (x ),..., u n ( x) удовлетворяют однородным условиям (2.3), и составляем функцию n C j u j ( x), y n (x ) = u 0 (x ) + (3.6) j = где С i ( i = 1, n ) – некоторые постоянные. Значения постоянных С i (i = 1, n ) подберем так, чтобы функция (3.6) доставляла экстремум функционалу (3.5).

Подс тавляя y ( x) = y n ( x) в (3.5), получаем квадратичную функцию переменных C 1,..., C n b 2 J ( yn ( x)) = K (x ) u0 (x ) + Ci ui (x) + ( x) u0 ( x) + Ci ui ( x) + n n a i =1 i = n n n + 2 g (x ) u 0 (x ) + Ci ui (x) dx + b u0 (b) + Ci ui (b) 2Tb u 0 (b) + C i ui (b) + i =1 i =1 i = n n + a u 0 (a) + C i ui (a) 2Ta u0 (a) + C i ui (a) + i =1 i = (3.7) n n + 2qb u0 (b) + C i u i (b) 2qa u 0 (a) + Ci u i (a ) = (C1, C 2,..., Cn ).

i =1 i = Необходимые условия экстремума функции (3.7), как известно из математического анализа, имеют вид:

= 0, k = 1, n, (3.8) C k Записав условия (3.8) в развернутом виде, для определения значений переменных C1, C 2,..., Cn получаем неоднородную систему линейных алгебраических уравнений n -го порядка n akj C j = bk, k = 1, n, (3.9) j = где a kj = (K ( x )u k u j + ( x)uk u j )dx + b u k (b)u j (b) + a u k (a )u j (a );

b a b bk = ( K ( x )u u + ( (x )u 0 + g (x ))uk )dx ( b (u 0 (b ) Tb ) + qb )u k (b) (3.10) 0k a ( a (u0 (a ) Ta ) qa )uk (a).

Решив систему (3.10) и подставив определяемые этим решением значения постоянных C1, C 2,..., Cn в (3.6), завершаем построение пробного решения y n (x).

Опишем теперь возможный алгоритм приближенного решения задачи (3.1), (3.2) методом Ритца, предполагая, что {y n (x)}0 сходится к Y (x ) при n.

1. Подготовительный шаг алгоритма. На этом шаге определяем значения параметров функционала (3.5) в соответствии с таблицей 3.1.

Выбираем функции u 0 ( x), u1 ( x),..., un (x ) и находим функцию R 0 ( x) = L(u0 ) g ( x ) = K ( x)u + K ( x)u 0 ( x) ( x )u 0 ( x) g ( x), т. е. невязку от подстановки u 0 ( x) в уравнение (3.1). Если x [a, b] : R 0 ( x) = 0, то u 0 ( x) = y 0 ( x) – искомое решение и вычисления заканчиваем. Если же R 0 ( x) 0, то переходим к следующему шагу алгоритма.

/ 2. Первый шаг алгоритма. Строим функцию y1 ( x ) = u0 ( x ) + C1u1 ( x), определив значение C1 из решения системы (3.9) при n = 1.

Находим невязку R (C1, x ) = L[ y1 ] g ( x ) = K ( x)(u 0 + C1 u1 ) + + K ( x)( u0 + C1u1 ) ( x)( u0 + C1u1 ) g ( x).

Если x [a, b] : R (C1, x) = 0, то Y (x ) = y1 (x ), и задача решена.

Если R (C1, x ) 0 на [a, b], то находим / max | y1 (x ) u 0 ( x) |= [ a,b ] или max | R (C1, x ) |= 12.

[ a,b ] Если 11 1 или 12 2, где 1 и 2 – заданные меры точнос ти приближенного решения, то полагаем Y (x ) y1 (x ) и вычисления заканчиваем.

Если же 11 1 или 12 2, то переходим к вычислениям на следующем шаге.

Таким образом, на m -м шаге ( m 1 ) алгоритма сначала строим функцию m y m (x ) = u 0 ( x) + C i ui (x ), i = определив значения C1, C 2,..., C m из решения системы (3.9) при n = m, а затем находим невязку m R (C1,..., C m, x) = L( ym ) g ( x ) = K (x ) u0 + C j u + j j = m m + K ( x) u0 + C j u j ( x) u 0 + C j u j g (x ).

j =1 j = Если x [a, b] : R (C1,..., C m, x) = 0, то Y (x ) = y m ( x). Если R (C1,..., C m, x ) 0, то / находим m1 = max | y m (x ) y m 1 ( x) | или m 2 = max | R (C1,..., C m, x )|. Если [ a,b ] [ a,b ] m1 1 или m 2 2, то Y (x ) y m ( x), если же m1 1 или m 2 2, то переходим к (m + 1) -му шагу, и т. д.

3.2. Построение систем пробных функций Некоторые методы подбора пробных функций были приведены в разделе 2.2.

Подчеркнем здесь, что если пробные функции выбираются на множестве многочленов, то их всегда можно найти методом неопределенных коэффициентов, причем неоднозначно. Одним из возможных наборов пробных функций u i (x ) ( i 1 ) будут многочлены u i ( x ) = ( x a) n + i 1 (b x) n (3.11) 1 или u i (x ) = ( x a) n (b x )n +i 1, (3.12) 1 где 1, если a1 = 0;

n1 = 2, если a1 0;

1, если b1 = 0;

n2 = 2, если b1 0.

Напомним также, что пробные функции u i (x ), i 1 можно выбрать из собственных функций соответствующей задачи на собственные значения, представляемых в виде u k ( x) = A k cos( k x) + B k sin( k x ). (3.13) Приведем теперь еще несколько примеров построения пробных функций для некоторых вариантов граничных условий (3.2).

Пример 1. Построить u 0 ( x ) и систему пробных функций для задачи (3.1) с краевыми условиями y (0) = 2, y (1) = 3.

Решение. Пусть u 0 = Bx + Cx 2. Тогда из граничных условий находим B = 2, B + 2C = 3, т. е. B = 2, C = 1 2. Итак, u0 ( x) = 2x + x 2.

Функции u k (x), k 1 будем искать в виде (3.13) u k ( x) = A k cos( k x) + B k sin( k x ).

Bk = 0, Удовлетворяя однородным граничным условиям, получим sin( k x) = 0. Отсюда имеем k = k, 2 = ( k ) 2, k = 1,2,.... Тогда, обозначая k C k = B k, находим u k ( x) = C k cos( kx). Таким образом пробное решение можно искать в виде y n (x ) = 2 x + x 2 + C k cos( kx).

n 2 k = Если же функции u k (x) ( k 1) взять в виде (3.11), то y n (x ) = 2 x + x 2 + C k (1 x ) 2 x k +1, n 2 k = если же в виде (3.12), то y n (x ) = 2 x + x 2 + C k (1 x ) k +1 x 2.

n 2 k = Пример 2. Построить u 0 ( x ) и систему пробных функций для краевой задачи с условиями y (0) = 2, y (1) = 4.

Решение. Положим u 0 ( x ) = A + Bx. Удовлетворяя граничным условиям, находим B = 2, A + B = 4, т. е. A = 2, B = 2, следовательно, u 0 ( x ) = 2 + 2 x.

Зададим u k ( x) = Ak cos( k x ) + B k sin( k x ).

Согласно однородным граничным условиям имеем 2k B k = 0, cos( k x ) = 0, k =, k = 1, 2, 3,...

2 Вводя обозначение C k = Ak, получаем 2k x, k 1.

u k (x ) = C k cos 2 Тогда 2k n x.

y n ( x ) = 2 + 2 x + Ck cos 2 k = Используя же пробные функции вида (3.11), (3.12), получаем n y n ( x) = 2 + 2x + Ck x k +1 (1 x) k = или n y n ( x ) = 2 + 2 x + Ck x 2 (1 x ) k.

k = Пример 3. Построить u 0 ( x ) и систему пробных функций для задачи с граничными условиями 2 y (0) = y (0) 1, y (1) = 3.

Решение. Ищем искомые функции в виде u 0 ( x ) = A + Bx, u k (x ) = Ak cos( k x ) + B k sin( k x), k 1.

Потребуем, чтобы функция u 0 ( x ) удовлетворяла неоднородным граничным условиям. Тогда 2B = A 1, A + B = 3, т. е. A = 7 / 3, B = 2 / 3, следовательно, u0 ( x) = + x.

Из соответс твующих однородных условий находим 2 k Bk Ak = 0, Ak cos( k ) + B k sin( k ) = 0.

Откуда получаем 2 k Bk = Ak и трансцендентное уравнение tg( k ) = 2 k (3.14) для определения собственных значений.

Последнее уравнение имеет счетное множес тво дейс твительных корней 1, 2,..., что подтверждает рисунок 3.1.

Корни уравнения (3.14) определяются приближенными численными методами, например, метод хорд, метод Ньютона, методом итерации или методом половинного деления. В системе MathCAD корни уравнений отыскиваются с помощью стандартной функции root (см. п. 8.2).

Таким образом, обозначая C k = Ak, имеем u (x ) = C k cos( k x) + sin( k x), 2 k тогда 72 n + x + C k cos( k x ) + sin( k x ).

yn ( x) = 2 k 33 k = Рис. 3.1. Геометрическая иллюстрация корней уравнения (3.14) Если использовать пробные функции вида (3.11), (3.12), то получаем 72 n y n ( x ) = + x + C k x k +1 (1 x ), 33 k = или 72 n yn ( x) = + x + C k x 2 (1 x) k.

33 k = 3.3. Задание к лабораторной работе Методом Ритца найти наиболее точное приближенное решение краевой задачи (2.16), построенное при помощи системы из n пробных функций – многочленов и системы из n пробных функций вида (3.13).

За меру точности выбрать (по указанию преподавателя) или max ym (x ) ym 1 ( x) = 1, [ a,b ] или max R (C1,..., C m, x ) = 2, [ a,b ] Варианты заданий приведены в таблице 2.1.

Лабораторная работа выполняется с использованием прикладной системы MathCAD, в которой реализуется алгоритм построения пробных решений y m (x) методом Ритца.

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

Числовые данные:

a, b – концы отрезка интегрирования [a, b ] ;

n – максимальное число параметров C1,..., C n в пробном решении.

Значение параметра n задает преподаватель.

Строчные данные:

аналитические выражения для пробных функций u 0 ( x ),..., um ( x).

В результате расчета программа выводит на экран дисплея значения коэффициентов C1,..., C n и таблицы всех пробных решений y1 (x ),..., y n ( x) и их невязок. Анализируя данные этих таблиц, надо найти обоснованный ответ на поставленную задачу лабораторной работы.

3.4. Выполнение работы в компьютерном классе 1. Прежде чем начать выполнение лабораторной работы на ЭВМ, внимательно ознакомьтесь с данной инструкцией.

2. При необходимости включите сами (или попросите лаборанта) питание компьютера. После того, как система загрузится, запускаем двойным щелчком левой кнопки мыши на рабочем столе программу Mathcad, если же ярлык отсутс твует, тогда открываем программу через кнопку «Пуск» (Программы Mathsoft Mathcad).


3. Узнайте у лаборанта расположение файла ODU.mcd и откройте его (File Open или, если программа русифицирована, Файл Открыть). При любой ошибке ввода программы нужно обратиться к лаборанту.

4. Прочитайте в начале файла задание на лабораторную работу и просмотрите пример выполнения работы, для которого исследование уже проведено. Программа файла ODU.mcd состоит из шести пунктов «1.

Постановка задачи», «2. Получение точного решения», «3. Получение приближенного решения методом Галеркина», «4. Получение приближенного решения вариационным методом Ритца», «5. Получение приближенного решения интегральным методом наименьших квадратов», «6. Выводы». Для выполнения лабораторной работы «Решение краевой задачи для линейного обыкновенного дифференциального уравнения второго порядка вариационным методом Ритца» необходимо использовать пункты 1, 2 и 4. Цели и задачи каждого из пунктов описаны ниже.

5. Для набора функций нужно либо воспользоваться всплывающим меню инструментов «Calculator», либо ввести ее с клавиатуры, используя следующие символы арифметических действий и стандартных функций: сложение – ‘+’;

вычитание – ‘–‘;

умножение – ‘*’;

деление – ‘/’;

возведение в степень – ‘^’;

квадратный корень – ‘\’;

синус – sin(x);

косинус – cos(x);

экспонента – exp(x);

натуральный логарифм – ln(x). При вводе числовых данных, являющихся десятичными дробями, целую и дробную части нужно разделять точкой (например, 0.5, 1.5 и т. д.).

6. Порядок выполнения работы Вам укажет программа подсказками и заданиями, выделенными красным цветом.

3.5. Порядок выполнения лабораторной работы Рекомендуется следующий порядок выполнения лабораторной работы.

1. Изучить разделы 3.1–3.4 настоящей главы и подготовить ответы на контрольные вопросы из раздела 3.7.

2. Пройти собеседование с преподавателем, получить допуск к выполнению работы на ПЭВМ, номер варианта задания лабораторной работы и значение параметра n.

3. В соответствии с полученным задания выполнить подготовительный шаг алгоритма метода Ритца и подготовить, если u 0 ( x) не является точным решением задачи, числовые и строчные исходные данные для расчетов на ПЭВМ.

4. Выполнить основную расчетную часть лабораторной работы в диалоге с ПЭВМ. В процессе диалога следует переписать с экрана дисплея значения коэффициентов C i пробных решений, а в конце диалога – итоговые таблицы пробных решений и их невязок.

5. Оформить и защитить отчет по лабораторной работе, который должен содержать титульный лист, математическую постановку задачи, результаты выполнения подготовительных расчетов, основные результаты расчетов на ПЭВМ и обоснованный итоговый отчет.

3.6. Тестирующий пример Методом Ритца найти на [0, 1] приближенное значение краевой задачи (2.17).

Сводим задачу (2.17) к задаче (3.1), (3.2), определяя, согласно (3.4), x K ( x) = exp ( 3)dt = e 3 x.

0 Получаем задачу (K ( x ) y ) (x ) y = g ( x ), y (0) + y (0) = 1, (3.16) y (1) + y (1) = 4, ( ) где K (x) = 3e 3 x ;

(x) = 2e3 x;

g(x) = 2x2 6x + 2 e3x. Напомним, что задача (3.16) имеет точное решение (2.18), значения которого представлены в таблице 2.2.

Для выполнения лабораторной работы используется тот же файл ODU.mcd, что и в разделе 2.5. В пункте «Постановка задачи» вводим числовые параметры и функции, входящие в задачу (2.17). В пункте «Получение точного решения в системе MathCAD» получаем таблицу 2.2 точного решения задачи.

В пункте «Получение приближенного решения вариационным методом Ритца» в качестве пробных функций используем те же функции, что и в разделе 2.5:

1 u0 ( x ) = 6 5 x;

u1 (x ) = 1 x + x 2 ;

u 2 ( x) = 1 x + x 3 ;

3 1 1 u 3 ( x ) = 1 x + x 4 ;

u 4 ( x ) = 1 x + x 5 ;

u5 ( x ) = 1 x + x 6 ;

5 6 причем R (u 0 ) 0.

/ Определяем параметры функционала (3.5). Так как a0 = a1 = b0 = b1 = 1, то в соответс твии с таблицей 3.1 имеем a b Ta = 2 = 1, Tb = 2 = 4, a0 b k (a ) = e = 1, a0 a = a b k (b ) = e = 0, 04979, b = b qa = qb = 0.

В результате расчета по программе ODU.mcd при n = 5 получим вектор коэффициентов С = (1,140938 2,564241 2,466128 0,133009 1,130778), следовательно, пробное решение имеет вид y 5 (x ) = u0 (x ) + 1,140938u1 (x ) 2,564241u 2 (x ) 2, 466128u3 (x ) + 0,133009u4 (x ) 1,130778u5 (x ).

Основные результаты расчета по программе ODU.mcd при n представлены в таблицах 3.2 и 3.3. В приложении A приведен пример при n = 5.

Анализируя таблицы, устанавливаем, что наилучшее приближение к точному решению Y (x ) дает решение y 5 (x ), для которого max Y ( x) y 5 (x ) 0,000031, [a,b ] max R (C1,..., C 5, x ) 0,017592, [a,b ] max y5 ( x) y 4 ( x) 0,000479.

[a,b ] Таблица 3.2.

Таблица пробных решений n=0 n=1 n=2 n=3 n=4 n= x 0,0 6,0 1,041255 0,835765 0,847800 0,846705 0, 0,1 5,5 1,020600 0,860600 0,865221 0,865260 0, 0,2 5,0 0,966887 0,893946 0,885499 0,886786 0, 0,3 4,5 0,880116 0,923337 0,904270 0,905463 0, 0,4 4,0 0,760287 0,936307 0,914450 0,914222 0, 0,5 3,5 0,607399 0,920389 0,906225 0,904175 0, 0,6 3,0 0,421453 0,863118 0,867057 0,864049 0, 0,7 2,5 0,202448 0,752026 0,781685 0,779613 0, 0,8 2,0 – 0,049615 0,574647 0,632119 0,633110 0, 0,9 1,5 – 0,334736 0,318515 0,397647 0,402690 0, 1,0 1,0 – 0,652915 – 0,028836 0,054831 0,061840 0, Таблица 3.3.

Таблица невязок пробных решений n=0 n=1 n=2 n=3 n=4 n= x 0,0 25,00 – 3,099555 1,276606 – 0,288980 0,045339 – 0, 0,1 24,58 – 1,569116 0,217300 0,038332 – 0,014853 0, 0,2 24,12 – 0,144793 – 0,490992 0,148810 – 0,015624 0, 0,3 23,62 1,173414 – 0,873203 0,115487 0,003438 – 0, 0,4 23,08 2,385504 – 0,954266 0,005948 0,017316 – 0, 0,5 22,50 3,491477 – 0,759113 – 0,117677 0,014424 0, 0,6 21,88 4,491334 – 0,312677 – 0,198705 – 0,004529 0, 0,7 21,22 5,385074 0,360109 – 0,185905 – 0,027676 – 0, 0,8 20,52 6,172698 1,234312 – 0,033496 – 0,033132 – 0, 0,9 19,78 6,854204 2,285000 0,298853 0,009867 – 0, 1,0 19,00 7,429595 3,487240 0,846022 0,139831 0, 3.7. Вопросы для самоконтроля 1. Опишите алгоритм решения краевой задачи для линейного дифференциального уравнения второго порядка аналитическим методом.

2. Каким образом уравнение (2.1) свести к равносильному уравнению типа (3.1)?

3. В чем основная идея вариационного подхода к решению краевой задачи (3.1), (3.2)?



Pages:   || 2 | 3 | 4 |
 





 
© 2013 www.libed.ru - «Бесплатная библиотека научно-практических конференций»

Материалы этого сайта размещены для ознакомления, все права принадлежат их авторам.
Если Вы не согласны с тем, что Ваш материал размещён на этом сайте, пожалуйста, напишите нам, мы в течении 1-2 рабочих дней удалим его.