Решение краевой задачи на графе методом Ритца

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

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

МИНОБРНАУКИ РОССИИ

ФЕДЕРАЛЬНОЕ ГОСУДАРСТВЕННОЕ БЮДЖЕТНОЕ ОБРАЗОВАТЕЛЬНОЕ УЧРЕЖДЕНИЕ ВЫСШЕГО ПРОФЕССИОНАЛЬНОГО ОБРАЗОВАНИЯ “ВОРОНЕЖСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ” (ФГБОУ ВПО «ВГУ»)

Факультет прикладной математики, информатики и механики

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

Дипломная работа

Направление/Специальность: 010501 Прикладная математика и информатика

Наименование специализации: Численные методы

Решение краевой задачи на графе методом Ритца

Допущена к защите в ГАК

Зав. кафедрой Т.М. Леденева, профессор, д.т.н.

Руководитель К.П. Лазарев, доцент, к.ф.-м.н.

Студент К.М. Горевалова

Воронеж 2012

Оглавление

Введение

1. Модельная задача уравнения колебаний струны

1.1 Вывод уравнения колебаний струны

1.2 Модельная задача деформации системы из трех струн

2. Вариационные методы

2.1 Поиск экстремума функционала

2.2 Минимизация функционала

2.3 Метод пробных функций

2.4 Метод Ритца

3. Решение модельной задачи

3.1 Подпространства сплайнов

3.2 Решение системы алгебраических уравнений

4. Описание программы и тестовые расчеты

Заключение

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

Приложение

Введение

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

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

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

1. Модельная задача уравнения колебаний струны

1.1 Вывод уравнения колебаний струны

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

Обозначим через u(x,t) отклонение струны от положения равновесия в момент t, предположив для определенности, что концы 0иlостаются закрепленными, так что

u(0,t)=u(l,t)=0.

Кинетическая энергия струны T, как сумма кинетических энергий ее частиц, выражается интегралом

, (1.1)

ритц колебание струна вариационный сплайн алгебраический

где ?dx есть масса элемента струны, отвечающего интервалу dx на оси х. Величина ?=?(х) есть плотность струны в точке х.

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

Определение 1.1. Струна есть одномерная механическая система, потенциальная энергия каждого участка которой пропорциональна его удлинению по сравнению с положением равновесия. [1]

Таким образом, мы имеем:

. (1.2)

Коэффициент р=р(х), фигурирующий в этом равенстве, называется модулем упругости струны (модулем Юнга).

Считая, что u_х -- настолько малая величина, что четвертой степенью u_х можно пренебречь, получим:

(1.3)

. (1.4)

Функция Лагранжа L = Т -- U имеет вид:

. (1.5)

Функционал Гамильтона теперь будет двойным интегралом:

. (1.6)

Напишем уравнение Эйлера - Остроградского:

. (1.7)

Если ? и p постоянны (т. е. струна однородна по плотности и упругости), , то мы получим уравнение:

, (1.8)

которое нас и интересовало.

Граничные условия, естественные с физической точки зрения, здесь можно взять следующие: при t=0 заданы значения функции u(х,0) (т. е. известна форма начального отклонения) и функции u_t(х,0) (известна начальная скорость каждой точки).

Покажем, что они определяют лишь единственное решение уравнения (1.7).

Если бы имелось два решения уравнения струны u⁽¹⁾(x,t) и u⁽²⁾(x,t), отвечающие одинаковым значениям

u⁽¹⁾(x,0)=u⁽²⁾(x,0)

и

u⁽¹⁾_(t)(x,0)=u⁽²⁾_(t)(x,0),

то их разность

u(x,t) = u⁽¹⁾(x,t) - u⁽²⁾(x,t)

была бы также решением, удовлетворяющим нулевым условиям

u(x,0)=0, u_t(x,0)=0.

Мы должны доказать, что u(x,t)?0. Используем для этого следующее соображение. Полная энергия струны

 (1.9)

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

u_t(x,0)=0,

u(x,0)=0,

так что

u_x(x,0)=0,

откуда при t=0 и E=0; а тогда в любой момент времени Е=0 и, следовательно,

u_t=u_x=0.

Отсюда следует, что и(х,t) постоянна; но так как

u(x,0)=0,

то и

и(х,t)=0

при каждом t.

Для струны справедлив закон сохранения энергии. Умножим уравнение струны (1.7) на и_t , и проинтегрируем по х:

. (1.10)

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

. (1.11)

Внеинтегральный член обращается в нуль, так как

u_t(0,t)=u_t(l,t)=0

и

u(0,t)=u(l,t)=0.

Таким образом, имеет место равенство:

, (1.12)

Откуда

(1.3)

и, следовательно, Е=const, что и требовалось доказать.

Рассмотрим некоторые частные случаи.

I. Стационарное линейное уравнение с учетом внешней нагрузки, сопротивления среды и при наличии пружины.

Такое уравнение имеет вид:

, (1.14)

А уравнение Эйлера:

 (1.15)

С условиями трансверсальности

, (1.16.1)

. (1.16.2)

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

, (1.17)

то есть

U= (1.18)

Тогда уравнение Эйлера имеет вид:

. (1.19)

При этом условия трансверсальности будут следующими:

, (1.16.1)

. (1.16.2)

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

, (1.20)

кинетическая энергия

, (1.21)

а функция Лагранжа примет вид:

(1.22)

или

. (1.23)

IV. Наконец, в следующем пункте более подробно рассмотрим нестационарную задачу на графе.

1.2 Модельная задача деформации системы из трех струн

Рассмотрим на примере задачу деформации системы из трех струн. Пусть точки О, A₁, A₂, A₃ (рис.1.1) находятся в одной плоскости. Рассмотрим механическую систему из трех струн, которые в положении равновесия совпадают с отрезками , , . Концы струн закреплены в точках A₁, A₂, A₃ и соединены между собой в точке O.

Рисунок 1.1. Механическая система из трех струн.

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

Отклонение -ой струны от положения равновесия обозначим , где x_i - абсцисса рассматриваемой точки струны, , -длина струны.

Для i-той струны (i=1,2,3) обозначим плотность внешней силы через , а натяжение - , аналогично рассмотренной задаче с одной струной.

В результате действия силы система деформируется и принимает следующий вид:

