Решение численными методами краевой задачи математической физики

Решение численными методами краевой задачи математической физики

1. Теоретико-аналитическая часть

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

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

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

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

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

— уравнение движения колебаний стержня

где ,

Граничные условия

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

Начальные условия

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

1.2 Вывод уравнения движения из основных законов физики

Стержень — упругое твёрдое тело, длина которого значительно превышает его поперечные размеры.

Рассмотрим стержень цилиндрической формы, на который действует вдоль оси стержня сила .

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

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

Рис. 1. Стержень

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

Найдём относительное удлинение участка стержня, ограниченного сечениями и .

Если абсцисса сечения в состоянии покоя, то смещение этого стержня в момент времени с точностью до бесконечно малых высшего порядка равно:

Отсюда ясно, что относительное удлинение стержня в сечении с абсциссой в момент времени выражается производной:

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

где — площадь поперечного сечения стержня, а модуль упругости (модуль Юнга) материала стержня.

Соответственно сила, действующая на сечение равна

Возьмем элемент стержня, заключённый между сечениями и. На этот элемент действуют силы и, приложенные в этих сечениях и направленные вдоль оси. Результирующая этих сил имеет величину и направлена также вдоль оси .

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

(1)

где объёмная плотность стержня, масса выделенного участка стержня

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

(2)

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

Если дополнительно предположить, что к стержню приложена внешняя сила, рассчитанная на единицу объёма и действующая вдоль оси стержня, то к правой части уравнения (1) добавится слагаемое и уравнение (1) примет вид:

(3)

(4)

это уравнение вынужденных продольных колебаний стержня.

1.3 Проверка задачи по критерию размерности

Вывод: размерности совпадают

1.4 Аналитическое решение задачи

Граничные условия:

Начальные условия:

Так как граничные условия ненулевые, использовать напрямую метод Фурье нельзя. С помощью введения новой переменной, приведём граничные условия к нулю:

тогда: граничные условия:

начальные условия:

частные производные:

.

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

граничные условия:

начальные условия:

В силу того, что задача неоднородна представим функцию в виде:

где функция будет описывать собственный колебаний стержня, а — вынужденные.

Собственные колебания

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

граничные условия:

начальные условия:

По методу Фурье решение можно представить в виде произведения функций, каждая из которых зависит только от одной переменной:

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

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

Решение дифференциального уравнения II порядка с постоянными коэффициентами ищем на основе характеристического уравнения:

Корни этого характеристического уравнения:

Следовательно, общее решение можно записать в виде:

Из граничных условий следует:

т.к., то

по условию и, т.к. следовательно,

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

Тогда собственные функции краевой задачи имеют вид:

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

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

Для нахождения решения этого уравнения составим характеристическое уравнение:

Корни этого характеристического уравнения:

Следовательно, общее решение можно записать в виде:

Каждому соответствует своё решение :

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

Пусть, тогда:

Для нахождения коэффициентов и используем начальные условия:

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

Итак,

Таким образом:

Вынужденные колебания

Рассмотрим задачу для, которая описывает вынужденные колебания стержня

граничные условия:

начальные условия:

Уравнение является неоднородным, но нулевые граничные условия позволяют строить решение в виде:

Пусть функция удовлетворяет условиям Дирихле на промежутке, тогда её можно представить в виде:

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

то есть:

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

Получили равенство двух линейных комбинаций, в этом случае коэффициенты равны:

Данное уравнение является неоднородным дифференциальным уравнением II порядка с постоянными коэффициентами. Общее решение такого уравнения складывается из суммы общего решения однородного уравнения и частного решения неоднородного уравнения II порядка

Общее решение неоднородного дифференциального уравнения II порядка с постоянными коэффициентами имеет вид:

Общее решение однородного дифференциального уравнения можно представить в виде:

Построение графиков приближенного решения при учете пяти и гораздо более пяти гармоник в среде MatLab. Для построения графика приближенного решения создан M-file (приложение 1), в котором записана функция, описывающая аналитическое решение задачи. Используя функцию из «analitic», построим графики приближенного решения, для крупной сетки и мелкой сетки при учёте 5 и 100 гармоник. Для это создадим M-file «a_reshenie» (приложение 2) и запустим его. Получим следующие графики:

Рис. 2. Аналитическое решение задачи для мелкой сетки при учёте 5 гармоник

Рис. 3. Аналитическое решение задачи для крупной сетки при учёте 5 гармоник

Рис. 4. Аналитическое решение задачи для мелкой сетки при учёте 100 гармоник

Из рисунков 2 и 4 видно, что количество гармоник не влияет на график функции, в отличие от масштаба сетки (рис. 2 и 3).

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

Рис. 5. Значения функции, вычисленной приближенно в узлах крупной сетки при учёте 5 гармоник

Найдём в заданной точке значение функции. Для это воспользуемся возможностями MatLab. А затем построим эту точку на графике аналитического решения для мелкой сетки при учёте 5 гармоник. Для этого создадим M-file «point» (приложение 3) и запустим его.

2. Дискретная модель

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

2.1 Построение дискретной модели и выбор сетки

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

Шаг по оси: ,

Шаг по оси: ,

где — длина стержня,

— рассматриваемый промежуток времени, — номера узлов по осям и соответственно,

— количество узлов по оси ,

— количество узлов по оси .

Получаем, что — граничные узлы; - внутренние узлы.

