Курсовой проект по курсы «Численные Методы» на тему ”Малые поперечные колебания неоднородной струны”

Казанский государственный университет

Факультет Вычислительной Математики и Кибернетики

Курсовой проект

по курсы «Численные Методы»

на тему

”Малые поперечные колебания неоднородной струны”.

Выполнила студентка гр. 956

Болдырева Юлия

Проверила Павлова М. Ф.

Казань, 2008

Вариант 3.

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

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

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

Задача на собственные значения [1]:

здесь -плотность струны.

Метод решения.

Задача (1)-(2) с помощью разностного методы сводится к дискретной задаче вида [2]:

Здесь - равномерная сетка на - сеточная функция, для

Вывод задачи на собственные значения.

,

, .

Решение ищется в виде произведения двух функций, одна из которых зависит только от x, другая - только от t (метод разделения переменных или метод Фурье):

.

Подставим в основное уравнение, получим: .

Равенство возможно лишь в том случае, если и левая и правая части не зависят ни от х, ни от t, т. е. представляют собой одну и ту же постоянную. Обозначим эту постоянную через – λ.

.

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

,

.

Выпишем задачу Штурма – Лиувилля (задачу на собственные значения) для второго уравнения:

, , , .

Вывод дискретной задачи (3)-(4).

Надо решить уравнение (1) для Можно приближенно заменить du/dx следующим образом:

(5)

где - расстояние между точками и . Если обозначить, соотвественно, i+1 точку сетки, значение функции в i-ом узле и i+1-ом узле то

(6)

получим правую разность.

При это разностное выражение стремится к du/dx.Но замена неоднозначна. Можно взять левую разность

(7)

или полусумму правой и левой – центральную разность:

(8)

Всюду здесь знак ~ означает соответствие или аппроксимацию. Говорят, что выражение аппроксимирует производную du/dx.

Аппроксимируем вторую производную тремя путями:

(9)

Уравнение (1) в соответствии с определением это есть разностное уравнение второго порядка. Его можно апроксимировать с помощью левой и правой разностей и записать в сделанных ранее обозначений в виде:

(10)

Рассмотрим уравнение

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

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

Складывая эти два выражения получим

В нашем случае задача о поперечных колебаниях струны – это краевая задача, когда дополнительные значения заданы при разных i, например, при i=0 задано (может быть также, например, задано при i=N ), т. е. требуется найти решение , 0<i<N, уравнения (1), если

(11)

где - заданные числа.

Множество точек (узлов) i=0,1,2,…,N называется сеткой. В граничных узлах сетки i=0 и i=N могут быть заданы не только значения функции, но и значения первой разности (как в нашем случае) или линейной комбинации функции и первой разности.

Общее условие можно записать в виде

(12)

В нашей задаче , т. е .

Подставляя

(13)

в первое из условий (12), получим

(14)

Случай означает, что в граничном узле i=0 задано значение функции (так называемое граничное условие первого рода). Если , то задано значение (граничное условие второго рода). В случае в точке i=0 задана линейная комбинация функции и первой разности (граничное условие третьего рода).

Наша задача – пример классической задачи с граничными условиями второго рода в точке и граничными условиями первого рода в точке .

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

(15)

Матрица системы уравнений (15) является трехдиагональной. Она имеет вид

(16)

Ее порядок равен N+1, если заданы краевые условия второго или третьего рода. Для системы уравнений (15), (11) имеем матрицу (N-1)-го порядка. У этой матрицы от нуля отличны только коэффициенты, стоящие на 3-х диагоналях – главной и двух соседних.

В нашем случае:

(17)

получаем (18)

Матрица системы A (порядок N+1=n):

(19)

Задача (3)-(4) является матричной задачей на собственные значения с симметричной матрицей порядка

(20)

где вектор - диагональная матрица, элементы которой

При вычислении первого собственного числа матрицы использовать метод обратной итерации:

пусть - произвольный вектор, образуем последовательности и

(21)

где

Условие выхода из итерационного процесса принять в виде

(22)

Замечание 1

На каждом шаге итерационного процесса вектор строится как решение системы которое вычисляется с помощью метода прогонки.

Метод прогонки

Метод прогонки = метод исключения Гаусса, который приводит к формулам прогонки, - эффективный метод решения для систем линейных уравнений с матрицами такого типа (A – трехдиагональная матрица).

Рассмотрим задачу

(23)

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

(24)

с неопределенными коэффициентами и . Выражение подставим в нашу систему