Рисунок 1.2. Деформированная система.

Найдем потенциальную энергию нашей системы трех струн. Эта энергия складывается из потенциальных энергий всех входящих в нее компонентов:

(1.24)

где ?_i =[0,l_i] - промежуток интегрирования, u_ix(x) - производная функции u_i(x), i=1,2,3.

Функция u(x)={u₁(x), u₂(x), u₃(x)} задает реальное положение струн в пространстве, и это реальное положение дает минимум в задаче минимизации функционала потенциальной энергии

(1.25)

в пространстве непрерывно-дифференцируемых функций v(x)={v₁(x), v₂(x), v₃(x)}, удовлетворяющих условиям:

v₁(0)= v₂(0)= v₃(0), (1.26)

v₁(l₁)=0, v₂(l₂)=0, v₃(l₃)=0. (1.27)

Будем выбирать произвольные непрерывно-дифференцируемые функции h_i(x) (i=1,2,3) и произвольный параметр ?, и строить новую функцию u_i(x)+?h_i(x), которая даст нам виртуальное расположение струн в пространстве.

Запишем энергию такой системы:

(1.28)

Очевидно, что при ?=0 виртуальное положение совпадает с реальным, потенциальная энергия которого минимальна. Отсюда следует, что выражение (1.28) имеет минимум при ?=0, поэтому производная этого выражения по ? в нуле равна нулю. Запишем:

 (1.29)

Выберем теперь функции h_i (i=1,2,3) так, чтобы на концах отрезка ?_i они были равны нулю, а h_j?0, i?j. Тогда для i=1,2,3:

(1.30)

Лемма 1.3. Лемма Дюбуа-Реймона. Пусть на [t₀, t₁] функции a(·) и b(·) непрерывны и пусть

(1.31)

для любой такой, что x(t₀)=x(t₁)=0. Тогда функция a(·) непрерывно дифференцируема и

. [3]

Пусть функции f и T_iu_ix (i=1,2,3) непрерывны, тогда воспользуемся леммой Дюбуа-Реймона. Т.к. интеграл (1.36) равен нулю , h(t₀)=h(t₁)=0 и функции f_i и T_iu_ix (i=1,2,3), то функция T_iu_ix и -f_i - (T_iu_ix)_x = 0 (1.32)

Таким образом, выражение (1.32) - это уравнение Эйлера для колебания струны.

Проинтегрируем теперь каждое слагаемое в выражении (1.29) по частям:

(1.33)

где д?_i - граница участка ?_i, i=1..3.

Предположим теперь, что функции h_i, i=1,2,3, в узле O принимают произвольное значение, а h_j?0, j?i, тогда:

(1.34)

Таким образом, т.к. u_1x(a)= u_2x(b)= u_3x(c)=0, то получим:

, (1.35)

откуда следует

 (1.36)

или

k=1,2,3. (1.37)

Итак, мы вывели нестационарное уравнение колебания:

(T_iu_ix)_x + f_i=0 (1.38)

или

(1.39)

Так же мы вывели условие, которое должны выполнятся в узлах соединения О:

(1.40)

Таким образом, краевая задача для описания формы отклонения системы имеет вид:

, , (1.41)

, (1.42)

, , (1.43)

. (1.44)

Эту задачу можно описать в операторном виде: Au=f или подробнее

Au(x)=f(x) (1.45)

Здесь х принадлежит графу, изображенному на рисунке 1.1 и состоящему из ребер ?₁, ?₂, ?₃; А - оператор, определенный формулой (1.46):

, . (1.46)

Оператор определен на дважды непрерывно дифференцируемых функциях

, (1.47)

удовлетворяющих условиям:

, (1.48)

, (1.49)

. (1.50)

2. Вариационные методы

2.1 Поиск экстремума функционала

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

Лемма 2.1. Если функционал F(y) дифференцируем при y=y₀, то при любом h функция F(y₀+ht), как функция от числа t, дифференцируема в обычном смысле по t и при t=0 ее производная равна L(y₀,h)=?F(y₀,h).

Линейный функционал L(y₀,h), определенный единственным образом, называется дифференциалом или, чаще, вариацией функционала F в точке y₀ и записывается ?F(y₀,h).

Рассмотрим функционал

(2.1)

