Решение краевых задач для обыкновенных дифференциальных уравнений методом конечных элементов

Размещено на http://www.allbest.ru/

МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ РОССИЙСКОЙ ФЕДЕРАЦИИ

Воронежский государственный технический университет

Физико-технический факультет

Кафедра высшей математики и физико-математического моделирования

КУРСОВАЯ РАБОТА

по дисциплине «Математика»

Тема: Решение краевых задач для обыкновенных дифференциальных уравнений методом конечных элементов

2010 г.

1. Теоретические сведения

1.1 Вариационный подход Ритца

Рассмотрим две задачи:

	(1.1)
	(1.2)

Эти задачи похожи: (1.1) является одномерным случаем более общей задачи (1.2). Уравнения (1.1) и (1.2) записаны в самосопряженной форме. Поставим задачам (1.1) и (1.2) в соответствие функционалы

	(1.3)
	(1.4)

Будем рассматривать пространство функций (пространство Соболева) с нормой

Это - функции с ограниченным интегралом.

Теорема: среди всех функций удовлетворяющих граничным условиям, решение задачи (1.1) придает наименьшее значение функционалу (1.3), а решение (1.2) - функционалу (1.4).

Доказательство.

Докажем это утверждение для одномерного случая, а доказательство для уравнений (1.2), (1.4) оставим в качестве упражнений.

Введем Поскольку а u(x) - дважды непрерывно дифференцируемая функция, то и о(0) = о(х) = 0.

	(1.5)

Третье слагаемое в (1.5) равно нулю в силу граничных условий для функции о; последнее слагаемое равно нулю, так как u - решение (1.1); второе слагаемое - неотрицательное. Следовательно, минимум функционала I₁(w) достигается, когда J(о) = 0, т. е. или, что то же самое, w(x) = u(x) .

Чуть сложнее эта теорема доказывается для двумерного случая, где надо воспользоваться теоремой Остроградского - Гаусса. Таким образом, решение соответствующей задачи в частных производных (1.2) или краевой задачи для ОДУ (1.1) сводится к задаче минимизации некоторого функционала.

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

1.2 Общая схема метода Ритца

Решение задачи (1.1) ищут в виде

	(1.6)

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

Подставив (1.6) в (1.3), получаем

	(1.7)

Находим минимум функционала (1.7) из условия

получаем систему из N линейных уравнений для определения коэффициентов C_k. Затем объявляем (1.6) решением задачи.

Точно так же поступаем и для функционала (1.4). Число уравнений в системе для определения коэффициентов тоже будет N (у базисной функции только один индекс!). Вид функционала будет аналогичен (1.7), но вместо интегралов по отрезку будут стоять двойные интегралы по рассматриваемой области пространства Щ, а вместо производных - градиенты.

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

Существуют два типа базиса: глобальный базис для метода Ритца и базис из функций с финитным носителем.

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

если вычисления проводятся точно.

Допустимый базис для применения в методе Ритца

q = 1, ..., N.

На отрезке [0, 1] допустимые базисы: где T_j(x) - j - й полином Чебышева;

Матрица системы линейных уравнений для определения коэффициентов разложения по базису метода Ритца получается заполненной. В случае использования "неудачных" базисов ее число обусловленности достаточно велико.

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

Рис. 1

Рассмотрим простейший вариант метода Ритца с использованием базиса функций с финитным носителем. Напомним, что носитель функции - множество точек x, для которых . Введем разбиение отрезка [0, X] точками x_j (сетку): 0 = x₀ < x₁ < ... < x_n = X. Строим базисные функции:



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

Пусть

тогда

а все остальные скалярные произведения равны нулю.

Также достаточно легко берутся интегралы, включающие в себя производные базисных функций. Носитель каждой такой базисной функции называется конечным элементом, а метод Ритца с использованием такого базиса - первый метод из семейства МКЭ. Иногда конечным элементом также называют саму базисную функцию с финитным носителем.

1.3 Пример построения схемы конечных элементов

Для уменьшения числа выкладок считаем, что

Рассмотрим несамосопряженный аналог задачи (1.1):

Найдем сопряженное уравнение:

Из этого соотношения легко получить условия, при которых . Теперь запишем разложение по базису:

подставляем в исходное дифференциальное уравнение:

Рассмотрим отдельно каждое слагаемое в левой части:

	(1.11)

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

В итоге после всех необходимых вычислений коэффициент перед C_k

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

Первые два слагаемые в правой части (1.11) зависят от граничных условий и относятся к правой части системы уравнений для определения C_k .

Рассмотрим теперь

Коэффициенты при C_k будут следующие:

Здесь опять предполагается, что функция q(x) кусочно - постоянная.

Последнее слагаемое с p(x) дает выражение

т.е. при "плохом" способе вычисляемых интегралов фактически получаем конечно - разностное соотношение, похожее на аппроксимацию Нумерова.

Вместе с тем, существует значительное отличие. Сеточная функция - это функция, заданная таблично. Решение (приближенное) МКЭ - это не сеточная функция, а элемент

1.4 Интерполяционный многочлен Лагранжа

Среди способов интерполирования наиболее распространен случай линейного интерполирования: приближение ищется в виде

,

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

, j=1,…,n. (1)

Метод решения задачи, при котором коэффициенты определяются непосредственным решением системы (1), называется методом неопределенных коэффициентов.

Наиболее изучен случай интерполирования многочленами

.(2)

Тогда

, i=1,…,n,

и система уравнений (1) имеет вид

j=1,…,n. (3)

Далее мы предполагаем, что все различные. Определитель этой системы отличен от нуля (определитель Вандермонда). Следовательно, система (3) всегда имеет решение, и притом единственное. Таким образом, доказано существование и единственность интерполяционного многочлена вида (2).

Непосредственное нахождение коэффициентов с помощью решения этой системы уже при сравнительно небольших n, например, при n=20, приводит к существенному искажению коэффициентов вычислительной погрешностью. Кроме того, уже сама запись многочлена в традиционной форме (2) часто приводит к большой вычислительной погрешности результата. При теоретических исследованиях, например при конструировании алгоритмов решения других задач, эти обстоятельства могут не играть роли. Однако при реальных вычислениях влияние вычислительной погрешности может быть недопустимо большим, и поэтому применяются другие виды интерполяционного многочлена и способы его записи.

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

Пусть есть символ Кронекера, определяемый соотношениями

Задача интерполирования будет решена, если удастся построить многочлены Ф_i(x) степени не выше n-1 такие, что Ф_i(x_j)= при i,j=1,…,n. Многочлен

будет искомым интерполяционным многочленом. В самом деле,

кроме того, g_n(x) - многочлен степени n-1. Поскольку =0 при , то делится на x-x_j при . Таким образом, нам известны n-1 делителей многочлена степени n-1, поэтому

.

Из условия Ф_i(x_i)=1 получаем

Интерполяционный многочлен (2), записанный в форме

(4)

называют интерполяционным многочленом Лагранжа.

Существуют другие формы записи того же интерполяционного многочлена(2), например рассматриваемая далее интерполяционная формула Ньютона и ее варианты. При точных(без округлений) вычислениях значения, получаемые по различным интерполяционным формулам, совпадают. Наличие же округлений приводит к различию в получаемых по этим формулам значений интерполяционных многочленов. Запись многочлена в форме Лагранжа, как правило, приводит к меньшей величине погрешности; запись же многочлена в форме Ньютона более наглядна и позволяет лучше проследить аналогию проводимых построений с основными построениями математического анализа. Кроме того, этим различным формам записи соответствует различное количество арифметических операций при вычислении с их помощью значений интерполяционного многочлена.

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

в точке х. вычисления можно проводить различными способами. Например, можно поступить следующим образом. Вычислить значение a1x и сложить с a₀. Далее вычислить значение a₂x² и сложить с полученным результатом и т.д. На j-м шаге, таким образом, вычисляется значение a_jx^j и складывается с уже вычисленной суммой a₀+a₁x+…+a_j-1x^j-1. Вычисление значения a_jx^j требует j операций умножения. Таким образом, описанный выше алгоритм требует для вычисления значения многочлена (1+2+…+n)=n(n+1)/2 операций умножения и n операций сложения. Количество арифметических действий в данном случае будет равно Ф₁=n(n+1)/2+n.