(25).

также (26)

Это уравнение выполнено для любых , если

(27)

Отсюда получаем рекуррентную формулу для :

(28)

и рекуррентную формулу для вычисления

(29).

Если коэффициенты и известны и известно значение , то двигаясь справа налево (от i+1 к i), мы определим последовательно все . Уравнения для и - нелинейные, они связывают значения этих функций в двух соседних точках. Для и задача решается слева направо, для - в противоположном направлении. Для каждой из функций надо решать задачу Коши. Чтобы найти начальные значения этих функций, используем граничные условия. Так как формула (24) справедлива при i=0,1,…,N-1, то при i=0 имеем

с другой стороны,

Поэтому (31,32 соответственно)

Таким образом, для функций получили задачи Коши: (28),(31) и (29),(32) (формулы прямой прогонки).

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

откуда, если то (33)

Таким образом, для определения получаем задачу Коши (24),(33) (формулы обратной прогонки). Формулы (24), (33) имеют смысл, если

(34)

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

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

В нашей задаче это выполняется, т. к.

Замечание 2

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

(35)

При этом нужно не забыть пронормировать полученный результат.

Тестовая задача.

Проведем тестирование метода при При этих данных известны все собственные числа задачи (3) –(4). Их значения равны

Будем изменять шаг, увеличивая количество разбиений (от 100 до 1000 через 30). Установив погрешность в (22) eps=0.0001

Графики зависимости погрешности вычисления собственных значений для тестовой задачи от выбранного шага

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

График зависимости нормы невязки от величины шага сетки для каждого собственного числа для тестовой задачи.

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

Графики зависимости собственных чисел от шага сетки для тестовой задачи.

Из графиков видно, что первого собственного значения алгоритм вычисляет довольно правильно уже при достаточно больших шагах (h=0.035), а вот для третьего и второго собственных значений требуется большее разбиение отрезка.

Графики собственных векторов для тестовой задачи при n=500 (шаг h = 0.0063)

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

Изменение погрешности для первого собственного значения в зависимости от изменения параметра в условии выхода из итерационного процесса

Условие выхода из итерационного процесса (22)

Построим графики для погрешности первого собственного значения для тестовой задачи для ,. Изменяем количество разбиений от 100 до 1000 через 30.

Количество итераций для первого собственного

Диапазоны итераций

1.для первого собственного значения

3

4

5

2.для второго собственного значения

4-5

6-8-9

12-13 -14

3.для третьего собственного значения

5-7

12-14

7,19-20-21

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

1.

2.

3.

Вычисление первых трех собственных чисел для неоднородной струны.

Исходные данные:

Вычисления проводим, изменяя количество разбиений от 100 до 1000 через 50.

График зависимости нормы невязки от величины шага сетки для каждого собственного числа.

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

Графики зависимости собственных чисел от шага сетки для тестовой задачи.

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

Графики собственных векторов для случая неоднородной струны.

Построим профили струны для исходных данных при разбиении отрезка на 500 частей (шаг h = 0.0063) – собственные вектора, соответствующие первым трем собственным значениям ( 0.3465, 3.4025, 9.5190).

Листинги с программой и результатами счета.

function ErrorEigenValue(lyamda, truelyamda)

%1/Графики зависимости погрешности вычисления собственных значений для тестовой задачи от выбранного шага

%2/График зависимости нормы невязки от величины шага сетки для каждого собственного числа.

%3/Графики зависимости собственных чисел от шага сетки.

k=1;

for j=100:30:1000

[truelyamda, shag(k)]=truetest(j);

[lyamda, shag(k),L, D,u]=test(j);

eigenvalue1(k)=lyamda(1);

eigenvalue2(k)=lyamda(2);

eigenvalue3(k)=lyamda(3);

eps1(k)=abs(lyamda(1)-truelyamda(1));

eps2(k)=abs(lyamda(2)-truelyamda(2));

eps3(k)=abs(lyamda(3)-truelyamda(3));

D=diag(D);

normanevyazki1(k)=norma((L-lyamda(1)*D)*u(:,1),diag(D));

normanevyazki2(k)=norma((L-lyamda(2)*D)*u(:,2),diag(D));

normanevyazki3(k)=norma((L-lyamda(3)*D)*u(:,3),diag(D));

k=k+1;

end

figure;%(1)

plot(shag, eps2,'b-');

xlabel('shag');

ylabel('eps1');