в пространстве D₁(a,b) непрерывных функций на отрезке [a, b], обладающих непрерывными производными первого порядка. Ядро f(x,y(x),y'(x)) предполагается непрерывной функцией, обладающей непрерывными производными до второго порядка, определенной в области:

.

Составим приращение функционала F(y):

(2.2)

Вариация такого функционала существует и имеет вид:

(2.3)

Поставим задачу найти те точки, в которых функционал F(y) достигает экстремального значения - своего минимума или максимума.

Лемма 2.2. Во всякой точке y₀, где дифференцируемый функционал F(y) достигает экстремума, первая вариация ?F(y₀,h) функционала F при любом приращении h равна нулю.

Всякая точка y₀, в которой первая вариация функционала обращается в ноль при любом h, называется стационарной точкой функционала. Найдя все стационарные точки функционала и проанализировав их, выберем интересующие нас точки экстремума.

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

 (2.4)

Из условия h(a)=h(b)=0 следует, что внеинтегральный член равен нулю, и выражение вариации принимает вид:

(2.5)

Так как h(x) - любое, то h(x)=0 в выражении (2.5). Поэтому

(2.6)

Мы получили уравнение Эйлера. Итак, если экстремум функционала F существует и достигается на функции у=y(x), обладающей производной второго порядка, то эта функция у=y(x) удовлетворяет уравнению Эйлера. [1]

2.2 Минимизация функционала

Если каждой из функций y(x) из некоторого множества функций Y сопоставлено число Ф[y(x)], то говорят, что на множество Y задан функционал.

Задача минимизации функционала формулируется следующим образом:

-найти функцию , на которой функционал достигает своей точной нижней грани на этом множестве:

, y(x), . (2.7)

Иногда эту задачу называют минимизацией функционала по аргументу. Не всякий функционал и не на всяком множестве имеет минимум. Например, решение задачи (2.7) не существует, если функционал не ограничен снизу на заданном множестве, а так же в некоторых других случаях.

Определение 2.1. Оператор А называется симметричным, если он определен на плотном множестве и если для любых элементов ? и ? из области определения этого оператора имеет место тождество

.

Определение 2.2. Симметричный оператор А называется положительным, если для любого элемента имеет место равенство

,

причем выполняется только тогда, когда u=0.

Пусть требуется решить операторное уравнение

, (2.8)

Пусть оператор А аддитивен, положителен и симметричен, так что

(y,Ay)>0 при y?0,

(z,Ay)=(Az,y),

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

Рассмотрим функционал

 (2.9)

где y=y(x),

, (2.10)

а вид (y, Ay) зависит от вида оператора А.

Утверждение 2.1. Решение операторного уравнения (2.8) эквивалентно задаче на минимум функционала (2.9), и наоборот.

Доказательство. В самом деле, запишем произвольную функцию в виде:

(2.11)

Подставляя это выражение в правую часть формулы, получим:

(2.12)

Если - решение уравнения (2.8), то второе слагаемое правой части последнего уравнения обращается в 0, последний же член неотрицателен из-за положительности оператора А. Значит,

,

то есть функционал (2.9) достигает минимума на решении операторного уравнения (2.8).

Наоборот, если в представлении

есть функция, на которой функционал достигает минимума, то первая вариация функционала на этой функции равна 0, следовательно

,

какова бы ни была функция z(x). Применяя это условие к (2.12) и одновременно полагая:

,

получим:

что выполняется только при

.

Это означает, что функция, на которой функционал (2.9) достигает минимума, является решением операторного уравнения (2.8). Утверждение доказано.

Теорема 2.1. Теорема о минимуме квадратичного функционала. Пусть А - положительный оператор в D_A, . Пусть уравнение (2.8) имеет решение , то есть выполняется

в H, . (2.13)

Тогда квадратичный функционал

(2.14)

принимает на u₀ минимальное значение в D_A , то есть для всех имеет место , причем только при .

С другой стороны, пусть среди всех элементов функционал (2.14) принимает свое минимальное значение на элементе u₀. Тогда u₀ является решением уравнения (2.8) в H, то есть справедливо (2.13).

Теорема 2.2. Теорема Рисса. Любой линейный ограниченный функционал F в гильбертовом пространстве H может быть представлен в виде:

Fu=(u,v),

где v - некоторый элемент Н, однозначно определяемый функционалом F. Кроме того, ||v||=||F||. [4]

Если , , то , и (2.14) может быть записан в виде:

.

В данном случае - линейный ограниченный функционал, поэтому

.

Для : , отсюда . В соответствии с теоремой Рисса, такое, что . Отсюда получим:

.

То есть при u=u₀ можно записать следующее выражение: .

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

. (2.15)

Если , то ,

, то .

В таком случае будет являться решением уравнения

, и из оценки (2.15) получим: .

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

Если n-ное приближение решения u₀ таково, что , то имеем:

 (2.16)

То есть, если известна постоянная С, то соотношение (2.16) позволяет оценить погрешность приближенного решения . Но в то же время вопрос сходимости к f тесно связан с вопросом удачного выбора базиса.

2.3 Метод пробных функций

Общая схема численного решения заключается в сведении задачи (2.7) к поиску минимума функции многих переменных.

Рассмотрим класс V_n пробных функций заданного вида:

(2.17)

содержащих n свободных параметров и принадлежащих множеству V_n. На этом классе функций рассматриваемый функционал будет функцией n переменных - свободных параметров:

; (2.18)

найдя минимум функции F_n(a) и соответствующие ему значения параметров , мы определим функцию , на которой функционал достигает своего минимума в классе V_n.

Определение 2.3. Будем называть функционал Ф[y(x)] непрерывным, если он непрерывно зависит от функции y(x), то есть если фиксировать y(x), то для любого ?>0 найдется такое ?(?)>0, что при выполнении

 будет выполняться равенство: .

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

, (2.19)

где f - непрерывная функция всех своих аргументов.

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

Теорема 2.3. а) если система функций полная, а функционал Ф[y(x)] непрерывен, то последовательность является минимизирующей,

б) если требования пункта (а) выполнены и функционал удовлетворяет дополнительному условию

,?,?>0 (2.20)

то последовательность сходится к решению задачи (2.7).

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

2.4 Метод Ритца

Рассмотрим абстрактную схему метода Ритца для нахождения приближения к обобщенному решению операторного уравнения.

При зададим элементы , каждый из которых принадлежит пространству H_A. Обозначим через H_N линейную оболочку элементов

.

Будем считать, что выполнены следующие условия:

1) при элементы

2) линейно независимы;

3) последовательность подпространств {H_N } предельно плотна в H_А , то есть для существуют такие элементы , что

,

где - оценка аппроксимации и при n>0.

Будем искать приближение к обобщенному решению u₀ уравнения

Au=f при каждом n в виде .

 Далее подробно рассмотрим нахождение коэффициентов a_i.

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

(2.21)

(2.22)

(функционал Тихонова) или функционалу (2.9).

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

, (2.23)

то на них квадратичный функционал будет квадратичной функцией параметров a_k. Задача на нахождение минимума квадратичной функции F(a) посредством дифференцирования по переменным a_k сводится к системе алгебраических линейных уравнений, ее нетрудно решить численно. Этот частный случай метода пробных функций называется методом Ритца.

Выбор функций  целесообразно связать с краевыми условиями для задач типа (2.8), которые обычно линейны. Пусть например это условия первого рода:

(2.24)

Выберем какую-нибудь гладкую функцию  так, чтобы она удовлетворяла этим краевым условиям:

(2.25.1)

Или

(2.25.2)

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

(2.26.1)

Или

(2.26.2)

В этом случае пробные функции (2.23) при любых коэффициентах a_k удовлетворяют неоднородным краевым условиям (2.24) и являются полными на множестве непрерывных функций, удовлетворяющих этим краевым условиям.

Рассмотрим задачу на минимум квадратичного функционала с вещественным симметричным положительным оператором А:

(2.27)

Подставим в функционал пробные функции Ритца и получим квадратичную форму свободных параметров:

(2.28)

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

(2.29)

Распишем более подробно выражение (2.29):

(2.30)

Данную систему можно переписать в виде матричного уравнения , где

 , (2.31)

