Решение численными методами краевой задачи математической физики
Решение численными
методами краевой задачи математической физики
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.
Численное решение задачи методом «бегущего» счёта
.1
«Ручной» счет методом «бегущего» счёта для крупной сетки
Пусть 
, 
, 
, 
, 
, 
, тогда из условия устойчивости получаем, что шаг по оси можно взять 
Выбранные значения параметров дают следующую точность вычисления
разностного решения:
Тогда разностная схема будет выглядеть следующим образом:
Преобразовав, получим:
При расчёте сетки необходимо произвести сглаживание и допустить,
что 
. Вычислим, чему равны остальные узлы сетки:
Таблица 1. Таблица значений вычисленной функции в узлах выбранной
сетки
|
|
|
|
|
|
|
|
|
0
|
0
|
0
|
0
|
0
|
0
|
|
|
0
|
0
|
0
|
0
|
0
|
0.0100
|
|
|
0
|
0
|
0
|
0.0025
|
0.0625
|
|
|
0
|
0
|
0
|
0.0006
|
0.0194
|
0.0994
|
|
|
0
|
0
|
0.0002
|
0.0058
|
0.0516
|
0.1516
|
|
|
0
|
0.0000
|
0.0017
|
0.0210
|
0.0973
|
0.2173
|
|
|
0
|
0.0005
|
0.0076
|
0.0504
|
0.1540
|
0.2940
|
|
|
0
|
0.0026
|
0.0225
|
0.0951
|
0.2197
|
0.3797
|
|
|
0
|
0.0090
|
0.0505
|
0.1527
|
0.2943
|
0.4743
|
|
|
0
|
0.0236
|
0.0937
|
0.2202
|
0.3785
|
0.5785
|
.2
Численное решение методом «бегущего» счёта в среде MatLab
Для численного решения задачи необходимо создать M-file «numerical» (приложение 4), в
котором функция вычисляет разностную задачу для крупной сетки. Для запуска
данной функции необходимо создать файл «n_reshenie». При этом мы получим
график.
А программа выведет нам значения функции в узлах сетки.
Аналогичные значения мы уже получили при расчёте разностной
задачи «вручную». Для более мелкой сетки создадим M-file «small_numerical» (приложение 6), и
запускающий функцию M-file «n_reshenie_sm» (приложение 7).
Список
литературы
1.
Вьюненко
Л.Ф., Бестужева А.Н. Применение численных методов для решения задач
электрического транспорта железных дорог (уравнения с частными производными).
Учебное пособие. СПб, 2003.
2.
Араманович
И.Г., Левин В.И. Уравнения математической физики. М., Наука, 1969.
3.
Владимиров
В.С. Уравнения математической физики. М., Наука, 1981.
4.
Смирнов
М.М. Уравнения в частных производных второго порядка. М., Наука, 1981.
5.
Михлин
С.Г. Курс математической физики. М., Наука, 1968.
Приложение
1
file
«analitic»
function u1=analitic (x, t, n)=t'*(x.^2)/2;k=1:n=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);
Приложение 2
file «a_reshenie»
clcall
%% Построение графиков при учёте 5 гармоник (n=5)
% Для мелкой сетки:=0:0.01:1;% 1 - один метр
t=0:0.01:1;% 1 - одна минута= analitic (x, t, 5);(1);
mesh(u1);
%
% Для крупной сетки, (аналогичной сетке численного решения
задачи):=0:0.2:1;% 1 - один метр=0:0.1:1;% 1 - одна минута
u2= analitic (x, t, 5);(2);
mesh(u2);
%% Построение графиков при учёте 100 гармоник (n=100)
% Для мелкой сетки:=0:0.01:1;% 1 - один метр
t=0:0.01:1;% 1 - одна минута= analitic (x, t, 5);(3);(u3);
Приложение 3
file «point»
clc
clear all
%% Построение графиков при учёте 5 гармоник (n=5)
% Для мелкой сетки:=0:0.01:1;% 1 - один метр
t=0:0.01:1;% 1 - одна минута=analitic (x, t, 5);(1);(u1);('x');('t');('u(x,
t)');
title ('Вывод заданной точки x=0.81 (м), t=45 (сек) на график
аналитического решения задачи')
%=0.81;% 81 сантиметр=45/60;% 45 секунд=analitic (x,
t, 5);
x=82;% номер узла по оси x=76;% номер узла по оси t
hold on(x, t, u1,'k*');off
%
Приложение 4
file «numerical»
function u=numerical (l, t, a, h);
h=0.2;=0.1;=1;=6;% по оси x=11;% по оси t
%=zeros (m, n);
%j=1:n(1, j)=0;
%i=1:m(i, 1)=0;(i, 2)=0;
%(m, 2)=0.01;
%j=3:ni=2:m-1(i, j)=0.25*(u (i+1, j-1)+u (i-1,
j-1))+1.5*u (i, j-1) - u (i, j-2);(m, j)=u (m-1, j)+0.02*j;=u'
Приложение 5
file «n_reshenie»
u=numerical
(1,1,1,0.2);(1);(u);('x');('t');('u(x, t)');
title ('Численное решение задачи для крупной сетки')
u
Приложение 6
file «small_numerical»
function u=small_numerical (l, t, a, h);
a=1;=h^2/(2*a^2);% условие устойчивости=round (1+l/h);% по
оси x
n=round (1+t/tau);% по оси t
%=zeros (m, n);
%j=1:n(1, j)=0;
%i=1:m(i, 1)=0;(i, 2)=0;
%(m, 2)=0.01;
%j=3:ni=2:m-1(i, j)=0.25*(u (i+1, j-1)+u (i-1,
j-1))+1.5*u (i, j-1) - u (i, j-2);(m, j)=u (m-1, j)+0.02*j;=u';
Приложение 7
file «n_reshenie_sm»
clcall=small_numerical
(1,1,1,0.05);(1);(u);('i');
ylabel('j');('Численное решение задачи для мелкой сетки')
уравнение стержень физика модель