Ясно, что количество арифметических операций, необходимых для вычисления значения Р_n(x), может быть уменьшено. Например, можно последовательно вычислить и запомнить значения х², х³,…,хⁿ. Для этого потребуется n-1 операций умножения. Далее вычисляем величины a_jx^j (j=1,1…,n). Это потребует n операций умножения. Складывая полученные значения (это требует n операций сложения), получаем P_n(x). В этом случае Ф₂=(2n-1)+n и уже при n>2 имеет место неравенство Ф₂<Ф₁.

Можно пойти еще дальше. Запишем Р_n(x) в виде

P_n=(…((a_nx+a_n-1)x+a_n-2)x+…)x+a₀.

Для вычисления значения во внутренних скобках a_nx+a_n-1 требуется одна операция умножения и одна операция сложения. Для вычисления значения в следующих скобках (а_nx+a_n-1)x+a_n-2 требуется опять одна операция умножения и одна операция сложения, так как а_nx+a_n-1 уже вычислено, и т.д. Таким образом, вычисление P_n(х) при помощи этого алгоритма потребует n операций умножения и n операций сложения, то есть Ф₃=n+n.

Ясно, что Ф₃<Ф₂<Ф₁ при n>2. Таким образом, вычисление P_n(х) по последнему алгоритму потребует меньше арифметических операций и , соответственно, меньше времени ЭВМ. Количество арифметических операций, которое требуется для получения результата, является одной из важнейших характеристик метода, по которой происходит сравнение методов.

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

Ф₁=О(n²), Ф₂, Ф₃=O(n),

n - степень многочлена.

В последнем случае (Ф₂, Ф₃=O(n)) говорят, что методы имеют одинаковый порядок количества арифметических операций.

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

Ф₁=1\2n²+0(n²), Ф₂=3n+0(n), Ф₃=2n.

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

Заметим, что значение многочлена P_n(х) определяется параметрами а₀,а₁,…,а_n и величиной х. поэтому в общем случае для вычисления P_n(х) потребуется не менее n арифметических операций, т.е. мы имеем оценку снизу для количества арифметических операций. Таким образом, второй и третий методы вычисления P_n(х) являются оптимальными по порядку, так как Ф₂, Ф₃=O(n) и для любого метода Ф?n.

1.5 Интерполяционный многочлен Эрмита

В узлах x_k, принадлежащих [a,b], k=0,1,…,m, среди которых нет совпадающих узлов, заданы значения финкции f(x_k) и ее производных f⁽ⁱ⁾(x_k) до порядка N_k-1 включительно, i=1,2,…,N_k-1. Таким образом, в каждой точке x_k, k=0,1,…,m, известны



и следовательно, всего известно N₀+N₁+…+N_m величин. Требуется построить алгебраический многочлен H_n(x) степени n=N₀+N₁+…+N_m-1, для которого

(1)

Многочлен H_n(x_k), удовлетворяющий условиям (1), называется интерполяционным многочленом Эрмита для функции f(x). число N_k называется кратностью узла x_k.

Докажем, что интерполяционный многочлен Эрмита существует и единственен. Условия интерполяции (1) представляют собой систему линейных алгебраических уравнений относительно коэффициентов a₀,a₁,…,a_n многочлена

Число уравнений этой системы равно числу неизвестных и равно N₀+N₁+…+N_m. Поэтому достаточно сказать, что однородная система

(2)

имеет только тривиальное решение a₀=a₁=…=a_n=0.

Группа условий (2) при фиксированном k и i=0,1,…,N_k-1 означает, что число x_k является корнем кратности N_k многочлена H_n(x). Таким образом, многочлен H_n(x) имеет всего с учетом кратности не менее N₀+N₁+…+N_m=n+1 корня на [a,b]. Поскольку степень H_n(x) равна n, этот многочлен тождественно равен нулю, следовательно, равны нулю его коэффициенты и однородная система уравнений (2) имеет единственное решение a₀=a₁=…=a_n=0. Неоднородная система (1) однозначно разрешима при любых правых частях.