В зависимости от выбранных базисных функций ?_i , i=1,2,..n, матрица Ф, описанная в формуле (2.31), может иметь различный образ. Она может быть как полностью заполненной, так и разреженной, например диагональной.

Решив матричное уравнение, получим коэффициенты . Подставим их в выражение для приближенного решения (2.23). Элемент

назовем приближенным решением уравнения Ay=f ((2.8)) по Ритцу.

Теорема 2.4. Пусть для подпространства H_N выполнены условия 1) и 2). Тогда из последовательности приближенных решений уравнения (2.8) по Ритцу можно выделить минимизирующую подпоследовательность для функционала

((2.9)),

которая сходится по норме пространства H_A , и следовательно, по норме H, к обобщенному решению уравнения (2.8).

Заметим, что для не квадратичных функционалов Ф[y] линейные по параметрам пробные функции (2.23) не дают никаких преимуществ, так как получающиеся функции параметров

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

3. Решение модельной задачи

3.1 Подпространства сплайнов

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

Пусть на отрезке [a, b] задано разбиение

?: a=x₀<x₁<…<x_k=b. (3.1)

Для целого k?0 обозначим через C^k=C^k[a,b] множество k раз непрерывно дифференцируемых на [a, b] функций, а через C^-1[a, b] - множество кусочно-непрерывных функций с точками разрыва первого рода.

Определение 3.1. Функция S_n,?(x) называется сплайном степени n дефекта ? (0???n+1) с узлами на стеке ?, если

1) на каждом отрезке [x_i,x_i+1] функция S_n,?(x) является многочленом степени n, то есть:

Сплайн S_n,?(x) имеет непрерывные производные до n-? порядка, а производные более высших порядков, вообще говоря, терпят разрывы в узлах сетки.

Простейшим примером сплайна является единичная функция Хэвисайда:

,

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

.

Функции ?(x) и являются сплайнами соответственно нулевой степени и степени n дефекта 1 с единственным узлом в нулевой точке.

Теорема 3.1. Функции

; ,

линейно независимы и образуют базис в пространстве S_n,?(?) размерности

(n+1) + ?(N-1).

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

Расширим сетку, введя дополнительно точки

x_-n<…<x_-1<a, b<x_N+1<…<x_n+N . (3.2)

Возьмем функцию

?_n(x,t)=(-1)ⁿ⁺¹(n+1)(x - t)ⁿ

и построим для нее разделенные разности (n+1) порядка по значениям аргумента

t=x_i,…,x_i+n+1.

В результате получаются функции переменной x:

, i=-n,…,N-1. (3.3)

Лемма 3.1. Сплайны ), i=-n,…,N-1, обладают следующими свойствами:

а)

б)

Лемма 3.2. Функции ) являются сплайнами степени n дефекта 1 с конечным носителем минимальной длины.

Теорема 3.2. Функции), i=-n,…,N-1, линейно независимы и образуют базис в пространстве сплайнов S_n,1(?).

Докажем эту теорему. Покажем сначала линейную независимость функций

, i=-n,…,N-1,

на всей действительной оси. Предположим противное, то есть что существуют такие постоянные c_-n,…,c_N-1 не все равные 0, что

. (3.4)

Выбирая , получаем, что

и, значит,

c_-n=0.

Беря затем

, находим, что

c_-n+1=0,

и т.д., то есть имеем:

c_i=0, i=-n,…,N-1.

Следовательно, функции линейно независимы на .

Предположим теперь, что соотношение (3.4) выполняется только на [a, b]. Это значит, что на отрезках [x_i, x_i+1] обращаются в 0 сплайны вида:

i=0,…,N-1.

Каждый из них отличен от нуля самое большее на интервале (x_i-n, x_i+n+1). Поэтому из предположения

 при

следует что

при и ,

а значит и на всей действительной оси. В силу линейной независимости функций на  имеем c_p=0, p=i - n,…,i; i=0,…,N - 1.

Таким образом, функции линейно независимы, и так как по теореме 3.1 размерность пространства , то они образуют базис в этом пространстве. Теорема доказана.

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

, (3.5)

эти функции называют нормализованными В-сплайнами.

Тождество для нормализованных В-сплайнов имеет вид:

(3.6)

Распишем формулы (3.5) более подробно:

а) (3.7)

б) (3.8)

в) (3.9)

г) (3.10)

Первые 4 функции, высчитанные по формулам 3.7-3.10, изображены на рис. 3.1:

Рисунок 3.1. Примеры нормализованных В-сплайнов.

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

Дадим определение графа.

Определение 3.2. Геометрическим графом Г называется объединение точек всех ребер и некоторого (может быть пустого) подмножества вершин, которые мы называем внутренними. Вершины, не вошедшие в Г, называются граничными. [2]

Пусть Е(Г) - множество ребер , J(Г), ?Г, V(Г) - это множество внутренних, граничных и всех вершин графа Г соответственно. Для каждой вершины a?V(?) обозначим через E(a) множество всех ребер, примыкающих к а. Аналогично, для каждого ребра ??E(?) через V(?) обозначим множество вершин, примыкающих к ?.

Пусть Г⁰ - объединение точек всех ребер графа. Будем рассматривать функции u: Г>R и u: Г⁰>R, u_A(·) - сужение функции u(·) на множество А.

Для примера возьмем граф, изображенный на рисунке 3.2. Для построения такого графа были выбраны 4 точки, лежащие в одной плоскости - А₁, А₂, А₃ и О.

Рисунок 3.2. Пример графа.

Как видно из рисунка 3.2, рассмотренный граф имеет 3 ребра, все ребра имеют одну общую вершину.

Введем сетку на каждом ребре графа. Пусть длины ребер графа соответственно равны l₁, l₂, l₃, оси координат расположены так, как изображено на рисунке 3.3, то есть оси абсцисс совпадают с ребрами графа и направлены от центра, общая точка всех трех ребер - точка О - принята за начало координат, ось ординат направлена перпендикулярно плоскости графа вверх.



Рисунок 3.3. Система координат.

Обозначим i-тое ребро графа за OA_i и разобьем каждое ребро на n_i частей, i=1,2,3, получив таким образом три сетки с шагами

.

В результате таких действий получим набор точек x_ij, где i=1,2,3 и символизирует условный номер ребра графа, а j изменяется от 0 до n_i, где n_i - это число точек разбиения на i-том ребре (рис.3.4).