2.2 Разностная схема и разностная задача

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

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

.

3. Численное решение задачи методом «бегущего» счёта

3.1 «Ручной» счет методом «бегущего» счёта для крупной сетки

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

Тогда разностная схема будет выглядеть следующим образом:

Преобразовав, получим:

При расчёте сетки необходимо произвести сглаживание и допустить, что. Вычислим, чему равны остальные узлы сетки:

Таблица 1. Таблица значений вычисленной функции в узлах выбранной сетки

3.2 Численное решение методом «бегущего» счёта в среде MatLab

Для численного решения задачи необходимо создать M-file «numerical» (приложение 4), в котором функция вычисляет разностную задачу для крупной сетки. Для запуска данной функции необходимо создать файл «n_reshenie». При этом мы получим график.

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

Аналогичные значения мы уже получили при расчёте разностной задачи «вручную». Для более мелкой сетки создадим M-file «small_numerical» (приложение 6), и запускающий функцию M-file «n_reshenie_sm» (приложение 7).

Вьюненко Л.Ф., Бестужева А. Н. Применение численных методов для решения задач электрического транспорта железных дорог (уравнения с частными производными). Учебное пособие. СПб, 2003.

Араманович И.Г., Левин В. И. Уравнения математической физики. М., Наука, 1969.

Владимиров В. С. Уравнения математической физики. М., Наука, 1981.

Смирнов М. М. Уравнения в частных производных второго порядка. М., Наука, 1981.

Михлин С. Г. Курс математической физики. М., Наука, 1968.

Приложение 1

M-file «analitic»

function u1=analitic (x, t, n)

u1=t'*(x.^2)/2;

for k=1:n

u1=u1+16/((pi*(2*k+1))^3)'.*(((-1)^(k+1))*sin ((2*k+1)*pi*t/2)+t)'*

*sin ((2*k+1)*pi*x/2);

end

Приложение 2

M-file «a_reshenie»

clc

clear all

%% Построение графиков при учёте 5 гармоник (n=5)

% Для мелкой сетки:

x=0:0.01:1;% 1 — один метр

t=0:0.01:1;% 1 — одна минута

u1= analitic (x, t, 5);

figure (1);

mesh (u1);

%

% Для крупной сетки, (аналогичной сетке численного решения задачи):

x=0:0.2:1;% 1 — один метр

t=0:0.1:1;% 1 — одна минута

u2= analitic (x, t, 5);

figure (2);

mesh (u2);

u2

%

%% Построение графиков при учёте 100 гармоник (n=100)

% Для мелкой сетки:

x=0:0.01:1;% 1 — один метр

t=0:0.01:1;% 1 — одна минута

u3= analitic (x, t, 5);

figure (3);

mesh (u3);

Приложение 3

M-file «point»

clc

clear all

%% Построение графиков при учёте 5 гармоник (n=5)

% Для мелкой сетки:

x=0:0.01:1;% 1 — один метр

t=0:0.01:1;% 1 — одна минута

u1=analitic (x, t, 5);

figure (1);

mesh (u1);

xlabel ('x');

ylabel ('t');

zlabel ('u (x, t)');

title ('Вывод заданной точки x=0.81 (м), t=45 (сек) на график аналитического решения задачи')

%

x=0.81;% 81 сантиметр

t=45/60;% 45 секунд

u1=analitic (x, t, 5);

u1

x=82;% номер узла по оси x

t=76;% номер узла по оси t

hold on

plot3 (x, t, u1,'k*');

hold off

%

Приложение 4

M-file «numerical»

function u=numerical (l, t, a, h);

h=0.2;

tau=0.1;

a=1;

m=6;% по оси x

n=11;% по оси t

%

u=zeros (m, n);

%

for j=1:n

u (1, j)=0;

end

%

for i=1:m

u (i, 1)=0;

u (i, 2)=0;

end

%

u (m, 2)=0.01;

%

for j=3:n

for i=2:m-1

u (i, j)=0.25*(u (i+1, j-1)+u (i-1, j-1))+1.5*u (i, j-1) — u (i, j-2);

end

u (m, j)=u (m-1, j)+0.02*j;

end

u=u'

Приложение 5

M-file «n_reshenie»

u=numerical (1,1,1,0.2);

figure (1);

surf (u);

xlabel ('x');

ylabel ('t');

zlabel ('u (x, t)');

title ('Численное решение задачи для крупной сетки')

u

Приложение 6

M-file «small_numerical»

function u=small_numerical (l, t, a, h);

a=1;

tau=h²/(2*a²);% условие устойчивости

m=round (1+l/h);% по оси x

n=round (1+t/tau);% по оси t

%

u=zeros (m, n);

%

for j=1:n

u (1, j)=0;

end

%

for i=1:m

u (i, 1)=0;

u (i, 2)=0;

end

%

u (m, 2)=0.01;

%

for j=3:n

for i=2:m-1

u (i, j)=0.25*(u (i+1, j-1)+u (i-1, j-1))+1.5*u (i, j-1) — u (i, j-2);

end

u (m, j)=u (m-1, j)+0.02*j;

end

u=u';

Приложение 7

M-file «n_reshenie_sm»

clc

clear all

u=small_numerical (1,1,1,0.05);

figure (1);

mesh (u);

xlabel ('i');

ylabel ('j');

title ('Численное решение задачи для мелкой сетки')

уравнение стержень физика модель