title('Zavisimost pogreshnosti 1 sobstvennogo znacheniya ot shaga dlya testovoi zadachi');

%hold on;

figure;%(2)

plot(shag, eps2,'b-');

xlabel('shag');

ylabel('eps2');

title('Zavisimost pogreshnosti 2 sobstvennogo znacheniya ot shaga dlya testovoi zadachi');

figure;%(3)

plot(shag, eps3,'b-');

xlabel('shag');

ylabel('eps3');

title('Zavisimost pogreshnosti 3 sobstvennogo znacheniya ot shaga dlya testovoi zadachi');

figure;%(4)

plot(shag, eigenvalue1,'b-');

xlabel('shag');

ylabel('eigenvalue1');

title('Zavisimost 1 sobstvennogo znacheniya ot shaga dlya testovoi zadachi');

figure;%(5)

plot(shag, eigenvalue2,'b-');

xlabel('shag');

ylabel('eigenvalue2');

title('Zavisimost 2 sobstvennogo znacheniya ot shaga dlya testovoi zadachi');

figure;%(6)

plot(shag, eigenvalue3,'b-');

xlabel('shag');

ylabel('eigenvalue3');

title('Zavisimost 3 sobstvennogo znacheniya ot shaga dlya testovoi zadachi');

figure;%(7)

plot(shag, normanevyazki1,'b-');

xlabel('shag');

ylabel('normanevyazki1');

title('Zavisimost nevyazki 1 sobstvennogo znacheniya ot shaga dlya testovoi zadachi');

figure;%(8)

plot(shag, normanevyazki2,'b-');

xlabel('shag');

ylabel('normanevyazki2');

title('Zavisimost nevyazki 2 sobstvennogo znacheniya ot shaga dlya testovoi zadachi');

figure;%(9)

plot(shag, normanevyazki3,'b-');

xlabel('shag');

ylabel('normanevyazki3');

title('Zavisimost nevyazki 3 sobstvennogo znacheniya ot shaga dlya testovoi zadachi');

function [lyamda, h,Anew, D,u]=test(KolRazbienii)

clc

m=3;%nahojdenie m sobstvennih znachenii

%c0=0.5;c1=0.1'%ro(x)=c0+c1(x)

c0=1;c1=0;%ro(x)=1

l=3.1415926535897;%dlina otrezka

eps=0.0001;%pogreshnost'=10 -4

n=KolRazbienii;%razbienie

u=zeros(n+1,m);%matrica sobstvennih vectorov

KolIteracii=zeros(m,1);%kolichestvo prodelannih iteracii

ksi1=0;

ksi2=1;

[lyamda, u,KolIteracii, h,Anew, D]=Zadacha(c0,c1,l, eps, n,m, ksi1,ksi2);

u;

KolIteracii;

%disp('wewewe');

SkalyarnoeProizvedenie(u(:,1),u(:,2));

function [truelyamda, h]=truetest(KolRazb)

%tochnoe reshenie testovoi zadachi

%ro(x)=1

m=3;%nahojdenie m sobstvennih znachenii

l=3.1415926535897;%dlina otrezka

n=KolRazb;%razbienie

h=l/n;

truelyamda=zeros(m,1);

truelyamda(1:m)=(4/(h^2))*(sin(((2*(1:m)-1)*h)/4)).^2;

function [lyamda, u,KolIteracii, h,Anew, D]=Zadacha(c0,c1,l, eps, n,m, ksi1,ksi2)

%nahojdenie m sobstvennih znachenii

%ro(x)=c0+c1(x)? ili celldata

h=l/n;%shag

w=(0:h:l)';%setka - indeksaciya

A=zeros(n-1,1);

C=zeros(n-1,1);

B=zeros(n-1,1);

for i=1:n-1

A(i,1)=-1/(h*h);

C(i,1)=-2/(h*h);

B(i,1)=-1/(h*h);

end%for i

D=zeros(n+1,1);

D(n+1)=(c0+c1*l)/2;%

D(1:n)=c0+c1*w(1:n);%

lyamda=zeros(m,1);

u=zeros(n+1,m);%matrica sobstvennih vectorov

NMax=1000;%maximal'noe kolichestvo iteracii

[lyamda(1),u(:,1),KolIteracii(1)]=MetodObratnoiIteracii(A, C,B, ksi1,ksi2,D, eps, NMax, u,0);