Таким образом, имеем набор точек на графе:

, (3.11)

где

, i=1,2,3, j=0,1,..,n_i . (3.12)



Рисунок 3.4. Сетка на графе.

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

, k=1,2,3; i=1,..,n_i - 1 (3.13)

В вершине О сплайн будет иметь особый вид, так как эта вершина инцидента сразу 3 ребрам - OA₁, OA₂ и OA₃, и сплайн должен располагаться одновременно на всех трех ребрах, поэтому он будет задаваться формулой (3.14):

 . (3.14)

Изобразив сплайны (3.9) и (3.10) на графе (3.2), получим следующую схему (рис. 3.5).

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

Рисунок 3.5. Схема сплайнов.

3.2 Решение системы алгебраических уравнений

В результате проведенной работы мы получили систему уравнений (1.56), где х принадлежит графу Г.

Au(x)=f(x) (1.56)

, . (1.57)

Оператор А имеет вид (1.57), пусть он аддитивен, положителен и симметричен, так что:

Согласно формуле (1.57),

. (3.15)

В главе 2 мы получили общую схему системы, с помощью которой будем решать исходную задачу. Система (2.22) имеет различные образы матриц в зависимости от того, какие функции будут взяты в качестве базисных. Наиболее удобными функциями являются линейные В-сплайны, так как в таком случае матрица Ф в системе (2.23) получится разреженной, почти трехдиагональной.

Рассмотрим подробнее выбор базисных функций.

Определение 3.5. Два элемента гильбертова пространства называются ортогональными, если их скалярное произведение равно нулю. Конечная или бесконечная последовательность попарно ортогональных элементов называется ортогональной системой. Если все элементы ортогональной системы нормированы, то система называется ортонормированной. [7]

Определение 3.6. Пусть Н - гильбертово пространство и {?_n} - ортонормированная в нем система. Мы будем называть эту систему полной в Н, если не существует элемента Н, кроме тождественного нуля, который был бы ортогонален ко всем элементам системы. [7]

Выберем функции ?_k(х), k=0,1,2,…, таким образом, чтобы они удовлетворяли однородным краевым условиям и при этом образовывали бы полную систему. Например, возьмем следующие функции:

, (3.16)

где k - номер ребра, i=1,..,n_i-1.

Для удобства дальнейшей работы перенумеруем сплайны, начиная от внутренней вершины вдоль ребер ?₁, ?₂, ?₃, и обозначим их через ?_i, i=1,..,n,

n = n₁ + n₂ + n₃ - 2.

После перенумерования базисные функции будут иметь вид:

(3.17)

где k - это номер текущего ребра графа, i = 1,2,.., n₁ + n₂ + n₃ - 2

Рисунок 3.6. Перенумерованные сплайны.

При этом система Ритца будет иметь вид Ay = b, где матрица А имеет портрет, представленный на рисунке 3.7.

Рисунок 3.7. Портрет матрицы.

Элементы матрицы, учитывая вид оператора А, указанный в формуле 1.57,

имеют следующий вид:

(3.18)

Для решения матричного уравнения, в котором матрица имеет портрет, изображенный на рисунке 3.7, воспользуемся специальным методом Гаусса. Будем последовательно исключать переменные, до тех пор, пока не получим верхнетреугольную матрицу. Так как матрица разрежена, то удобно хранить только ненулевые элементы. Здесь можно выделить несколько векторов: главная диагональ (1), диагонали над (2) и под (3) главной, а так же единичные ненулевые элементы над (4, 5) и под (6, 7) главной диагональю. При исключении элементов (3), (6) и (7) автоматически будут заполняться элементы, находящиеся непосредственно под элементами (4) и (5), поэтому необходимо хранить так же и эти вектора. Общая схема представлена на рисунке 3.8.

Рисунок 3.8. Схема хранения элементов матрицы.

4. Описание программы и тестовые расчеты

Метод, приведенный в теоретической части данной работы, реализован с помощью языка программирования Delphi 7.0. Разработанная программа позволяет найти приближенное решение краевой задачи на графе и для модельных задачи с известным точным решением вычислить погрешность приближённого решения.

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

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

На выходе пользователь получает несколько видов данных:

а) графики точного и полученного решений для каждой из трех струн, на графиках возможно масштабирование;

б) текстовый файл с сетками по каждому ребру графа;

в) текстовый файл с получившейся матрицей;

г) текстовый файл с правой частью матричного уравнения;

д) текстовый файл с получившимися коэффициентами для метода Ритца.

Основные процедуры и функции.

procedure TForm1.Button1Click(Sender: TObject); - ввод данных с формы. Вводятся длины струн, количество точек разбиения на каждой из струн

function u_toch(x: real; k: integer): real; - точное решение, найденное предварительно точными методами, необходимо для наглядного представления результатов работы программы; входные параметры : х - точка, значение в которой необходимо вычислить, k - номер рассматриваемой струны; возвращает точное значение в точке х

function f(x: real; k: integer): real; - входные данные, плотность струны; входные параметры : х - точка, значение в которой необходимо вычислить, k - номер рассматриваемой струны; возвращает значение в точке х

function Tt(t: real; k: integer): real; - входные данные, сила натяжения струны; входные параметры : х - точка, значение в которой необходимо вычислить, k - номер рассматриваемой струны; возвращает значение в точке х

function splain(t: real; k, i: integer): real; - описание сплайна 1 степени на струне номер k; входные параметры : t - точка, значение в которой необходимо вычислить, k - номер рассматриваемой струны, i - номер сплайна (соответствует номеру узла, являющегося вершиной) ; возвращает значение в точке t

function Аsplain(t: real; k, i: integer): real; - описание производной по t сплайна 1 степени на струне номер k; входные параметры : t - точка, значение в которой необходимо вычислить, k - номер рассматриваемой струны, i - номер сплайна (соответствует номеру узла, являющегося вершиной) ; возвращает значение в точке t

function skal_pr(p, r, k: integer): real; - вычисление скалярного произведения p-того и r-того сплайнов на k-той струне

function integral(p, k: integer): real; - вычисление правой части; p - номер сплайна, k - номер струны, возвращает значение