Поскольку значения f⁽ⁱ⁾(x_k), k=0, 1, …, m, i=0, 1, …, N_k-1, входят только в правую часть системы (1), коэффициенты a_j многочлена H_n(x) выражаются линейно через значения f⁽ⁱ⁾(x_k), и этот многочлен можно представить в виде линейной комбинации

где c_ki(x) - многочлены степени n.

Ввиду громоздкости выражений для c_ki(x) мы их не приводим.

Получим представление для погрешности интерполирования

Для этого рассмотрим вспомогательную функцию

(3)

где K - постоянная и

(4)

Постоянную К выберем так, чтобы в точке интерполирования х выполнялось условие g(x)=0, т.е. положим

Узлы x_k являются корнями кратности N_k функции g(s), k=1,2,…,m. Кроме того, точка x, принадлежащая [a,b], является корнем g(s). Таким образом, функция g(s) имеет с учетом кратности N₀+N₁+…+N_m+1=n+2 корня на отрезке [a,b]. По теореме Ролля производная g?(s) имеет, по крайней мере, один нуль между двумя соседними корнями функции g(s). Следовательно, g?(s) имеет не менее m+1 корня на [a,b] в точках, не совпадающих ни с одной из точек x₀,x₁, …, x_m, x. Кроме того g?(s) имеет в точке x_k корень кратности N_k-1, k=0, 1, …, m. Таким образом, g?(s) имеет с учетом кратности не менее

(N₀-1)+…+(N_m-1)+(m+1)=N₀+N₁+…+N_m=n+1

корней на [a,b]. Аналогично g??(s) имеет не менее n корней и т.д. Производная g⁽ⁿ⁺¹⁾(s) по крайней мере один раз обращается в нуль на [a,b], т.е. существует точка о, принадлежащая [a,b], в которой g⁽ⁿ⁺¹⁾(о)=0. Из (3) имеем

g⁽ⁿ⁺¹⁾(s)=f⁽ⁿ⁺¹⁾(s) - Kщ⁽ⁿ⁺¹⁾(s).

Так как щ(s) - многочлен степени n+1 со старшим коэффициентом 1, имеем щ⁽ⁿ⁺¹⁾(s)=(n+1)! Поэтому из условия g⁽ⁿ⁺¹⁾(о)=0 получаем, учитывая выражение для К, следующее представление для погрешности интерполирования:

.(5)

1.6 Иерархические многочлены

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



Рис. 2 - Одномерные элементы и ассоциируемые с ними иерархические базисные функции и аппроксимации: линейная (а), квадратичная (б) и кубическая (в)

многочлен линейный дифференциальный уравнение

Ясно, что линейное представление на этом элементе может быть получено только с помощью базисных функций, приведенных на рисунке 3, и не может быть улучшено, поскольку отождествление общих узлов смежных элементов гарантирует С⁰-гладкость.

(1)

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

(2)

c коэффициентами, выбранными таким образом, что при о = ±1, то требуемая С⁰-гладкость аппроксимации ц между элементами будет сохранена. Сформулированным условиям удовлетворяет симметричная парабола, соответствующим образом масштабированная.

Таким образом, для аппроксимации квадратным многочленом на элементе e (рис.2,б) можно записать

(3)

(4)

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

Аналогично для кубического элемента к квадратичному представлению (3) требуется добавить слагаемое , где - кубический многочлен вида

,(5)

принимающий нулевое значение при о = ±1 (т.е. в узлах 0 и 1). Из бесконечного числа возможностей выберем кубический многочлен, изображенный на рис.2, в. Он принимает нулевое значение в центре элемента, причем в этой точке . Можно сразу положить

(6)

поскольку такой кубический многочлен обладает требуемыми свойствами. Теперь параметр обозначает отклонение наклона в центре элемента от наклона предшествующей аппроксимации.

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

(7)

однако теперь придать физический смысл параметрам затруднительно (да и в этом нет особой необходимости).

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

(8)

где р(?2) - степень используемого многочлена [ и - те же, что и выше]. Это делает систему базисных функций элемента

(9)