Anew=diag([1 - C' 1],0)+diag([A' - ksi2],-1)+diag([-ksi1 B'],1);

for i=2:m

[lyamda(i,1),u(:,i),KolIteracii(i,1)]=MetodObratnoiIteracii(A, C,B, ksi1,ksi2,D, eps, NMax, u,i-1);

end%for i

function [lyamda, u,KolIteracii]=MetodObratnoiIteracii(A, C,B, ksi1,ksi2,D, eps, NMax, u,kol)

%y - proizvolnii vector (posle ortogonalizacii???)

%lyamda -sobstvennoe chislo, u - sootvetstvuyushii sobstvennii vector

%eps=tochnost',NMax=maximalnoe kolichestvo itercacii

%u - sistema ortonormirovannih vectorov

%kol - kolichestvo ortogonal'nih vectorov v sisteme u

n=numel(D);

y=rand(n,1);%proizvolnii vector

y(1,1)=0;

x=running(A, C,B, ksi1,ksi2,D, y);

x = orto(x, u,D, kol);%ortogonalizaciya k postroennim i-1 vectoram iz u

y=x/norma(x, D);

Anew=diag([1 - C' 1],0)+diag([A' - ksi2],-1)+diag([-ksi1 B'],1);

lyamda0=SkalyarnoeProizvedenie(Anew*y, y);

for ij=1:NMax

x=running(A, C,B, ksi1,ksi2,D, y);

%y=x/norma(x, D);

x=orto(x, u,D, kol);%ortogonalizaciya k postroennim i-1 vectoram iz u

y=x/norma(x, D);

lyamda=SkalyarnoeProizvedenie(Anew*y, y);

if abs(lyamda-lyamda0)<=(eps*abs(lyamda0))

KolIteracii=ij;

break;

end

lyamda0=lyamda;

end%for ij

KolIteracii=ij;

u=y;

function x = running(A, C,B, ksi1,ksi2,D, y)

%metod progongi

%x:reshenie sistemi Ax=Dy

%(n*n)*(n*1)=(n*n)*(n*1)

%A - texdiagonalnaya matrica, D-diagonal'naya, y-stolbec

my=length(y);

x=zeros(my,1);

f=D.*y;%

alfa=zeros(my-1,1);

beta=zeros(my-1,1);

alfa(1,1)=ksi1;

beta(1,1)=0;%m'u1

for i=1:(my-2)

alfa(i+1)=B(i)/(C(i)-A(i)*alfa(i));

beta(i+1)=(f(i)+A(i)*beta(i))/(C(i)-A(i)*alfa(i));

end%for i

x(my)=ksi2*beta(my-1)/(1-alfa(my-1)*ksi2);

for k=my-1:-1:1

x(k,1)=alfa(k,1)*x(k+1,1)+beta(k,1);

end%for k

function n=norma(x, D)

%vichislyaet ||x||D

n=(SkalyarnoeProizvedenie(D.*x, x))^0.5;

function s=SkalyarnoeProizvedenie(v1,v2)

%vichislenie skalyarnogo proizvedeniya

%length(vector1)=length(vector2),vector1,vector2-stolbci

s=sum(v1.*v2);

function new = orto(beginvector, umatrica, D,index)

%orthonalizaciya nachal'nogo vectora po postroennim vectoram

%D-diagonal'naya matrica

%umatrica-postroennie vectora

%index-skol'ko vectorov postroenno

%beginvector - vector nachal'nogo priblijeniya

new=beginvector;

for j=1:index

new = new-(SkalyarnoeProizvedenie(beginvector, D.*umatrica(:,j)))*umatrica(:,j);

end%for j

Литература

[1] А. Н.Тихонов, А. А.Самарский Уравнения математической физики. – М.: Наука, 1972.

[2] А. А.Самарский Теория разностных схем. – М.: Наука, 1983.

[3] Б. Парлетт Симметричная проблема собственных значений. Численные методы. – М.: Мир, 1983.

Основные порталы, построенные редакторами

Домашний очаг

ДомДачаСадоводствоДетиАктивность ребенкаИгрыКрасотаЖенщины(Беременность)СемьяХобби
Здоровье: • АнатомияБолезниВредные привычкиДиагностикаНародная медицинаПервая помощьПитаниеФармацевтика
История: СССРИстория РоссииРоссийская Империя
Окружающий мир: Животный мирДомашние животныеНасекомыеРастенияПриродаКатаклизмыКосмосКлиматСтихийные бедствия

Справочная информация

ДокументыЗаконыИзвещенияУтверждения документовДоговораЗапросы предложенийТехнические заданияПланы развитияДокументоведениеАналитикаМероприятияКонкурсыИтогиАдминистрации городовПриказыКонтрактыВыполнение работПротоколы рассмотрения заявокАукционыПроектыПротоколыБюджетные организации
МуниципалитетыРайоныОбразованияПрограммы
Отчеты: • по упоминаниямДокументная базаЦенные бумаги
Положения: • Финансовые документы
Постановления: • Рубрикатор по темамФинансыгорода Российской Федерациирегионыпо точным датам
Регламенты
Термины: • Научная терминологияФинансоваяЭкономическая
Время: • Даты2015 год2016 год
Документы в финансовой сферев инвестиционнойФинансовые документы - программы

Техника

АвиацияАвтоВычислительная техникаОборудование(Электрооборудование)РадиоТехнологии(Аудио-видео)(Компьютеры)

Общество

БезопасностьГражданские права и свободыИскусство(Музыка)Культура(Этика)Мировые именаПолитика(Геополитика)(Идеологические конфликты)ВластьЗаговоры и переворотыГражданская позицияМиграцияРелигии и верования(Конфессии)ХристианствоМифологияРазвлеченияМасс МедиаСпорт (Боевые искусства)ТранспортТуризм
Войны и конфликты: АрмияВоенная техникаЗвания и награды

Образование и наука

Наука: Контрольные работыНаучно-технический прогрессПедагогикаРабочие программыФакультетыМетодические рекомендацииШколаПрофессиональное образованиеМотивация учащихся
Предметы: БиологияГеографияГеологияИсторияЛитератураЛитературные жанрыЛитературные героиМатематикаМедицинаМузыкаПравоЖилищное правоЗемельное правоУголовное правоКодексыПсихология (Логика) • Русский языкСоциологияФизикаФилологияФилософияХимияЮриспруденция

Мир

Регионы: АзияАмерикаАфрикаЕвропаПрибалтикаЕвропейская политикаОкеанияГорода мира
Россия: • МоскваКавказ
Регионы РоссииПрограммы регионовЭкономика

Бизнес и финансы

Бизнес: • БанкиБогатство и благосостояниеКоррупция(Преступность)МаркетингМенеджментИнвестицииЦенные бумаги: • УправлениеОткрытые акционерные обществаПроектыДокументыЦенные бумаги - контрольЦенные бумаги - оценкиОблигацииДолгиВалютаНедвижимость(Аренда)ПрофессииРаботаТорговляУслугиФинансыСтрахованиеБюджетФинансовые услугиКредитыКомпанииГосударственные предприятияЭкономикаМакроэкономикаМикроэкономикаНалогиАудит
Промышленность: • МеталлургияНефтьСельское хозяйствоЭнергетика
СтроительствоАрхитектураИнтерьерПолы и перекрытияПроцесс строительстваСтроительные материалыТеплоизоляцияЭкстерьерОрганизация и управление производством

Каталог авторов (частные аккаунты)

Авто

АвтосервисАвтозапчастиТовары для автоАвтотехцентрыАвтоаксессуарыавтозапчасти для иномарокКузовной ремонтАвторемонт и техобслуживаниеРемонт ходовой части автомобиляАвтохимиямаслатехцентрыРемонт бензиновых двигателейремонт автоэлектрикиремонт АКППШиномонтаж

Бизнес

Автоматизация бизнес-процессовИнтернет-магазиныСтроительствоТелефонная связьОптовые компании

Досуг

ДосугРазвлеченияТворчествоОбщественное питаниеРестораныБарыКафеКофейниНочные клубыЛитература

Технологии

Автоматизация производственных процессовИнтернетИнтернет-провайдерыСвязьИнформационные технологииIT-компанииWEB-студииПродвижение web-сайтовПродажа программного обеспеченияКоммутационное оборудованиеIP-телефония

Инфраструктура

ГородВластьАдминистрации районовСудыКоммунальные услугиПодростковые клубыОбщественные организацииГородские информационные сайты

Наука

ПедагогикаОбразованиеШколыОбучениеУчителя

Товары

Торговые компанииТоргово-сервисные компанииМобильные телефоныАксессуары к мобильным телефонамНавигационное оборудование

Услуги

Бытовые услугиТелекоммуникационные компанииДоставка готовых блюдОрганизация и проведение праздниковРемонт мобильных устройствАтелье швейныеХимчистки одеждыСервисные центрыФотоуслугиПраздничные агентства

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