procedure sost_sist; - процедура для составления матрицы оператора; работает с глобальными параметрами, заполняет matr[i,j] элементами skal_pr(i,j,k) и вектор правой части

procedure resh_sist; - процедура для решения матричного уравнения методом Гаусса, вычисляет значения коэффициентов метода Ритца

procedure print_vect(vectors: string); - записывает в указанный в параметре файл значения всех векторов на момент вызова, нужна для промежуточного контроля

procedure TForm1.Button2Click(Sender: TObject); - вызывает функции для составления и решения системы линейных уравнений

procedure TForm1.Button3Click(Sender: TObject); - чертит графики точных и полученных решений.

Тестовый пример

Краевая задача для описания формы отклонения системы имеет вид (1.41)-(1.44):

, , (1.41)

, (1.42)

, , (1.43)

. (1.44)

Задачу можно описать в операторном виде: Au=f .

Таким образом, входными данными при решении подобной задачи будут:

1) Сила натяжения (на каждой струне)

function Tt(x: real; k: integer): real;

begin

case k of

1: Tt:=1;

2: Tt:=1;

3: Tt:=1;

end;

end;

2) Плотность сил

function f(k: integer): real;

begin

case k of

1: f:=-4/3;

2: f:=7/6;

3: f:=2/27;

end;

end;

Длины струн (задаются с клавиатуры):

Пример работы программы:



Рис.4.1. Пример работы программы при параметрах

l₁=1, l₂=2, l₃=3, n₁=n₂=n₃=4, koef=1000,

изображена струна 1

Рис. 4.2. Пример работы программы при параметрах

l₁=1, l₂=2, l₃=3, n₁=n₂=n₃=8, koef=1000

изображена струна 2



Рис. 4.3. Пример работы программы при параметрах

l₁=1, l₂=2, l₃=3, n₁=n₂=n₃=16, koef=1000

изображена струна 3

Таблица 4.1. Таблица погрешностей

Число точек разбиения	Абсолютная погрешность
4, 4, 4	5.20883837858082E-0002
8, 8, 8	2.96830985896908E-0002
16, 16, 16	1.57447173090176E-0002
32, 32, 32	8.09155296842234E-0003
64, 64, 64	4.66689785402547E-0003

Заключение

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

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

1. Шилов Г.Е. Математический анализ/ Г.Е. Шилов. - М.: Гос. изд-во физ.-мат. литературы, 1961. - 436 с.

2. Гайдай В.А. Моделирование и исследование сложных систем, параметризованных геометрическим графом: дис. канд. физ.-мат. наук: 05.13.18/ В.А. Гайдай. - Воронеж, 2009. - 79 с.

3. Алексеев В.М. Оптимальное управление/ В.М. Алексеев, В.М. Тихомиров, С.В. Фомин. - М.: Наука, 1979. - 432 с.

4. Ректорис К. Вариационные методы в математической физике и технике/ К. Ректорис. - М.: Мир, 1985. - 590 с.

5. Завьялов Ю.С. Методы сплайн-функций/ Ю.С. Завьялов, Б.И. Квасов, В.Л. Мирошниченко. - М.: Наука, 1980. - 352 с.

6. Уилсон Р. Введение в теорию графов/ Р. Уилсон. - М.: Мир, 1977. - 208 с.

7. Михлин С.Г. Вариационные методы в математической физике/ С.Г. Михлин. - М.: Наука, 1970. - 512 с.

8. Островский А.М. Решение уравнений и систем уравнений/ А.М. Островский. - М.: Изд-во ин. лит-ры, 1963. - 220 с.

9. Самарский А.А. Методы решения сеточных уравнений/ А.А. Самарский, Е.С. Николаев. - М.: Наука, 1978. - 592 с.

10. Березин И.С. Методы вычислений/ И.С. Березин, Н.П. Жидков. - М.: Гос. изд-во физ.-мат. литературы, 1959. - 620 с.

11. Самарский А.А. Численные методы/ А.А. Самарский, А.В. Гулин. - М.: Наука, 1989. - 432 с.

12. Алексеев В.М. Сборник задач по оптимизации/ В.М. Алексеев, Э.М. Галлеев, В.М. Тихомиров. - М.: Физматлит, 2007. - 255 с.

Приложение

Листинг

unit Unit1;

interface

uses

Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,

Dialogs, StdCtrls, ComCtrls, Series, TeEngine, ExtCtrls, TeeProcs, Chart;

type

TForm1 = class(TForm)

Label1: TLabel;

Label2: TLabel;

Label3: TLabel;

Label4: TLabel;

Button1: TButton;

Button2: TButton;

Button3: TButton;

Button4: TButton;

Edit1: TEdit;

Edit2: TEdit;

Edit3: TEdit;

Label5: TLabel;

Label6: TLabel;

Label7: TLabel;

Label8: TLabel;

Edit4: TEdit;

Edit5: TEdit;

Edit6: TEdit;

PageControl1: TPageControl;

TabSheet1: TTabSheet;

TabSheet2: TTabSheet;

TabSheet3: TTabSheet;

Chart1: TChart;

Series1: TLineSeries;

Series2: TLineSeries;

Chart2: TChart;

Chart3: TChart;

Series3: TLineSeries;

Series4: TLineSeries;

Series5: TLineSeries;

Series6: TLineSeries;

Label9: TLabel;

Edit7: TEdit;

Label10: TLabel;

procedure Button4Click(Sender: TObject);

procedure Button1Click(Sender: TObject);

procedure Button3Click(Sender: TObject);

procedure Button2Click(Sender: TObject);

private

{ Private declarations }

public

{ Public declarations }

end;

var

Form1: TForm1;

h: array [1..3] of real;

n: array [1..3] of integer;

l: array [1..3] of real;

x1: array [-1..101] of real;

x2: array [-1..101] of real;

x3: array [-1..101] of real;

z: array [-1..303] of real;

koef: integer;

matr: array [1..1000,1..1000] of real;

pr_ch, cf: array[1..1000] of real; // cf - коэффициенты разложения в методе Ритца

razb, razb_z, matrix, pravch, prav_ch_int: textfile;

implementation

{$R *.dfm}

procedure TForm1.Button4Click(Sender: TObject);

begin

close;

end;

procedure TForm1.Button1Click(Sender: TObject);

var

i: integer;

begin