Нетрудно видеть, что все производные от второго и более высоких порядков принимают нулевое значение при о=0, за исключением равной в этой точке единице. Следовательно, при использовании базисных функций вида (8) входящие в аппроксимацию параметры можно отождествить со значениями соответствующих производных:



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

2. Практическое задание

Задание.

Решить задачу одномерной теплопроводности

где к=1, Q=1 при 0?х?1/2 и Q=0 при1/2?х?1.

Решение.

Решим эту задачу методом конечных элементов в системе «Maple». Для каждого элемента строим аппроксимацию, используя многочлены Лагранжа:

Строим аппроксимацию:

>k0:=0: k1:=0.25: k2:=0.5:

>z1:=x->(x-k1)*(x-k2)*a0/((k0-k1)*(k0-k2))+(x-k0)*(x-k2)*a1/((k1-k0)*(k1-k2))+(x-k0)*(x-k1)*a2/((k2-k0)*(k2-k1));

> x0:=0.5: x1:=0.75: x2:=1.:

>z2:=x->(x-x1)*(x-x2)*a2/((x0-x1)*(x0-x2))+(x-x0)*(x-x2)*a3/((x1-x0)*(x1-x2))+(x-x0)*(x-x1)*a4/((x2-x0)*(x2-x1));

Для каждого конечного элемента составляем функционалы(вид этих функционалов - как в методе Ритца). Для первого элемента:

>F1:=z1->(1/2.)*int(diff(z1(x),x)^2-2*z1(x),x=0..0.5);

>W1:=F1(z1);

Для второго элемента:

>F2:=z2->(1/2.)*int(diff(z2(x),x)^2,x=0.5..1)-z2(1);

>W1:=F1(z1);

Cоставляем суммарный функционал:

>W:=W1+W2;

Из условия, что суммарный функционал должен достигать экстремума, получаем систему уравнений:

>p0:=diff(W,a0);p1:=diff(W,a1);p2:=diff(W,a2);p3:=diff(W,a3);p4:=diff(W,a4);

C учетом того, что (первое уравнение, которое соответствует нулевому узлу, при этом вычеркиваем):

> a0:=0;

> eqns:={p1=0,p2=0,p3=0,p4=0};

> r:=solve(eqns,{a1,a2,a3,a4});

r:= {a2 = .6249877528, a1 = .3437438764, a3 = .8749831600, a4 = 1.124981629}, {a4 = 1.125018371, a3 = .8750168405, a1 = .3437561238, a2 = .6250122476}

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

Для сравнения полученных и точных значений составим таблицу:

х		ц
0	0	0
0.25	0.3437438	0.34375
0.5	0.6249877	0.625
0.75	0.8749831	0.875
1	1.1249816	1.125

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

>z1:=unapply(subs(r,z1(x)),x);

>z2:=unapply(subs(r,z2(x)),x);

>y1:=-x^2/2+1.5*x;

>y2:=x+1/8;

>g0:=plot(y1(x),x=0..0.5,color=green, thickness=4):

>g1:=plot(z2(x),x=0.5..1,color=green, thickness=4):

> g2:=plot(z1(x),x=0..0.5):

> g3:=plot(z2(x),x=0.5..1):

> plots[display](g0,g1,g2,g3);



Рис. 3

Список литературы

1. Зенкевич О., Морган К. Конечные элементы и аппроксимация: Пер. с англ. - М.: Мир, 1986. - 318 с.

2. Самарский А.А., Гулин А.В. Численные методы: Учеб. пособие для вузов. - М.: Наука. Гл. ред. физ-мат. лит., 1989. - 432 с.

3. Численное решение краевых задач для обыкновенных дифференциальных уравнений методами конечных разностей и конечных элементов: Методические указания для самостоятельной работы студентов физико-технического факультета очной формы обучения по курсам “Математика”, “Методы математической физики”, “Математические методы моделирования физических процессов” / Воронеж. гос. техн. ун-т; Сост. С.А. Кострюков, В.Н. Нечаев, В.В. Пешков, Г.Е. Шунин. Воронеж, 2004. 35 с.

4. Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы с. 39-43.

Размещено на Allbest.ru