l[1]:=strtofloat(edit1.Text);

l[2]:=strtofloat(edit2.Text);

l[3]:=strtofloat(edit3.Text);

n[1]:=strtoint(edit4.Text);

n[2]:=strtoint(edit5.Text);

n[3]:=strtoint(edit6.Text);

koef:=strtoint(edit7.Text);

for i:=1 to 3 do h[i]:=l[i]/n[i];

assignfile(razb,'razb.txt');

rewrite(razb);

assignfile(razb_z,'razb_z.txt');

rewrite(razb_z);

for i:=0 to n[1] do

begin

x1[i]:=h[1]*i;

z[i]:=x1[i];

writeln(razb,x1[i]);

writeln(razb_z,z[i]);

end;

writeln(razb,'');

writeln(razb_z,'');

x1[-1]:=x1[0]-h[1];

for i:=0 to n[2] do

begin

x2[i]:=h[2]*i;

writeln(razb,x2[i]);

end;

writeln(razb,'');

x2[-1]:=x2[0]-h[2];

for i:=0 to n[3] do

begin

x3[i]:=h[3]*i;

writeln(razb,x3[i]);

end;

x3[-1]:=x3[0]-h[3];

for i:=1 to n[2] do

begin

z[i+n[1]]:=x2[i];

writeln(razb_z,z[i+n[1]]);

end;

writeln(razb_z,'');

for i:=1 to n[3] do

begin

z[i+n[1]+n[2]]:=x3[i];

writeln(razb_z,z[i+n[1]+n[2]]);

end;

writeln(razb_z,'');

closefile(razb);

closefile(razb_z);

end;

function u_toch(x: real; k: integer): real;

begin

case k of

1: u_toch:=1/3-x+(2*x*x)/3;

2: u_toch:=1/3+x-(7*x*x)/12;

3: u_toch:=1/3-(x*x)/27;

end;

end;

function Tt(x: real; k: integer): real;

begin

case k of

1: Tt:=1;

2: Tt:=1;

3: Tt:=1;

end;

end;

function f(k: integer): real;

begin

case k of

1: f:=-4/3;

2: f:=7/6;

3: f:=2/27;

end;

end;

function splain(t: real; k, i: integer): real; // i- номер струны k- номер узла

begin

case i of

1:

if ((t> x1[k-1]) and (t<= x1[k])) then splain:=(t-x1[k-1])/(x1[k]-x1[k-1])

else

if ((t> x1[k]) and (t<= x1[k+1]))then splain:=(-t+x1[k+1])/(x1[k+1]-x1[k])

else splain:=0;

2:

if ((t> x2[k-1]) and (t<= x2[k])) then splain:=(t-x2[k-1])/(x2[k]-x2[k-1])

else

if ((t> x2[k]) and (t<= x2[k+1]))then splain:=(-t+x2[k+1])/(x2[k+1]-x2[k])

else splain:=0;

3:

if ((t> x3[k-1]) and (t<= x3[k])) then splain:=(t-x3[k-1])/(x3[k]-x3[k-1])

else

if ((t> x3[k]) and (t<= x3[k+1]))then splain:=(-t+x3[k+1])/(x3[k+1]-x3[k])

else splain:=0;

end;

end;

function Asplain(t: real; k, i: integer): real;

begin

case i of

1:

if ((t> x1[k-1]) and (t<= x1[k])) then Asplain:=1/h[1]

else

if ((t> x1[k]) and (t<= x1[k+1]))then Asplain:=-1/h[1]

else Asplain:=0;

2:

if ((t> x2[k-1]) and (t<= x2[k])) then Asplain:=1/h[2]

else

if ((t> x2[k]) and (t<= x2[k+1]))then Asplain:=-1/h[2]

else Asplain:=0;

3:

if ((t> x3[k-1]) and (t<= x3[k])) then Asplain:=1/h[3]

else

if ((t> x3[k]) and (t<= x3[k+1]))then Asplain:=-1/h[3]

else Asplain:=0;

end;

end;

// скалярное произведение. k - номер струны

function skal_pr(p, r, k: integer): real;

var

i, j: integer;

hh: array [1..3] of real;

skal_pr_1: real;

skal_pr_2, skal_pr_3: real;

begin

for i:=1 to 3 do hh[i]:=h[i]/koef;

skal_pr_1:=0;

case k of

0: begin

if (p=0) and (r=0) then

begin

skal_pr_1:=(Tt*Asplain(0,0,1)*Asplain(0,0,1)+

Tt*Asplain(x1[1],0,1)*Asplain(x1[1],0,1))/2;

for j:=1 to koef do

skal_pr_1:=skal_pr_1 + Tt*Asplain(j*hh[1],0,1)*Asplain(j*hh[1],0,1);

skal_pr_1:=skal_pr_1*hh[1];

skal_pr_2:=(Tt*Asplain(hh[2],0,2)*Asplain(hh[2],0,2)+

Tt*Asplain(x2[1],0,2)*Asplain(x2[1],0,2))/2;

for j:=2 to koef do

skal_pr_2:=skal_pr_2 + Tt*Asplain(j*hh[2],0,2)*Asplain(j*hh[2],0,2);

skal_pr_2:=skal_pr_2*hh[2];

skal_pr_3:=(Tt*Asplain(hh[3],0,3)*Asplain(hh[3],0,3)+

Tt*Asplain(x3[1],0,3)*Asplain(x3[1],0,3))/2;

for j:=2 to koef do

skal_pr_3:=skal_pr_3 + Tt*Asplain(j*hh[3],0,3)*Asplain(j*hh[3],0,3);

skal_pr_3:=skal_pr_3*hh[3];

skal_pr_1:=skal_pr_1+skal_pr_2+skal_pr_3;

end;

end;

1: begin

skal_pr_1:=( Tt*Asplain(0,p,k)*Asplain(0,r,k) +

Tt*Asplain(x1[n[1]],p,k)*Asplain(x1[n[1]],r,k) )/2;

for j:=1 to koef*n[1]-1 do

skal_pr_1:=skal_pr_1 + Tt*Asplain(j*hh[1],p,k)*Asplain(j*hh[1],r,k);

skal_pr_1:=skal_pr_1*hh[1];

end;

2: begin

skal_pr_1:=( Tt*Asplain(0,p,k)*Asplain(0,r,k) +

Tt*Asplain(x2[n[2]],p,k)*Asplain(x2[n[2]],r,k) )/2;

for j:=1 to koef*n[2]-1 do

skal_pr_1:=skal_pr_1 + Tt*Asplain(j*hh[2],p,k)*Asplain(j*hh[2],r,k);

skal_pr_1:=skal_pr_1*hh[2];

end;

3: begin

skal_pr_1:=( Tt*Asplain(0,p,k)*Asplain(0,r,k) +

Tt*Asplain(x3[n[3]],p,k)*Asplain(x3[n[3]],r,k) )/2;

for j:=1 to koef*n[3]-1 do

skal_pr_1:=skal_pr_1 + Tt*Asplain(j*hh[3],p,k)*Asplain(j*hh[3],r,k);

skal_pr_1:=skal_pr_1*hh[3];

end;

end;

skal_pr:=skal_pr_1;

end;

// правая часть ; p - номер сплайна, k - номер струны

function integral(p, k: integer): real;

var

i: integer;

int0, int: real;

hh: real;

begin

if p<>0 then

begin

hh:=h[k]/koef;

int0:=( f(k)*splain(0,p,k) + f(k)*splain(l[k],p,k) )/2;

for i:=1 to n[k]*koef-1 do int0:=int0 + f(k)*splain(i*hh,p,k);

int0:=int0*hh;

end;

if p=0 then

begin

int0:=0;

int:=0;

for k:=1 to 3 do

begin

hh:=h[k]/koef;

int0:=int0 +(f(k)*splain(0,p,k) + f(k)*splain(l[k],p,k))/2;

for i:=1 to n[k]*koef-1 do int0:=int0 + f(k)*splain(i*hh,p,k);

int0:=int0*hh;

int:=int+int0;

end;

int0:=int;

end;

integral:=int0;

end;

// составление матрицы

procedure sost_sist;

var

i, j, nn: integer;

vspom: array [1..100] of real;

begin

nn:=n[1]+n[2]+n[3]-2;

for i:=1 to nn do

for j:=1 to nn do matr[i,j]:=0;

assignfile(matrix,'matrix.txt');

rewrite(matrix);

write(matrix,'koef=');

writeln(matrix,inttostr(koef));

writeln(matrix,'');

matr[1,1]:=skal_pr(0,0,0);

for i:=2 to n[1] do

begin

for j:=2 to n[1] do if abs(i-j)<2 then matr[i,j]:=skal_pr(i-1,j-1,1);

end;

for i:=1 to n[2]-1 do

begin

for j:=1 to n[2]-1 do if abs(i-j)<2 then matr[i+n[1],j+n[1]]:=skal_pr(i-1,j-1,2);

end;

for i:=1 to n[3]-1 do

begin

for j:=1 to n[3]-1 do if abs(i-j)<2 then matr[i+n[1]+n[2]-1,j+n[1]+n[2]-1]:=skal_pr(i-1,j-1,3);

end;

matr[1,2]:=skal_pr(0,1,1);

matr[1,n[1]+1]:=skal_pr(0,1,2);

matr[1,n[1]+n[2]]:=skal_pr(0,1,3);

matr[2,1]:=matr[1,2];

matr[n[1]+1,1]:=matr[1,n[1]+1];

matr[n[1]+n[2],1]:=matr[1,n[1]+n[2]];

for i:=1 to nn do

begin

for j:=1 to nn-1 do write(matrix,matr[i,j]:10:4);

writeln(matrix,matr[i,nn]:10:4);

end;

closefile(matrix);

assignfile(pravch,'pravch.txt');

rewrite(pravch);

assignfile(prav_ch_int,'prav_ch_tolko_int.txt');

rewrite(prav_ch_int);

j:=1; // j - номер струны

for i:=0 to n[j]-1 do

begin

writeln(prav_ch_int,integral(i,j):7:4);

pr_ch[i+1]:=integral(i,j);

writeln(pravch,i,' ',pr_ch[i+1]:7:4);

end;

writeln(prav_ch_int,'');

writeln(pravch,'');

j:=2; // j - номер струны

for i:=1 to n[j]-1 do

begin

writeln(prav_ch_int,integral(i,j):7:4);

pr_ch[i+n[1]]:=integral(i,j);

writeln(pravch,i,' ',pr_ch[i+n[1]]:7:4);

end;

writeln(prav_ch_int,'');

writeln(pravch,'');

j:=3; // j - номер струны

for i:=1 to n[j]-1 do

begin

writeln(prav_ch_int,integral(i,j):7:4);

pr_ch[i+n[1]+n[2]-1]:=integral(i,j);

writeln(pravch,i,' ',pr_ch[i+n[1]+n[2]-1]:7:4);

end;

writeln(prav_ch_int,'');

writeln(pravch,'');

closefile(pravch);

closefile(prav_ch_int);

end;

procedure resh_sist; // решение системы

var

i, j, k, nn: integer;

a, b, c, d, e, v, w, y: array [1..100] of real;

vec1, vec2, vec3: textfile;

procedure gauss;

var

i, j: integer;

x: real;

begin

for i:=1 to n[1]-2 do

begin

c[i]:=c[i]/b[i];

d[i]:=d[i]/b[i];

e[i]:=e[i]/b[i];

y[i]:=y[i]/b[i];

b[i]:=1;

b[i+1]:=b[i+1]-c[i]*a[i];

d[i+1]:=d[i+1]-d[i]*a[i];

e[i+1]:=e[i+1]-e[i]*a[i];

y[i+1]:=y[i+1]-y[i]*a[i];

a[i]:=0;

v[i+1]:=v[i+1]-c[i]*v[i];

b[n[1]+1]:=b[n[1]+1]-d[i]*v[i];

e[n[1]+1]:=e[n[1]+1]-e[i]*v[i];

y[n[1]+1]:=y[n[1]+1]-y[i]*v[i];

v[i]:=0;

w[i+1]:=w[i+1]-c[i]*w[i];

w[n[1]+1]:=w[n[1]+1]-d[i]*w[i];

b[n[1]+n[2]]:=b[n[1]+n[2]]-e[i]*w[i];

y[n[1]+n[2]]:=y[n[1]+n[2]]-y[i]*w[i];