Решение двумерного уравнения колебаний с помощью программной системы Matlab

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

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

МИНИСТЕРСТВО ОБРАЗОВАНИЯ РЕСПУБЛИКИ БЕЛАРУСЬ

Учреждение образования

«Гомельский государственный университет имени Франциска Скорины»

Физический факультет

Кафедра теоретической физики

Курсовая работа

Решение двумерного уравнения колебаний с помощью программной системы Matlab

Исполнитель

Н. Ю. Левчишина

Научный руководитель

Е.А. Дей

Гомель 2014



Реферат

Курсовая работа 42 страницы, 22 рисунка, 4 источника, 6 приложений

Ключевые слова: уравнение колебаний струны, одномерное уравнение теплопроводности, численное решение, основное уравнение электростатики, среда Matlab, команда pdepe, графический интерфейс

Объект исследования: уравнения математической физики.

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

Метод исследования: метод конечных разностей.

Полученные результаты: в результате был разработан интерфейс позволяющий решать уравнение теплопроводности.

Область применения: решение задач математической физики.

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



Содержание

Введение

1. Основные уравнения математической физики

1.1 Уравнение колебаний

1.2 Уравнение теплопроводности

1.3 Основное уравнение электростатики

2. Решение уравнений в частных производных с помощью Matlab

2.1 Численное решение граничной задачи для дифференциального уравнения 2-го порядка методом конечных разностей

2.2 Решение уравнения колебаний мембраны методом конечных разностей

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

3. Создание графического интерфейса в среде Matlab

3.1 Графические окна системы Matlab и элементы управления

3.2 Создание основных элементов управления

3.3 Callback-функции

4. Решение прикладных задач средствами Matlab

4.1 Создание интерфейса

4.2 Зависимость решений уравнения теплопроводности от заданных параметров

Заключение

Список использованных источников

Приложение А. Пример m-файла численного решения граничной задачи для дифференциального уравнения 2-го порядка методом конечных разностей

Приложение Б. Пример m-файла решения уравнения колебаний мембраны методом конечных разностей

Приложение В. Пример m-файла решения уравнения теплопроводности с помощью встроенной функции pdepe

Приложение Г. Текс m-файла создающего интерфейс

Приложение Д. Текс m-файла создающего интерфейс позволяющий решать уравнение теплопроводности

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



Введение

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

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

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

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

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

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

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

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



1. Основные уравнения математической физики

1. 1 Уравнение колебаний

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

. (1.1)

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

Рисунок 1.1 Элемент струны

По формуле квадрата длины элемента дуги имеем:

d

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

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

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

. (1.2)

Легко видеть, что

и .

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

.

Следовательно, правая часть уравнения (1.2) сводится к выражению

.

Разложим теперь в ряд Тейлора

.

Тогда с точностью до малых второго порядка

. (1.3)

Подставляя (1.3) в (1.2), приходим к равенству

. (1.4)

Введя обозначение

, (1.5)

получаем уравнение свободных колебаний струны:

. (1.6)

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



, . (1.7)

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

, . (1.8)

Таким образом, колебания струны описываются одномерным волновым уравнением (1.6). Если поперечные колебания совершает натянутая упругая пленка (мембрана), то соответствующее волновое уравнение является

двумерным:

. (1.9

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

, .

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

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

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

(1.10)

или (в сокращенной записи):

.

Здесь -- потенциал скоростей движения точек среды, - скорость звука в данной среде.

Начальные условия записывают в такой форме:

, .

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

.

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

,

,

где скорость света в вакууме.



1.2 Уравнение теплопроводности

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

, (1.11)

где -- коэффициент (внутренней) теплопроводности вещества.

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

. (1.12)

Происхождение знака минус в (1.11) и (1.12) понятно: градиент направлен в сторону возрастания температуры, а тепло течет к более холодным точкам тела.

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

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

(1.13)

Часть этого тепла идет на повышение температуры элемента , а остальная доля из-за теплопроводности уйдет в окружающие слои тела.

Сначала определим . Количество тепла, необходимое для повышения температуры бесконечно малого элемента от до , равно:

,

где удельная теплоемкость тела, его плотность.

Интегрируя это равенство по объему , получим:

(1.14)

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

ограничивающую объем , протекает количество теплоты

Поэтому

(1.15)

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

. (1.16)

В этом равенстве стоят интегралы по разным переменным. Поэтому применим к последнему из них теорему Остроградского-Гаусса:

Тогда равенство (1.16) примет вид:

или

Поскольку это соотношение справедливо для произвольного объема , то должно быть равным нулю само подынтегральное выражение:

.

Подставляя сюда значение из (1.12), получаем:

. (1.17)

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

.

Учитывая, что , приходим к следующему дифференциальному уравнению распространения тепла:

. (1.18)

Рассмотрим частные случаи этого уравнения.

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

то уравнение (1.18) принимает более простой вид:

, (1.19)

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

2. Установившийся поток тепла.

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

, (1.20)

где .

3. Установившийся поток тепла без тепловыделения.

В этом случае и и , поэтому распределение температуры в теле подчиняется уравнению Лапласа:

. (1.21)

Уравнение нестационарной теплопроводности

,

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

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

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

. (1.22)

Что касается граничных, или краевых, условий, то при решении физических задач они бывают трех видов.

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

. (1.23)

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

,

где единичный вектор нормали к поверхности.

Так как согласно (1.12) , то граничные условия имеют дифференциальный характер:

. (1.24)

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

. (1.25)

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

.

Множитель пропорциональности называется коэффициентом теплоотдачи.

Таким образом, тепловой поток , вытекающий из тела наружу, равен:

. (1.26)

С другой стороны, такой же тепловой поток должен подводиться изнутри путем теплопроводности. Поэтому согласно (1. 12)

. (1.27)

Приравнивая правые части (1.26) и (1.27) получаем:

.



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

.

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

1.3 Основное уравнение электростатики

Основным физическим законом электростатического поля является теорема Гаусса.

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

(1.28)

В общем случае электрические заряды распределены по объему с некоторой плотностью . Поэтому вместо суммы в правой части (1.28) появляется интеграл:

. (1.29)

В обеих частях этого равенства интегрирование производится по разным переменным. Поэтому, применив к (1.29) теорему Остроградского-Гаусса:

получаем:

(1.30)

Так как объем в (1.30) является произвольным, то равно нулю само подынтегральное выражение, и мы переходим к дифференциальной форме теоремы Гаусса:

, (1.31)

представляющей собой в теории электричества третье уравнение Максвелла.

Из электродинамики известно, что электростатическое поле потенциально:

, (1.32)

где электрический потенциал. Учитывая это, можно уравнение (1.32) записать в таком виде:

. (1.33)



Это основное дифференциальное уравнение электростатики уравнение Пуассона.

В системе единиц СИ это уравнение записывается проще:

. (1.34)

С его помощью можно установить вид скалярного поля потенциала , если известно распределение зарядов в пространстве .

Если же в некоторой области зарядов нет (), то потенциал удовлетворяет уравнению Лапласа:

. (1.35)



2. Решение уравнений в частных производных с помощью Matlab

2.1 Численное решение граничной задачи для дифференциального уравнения 2-го порядка методом конечных разностей

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

, (2.1)

где - координата;

- искомая функция;

, - некоторые непрерывные функции координаты.

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

. (2.2)

Зададим на отрезке равномерную координатную сетку с шагом :

. (2.3)

Граничные условия первого рода (условия Дирихле) для рассматриваемой задачи могут быть представлены в виде

, (2.4)

, (2.5)



где , - координаты граничных точек ;

, - некоторые константы.

Граничные условия второго рода (условия Неймана) для рассматриваемой задачи могут быть представлены в виде

, (2.6)

. (2.7)

Проводя дискретизацию граничных условий Дирихле на равномерной координатной сетке (2.3) с использованием метода конечных разностей, получим

, (2.8)

, (2.9)

где, - значения функции в точках , соответственно.

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

, (2.10)

. (2.11)

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

, (2.12)

где , - значения функций , в точке сетки с координатой .

Таким образом, в результате дискретизации получим систему линейных алгебраических уравнений размерностью , содержащую уравнения вида (2.12) для внутренних точек области и уравнения (2.8) или (2.10) и (2.9) или(2.11) для двух граничных точек.

Приведем один из вариантов m-файла для численного решения уравнения (2.2) с граничными условиями (2.4) - (2.7) на координатной сетке (2.3). Текст m-файла см. в приложении A.

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

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

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

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

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

Например, при запуске m-файла и вводе исходных данных

x0=0;

xn=5;

n=60;

f='2*sin(x.^2)+cos(x.^2)';

v1=1;

g1=0;

v2=1;

g2= -0.5;

на экране в отдельных графических окнах появятся график функции правой части уравнения Пуассона (рисунки 2.1) и график искомой функции (рисунок 2.2). Точками на графиках будет отображаться координатная сетка.

Рисунок 2.1 График функции правой части уравнения (2.2)

Рисунок 2.2 График искомой функции уравнения (2.2)

2.2 Решение уравнения свободных колебаний мембраны методом конечных разностей

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

Волновое уравнение относится к уравнениям гиперболического типа и в двухмерном случае имеет вид

, (2.13)

где - время;

, - координаты;

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

Зададим на отрезке равномерную координатную сетку с шагом

, (2.14)

на отрезке равномерную координатную сетку с шагом

, (2.15)

на отрезке - равномерную сетку с шагом

. (2.16)

Векторы, заданные выражениями (2.14) - (2.16), определяют на прямоугольной области равномерную пространственно-временную сетку:

. (2.17)

Граничные условия первого рода (Дирихле) для рассматриваемой задачи представим в виде

, (2.18)

, (2.19)

, (2.20)

, (2.21)

где , - координаты граничных точек области , ;

, - координаты граничных точек области , ;

, , , - некоторые непрерывные функции соответствующих координат.

Граничные условия второго рода (Неймана) для рассматриваемой задачи представим в виде

, (2.22)

, (2.23)

, (2.24)

. (2.25)

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

, (2.26)

, (2.27)

где - начальный момент времени;

- конечный момент времени;

, - некоторые непрерывные функции соответствующих координат.

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

, (2.28)

. (2.29)

Проводя дискретизацию граничных условий Дирихле на равномерной сетке (2.17) с использованием метода конечных разностей, получим

, (2.30)

, (2.31)

, (2.32)

, (2.33)

где ,,, - значения функции в точках , ,, соответственно.

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

, (2.34)

, (2.35)

, (2.36)

. (2.37)

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

, (2.38)

, (2.39)

где - значения функции в точке .

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

, (2.40)

. (2.41)

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

,

. (2.42)

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

Приведем один из вариантов функции с комментариями для решения системы (2.30) - (2.42) на равномерной сетке (2.17). Текст m-файла см. в приложении Б.

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

Функция может запускаться на выполнение из командной строки или из m-файла. Запуск на выполнение функции осуществляется директивой

.

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

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

Рисунок 2.3 Меню для повторного просмотра результатов

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



а

б

в

г

д

е

Рисунок 2.4 Волновой процесс u(x, y, t) в различные моменты времени:

а) t=0; б) t=0,04; в) t=0,08; г) t=0,12; д) t=0,16; е) t=0,2

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

Общий вид системы относительно неизвестной вектор-функции может быть представлен следующим уравнением:

. (2.43)

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

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

. (2.44)

Краевые условия задаются в следующем виде:

. (2.45)

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

.

Здесь параметр определяет систему координат: 0 отвечает декартовым координатам, 1 - цилиндрическим, 2 - сферическим.

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

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

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

Решение выводится в трехмерный массив , причем первый индекс () отвечает временному слою , второй индекс () определяет номер узла , а третий индекс дает номер компоненты вектора решения. Таким образом, элемент дает значение решения

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

Начальные условия даются функцией:

Краевые условия обрабатываются функцией:

Здесь векторы и дают значения вектора на разных концах интервала , - время, а выходные параметры представлены векторами на левом и правом () концах интервала.

Рассмотрим применение команды для решения уравнения теплопроводности с источником ( - параметр):

.

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

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

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

Пример реализации m-файла приведен в приложении В.

На рисунке 2.5 приведены графики, полученные в результате запуска функции .

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

Уравнение решалось численно с помощью разработанной программы при значениях параметров:

a=0;

b=3*pi/2;

n=16;

p=1;

tt=0:.2:3;

C = 1;

F = DuDx;

S = -p*exp(-t);

u0 = p+sin(x);

Pa = ua-p*exp(-t);

Qa = 0;

Pb = 0;

Qb = 1;

Рисунок 2.5 Решение уравнения теплопроводности



3. Создание графического интерфейса в среде Matlab

3.1 Графические окна системы Matlab и элементы управления

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

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

Создадим на поверхности графического окна командную кнопку:

hF1 = figure;

uicontrol ( hF1,' Style','pushbutton',...

'String', 'MyButton1',...

'Position', [ 10 10 70 30 ] );

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

В функции uicontrol, создающей элемент управления, самым важным параметром после описателя родительского окна является свойство 'Style', так как оно задает тип управляющего элемента. Задав для этого свойства значение 'pushbutton', мы создали именно кнопку.

Имена двух других свойств говорят сами за себя: String задает надпись на поверхности кнопки (в данном случае это MyButton1), a Position имеет значением вектор-строку из четырех чисел и задает положение управляющего элемента относительно левого нижнего угла графического окна. Если более конкретно, то положение левого нижнего угла кнопки относительно левого нижнего угла графического окна задают первые два элемента числовой строки. Третий же элемент этой строки задает ширину кнопки, а четвертый высоту кнопки.

3.2 Создание основных элементов управления

Для ввода и редактирования входной информации (числовой и текстовой) предназначены текстовые поля с возможностью редактирования.

Они создаются командой uicontrol (hParent, 'Style', 'edit', ... ), где многоточие означает дополнительные свойства, одно из которых, 'Position', нужно указывать практически обязательно, иначе управляющий элемент будет расположен в неудачном месте.

Создадим на поверхности графического окна текстовое поле с возможностью редактирования:

uicontrol( hF1, 'Style', 'edit',...

'Position', [ 100 10 70 30 ],...

'BackgroundColor', 'white',...

'HorizontalAlignment', 'left' ) ;

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

Заметим здесь, что свойство 'BackgroundColor' нельзя произвольно задавать для кнопок (требование операционной системы Windows).

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

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

uicontrol( hF1, 'Style', 'text',...

'Position', [ 100 50 70 30 ],...

'BackgroundColor', 'white',...

' String', `Hello, World!',...

'HorizontalAlignment', 'center' );

Как видно из приведенного фрагмента, для создания нередактируемого элемента свойство ' Style' должно иметь значение ' text'.

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

Далее создадим элемент управления, называемый флажком (checkbox):

uicontrol( hF1, 'Style', 'checkbox',...

' Position', [10 100 150 30],...

' String', 'Important option',...

'HorizontalAlignment', 'left' );

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

Далее, создадим список, содержащий несколько имен функций:

uicontrol( hF1, 'Style', 'listbox',...

'Position', [200 10 100 80],...

'String', ('sin', 'cos', 'tan'),...

'HorizontalAlignment', 'left' );

Так как список (свойство 'Style' должно быть 'listbox') содержит сразу несколько строковых значений для свойства ' String', а должно быть одно значение, то это единственное значение является массивом ячеек (каждый элемент этого массива есть текстовая строка, то есть массив типа char).

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

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

Теперь, поместим графический объект axes, предназначенный для показа графиков. Он не является графическим объектом типа uicontrol. Создается объект axes с помощью функции axes (конструктор объекта):

axes( 'Parent', hF1, 'Color', [ 1 1 1],...

'Units', 'points',...

'Position', [ 12 33 40 22 ],...

'FontSize', 6 );

где ' FontSize' размер шрифта, используемого для надписей на осях координат объекта axes.

В результате получается следующая картинка:

Рисунок 3.1 Пример оформления интерфейса

Все команды, которые создают само графическое окно, а также создают объект axes и элементы пользовательского интерфейса, удобно поместить в один m-файл. В итоге от пользователя потребуется минимальное количество действий - набрать в командном окне системы Matlab имя этого m-файла и нажать клавишу Enter.

Элемент управления togglebutton является кнопкой, отличающейся от элемента pushbutton своим поведением. Этот вариант кнопки может стабильно находиться в одном из двух состояний: нажатом или отжатом. С точки зрения внешнего вида он соответствует обычной командной кнопке. Вот код, который создает в новом графическом окне (Figure No. 2) кнопку togglebutton:

hF2 = figure;

uicontrol ( hF2,' Style',' togglebutton',...

'String', 'MyButton1',...

'Position', [ 10 10 70 30 ]);

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

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

uicontrol(hF2,'Style','radiobutton',...

'String','Third',...

'Position', [120 10 70 30]);

uicontrol(hF2,'Style','radiobutton',...

'String','Second',...

'Position', [120 50 70 30]);

uicontrol(hF2,'Style','radiobutton',...

'String','First',...

'Position', [120 90 70 30]);

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

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

uicontrol ( hF2,' Style',' frame' ,...

'Position', [ 115 5 80 120 ] );

Создадим элемент который обычно называют ползунком (по-английски - slider). Этот элемент управления имеет стиль slider:

uicontrol ( hF2,' Style',' slider',...

Position', [ 10 90 70 30 ] );

Графический элемент управления slider позволяет мышью перемещать свой ползунок, относительное положение которого задает некоторое входное число. Самое левое положение ползунка соответствует свойству Min этого элемента управления, а самое правое положение свойству Мах. Оба этих свойства нужно заранее приписать элементу управления для задания диапазона изменения входной величины. Если явно не задать числовых значений для этих свойств, то по умолчанию будут приняты величины соответственно 0 и 1.

Теперь наше иллюстративное графическое окно «Figure 2» принимает следующий вид:

Рисунок 3.2 Пример оформления интерфейса



3.3 Callback-функции

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

Отметим, что из семи созданных текстовых полей (без возможности редактирования) шесть применяются в виде поясняющих надписей, расположенных над редактируемыми текстовыми полями ввода, над списком и над объектом axes. Седьмое же текстовое поле (его описатель hTxt6) не содержит никаких надписей, так как мы не задаем явно его свойство String. Оно по умолчанию равно пустой строке.

Возникает законный вопрос, зачем же нужно такое текстовое поле, если в нем нет никакого текста? Ответ на этот вопрос заключается в том, что графический элемент управления стиля text можно использовать для нанесения на поверхность графического окна системы Matlab различных элементов оформления, например разграничительных линий. Именно такую вертикальную линию мы и формируем, создавая текстовое поле с помощью кода

hTxt6 = uicontrol ( hFigl,'Style',' text',...

' BackgroundColor', [0 0 0], ...

'Position', [430 5 1 900] );

в котором задаем цвет фона управляющего элемента отличным от цвета фона графического окна (цвет [0 0 0] это черный цвет, а цвет графического окна по умолчанию серый). Кроме того, мы задаем ширину этого элемента в один пиксел.

В результате получается графическое окно, изображенное на рисунке 3.3:

Рисунок 3.3 Пример оформления интерфейса

Это окно предназначено для численного решения дифференциального уравнения

;

;

;

где функция выбирается из списка, расположенного в правом нижнем углу графического окна. Остальные параметры задаются с клавиатуры в редактируемых полях. Кнопка GO предназначена для начала расчетов по методу Эйлера и показа результата в виде графика на поверхности объекта axes. Кнопка CLEAR позволяет очистить объект axes от содержимого.

Как мы видим из рисунка, предусмотрены расчеты со стандартными функциями sin(x) и cos(x), а также с собственными функциями , и . Последние должны быть записаны в m-файлы MyFunct1 .m, MyFunct2 .m, f3.m.

Число N вводится перед расчетами с клавиатуры и определяет шаг интегрирования, который равен

.

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

Функции, которые связываются с графическими элементами управления, называют callback-функциями (их не вызывают наши m-функции, они вызываются средой Matlab у нас «за спиной», отсюда и происхождение их названия).

Чтобы связать наши командные кнопки с callback-функциями, которые будут вызываться средой Matlab при «нажатии» на эти кнопки, при их создании функцией uicontrol нужно прописать свойство ' Callback', указав ему в качестве значения имя соответствующей m-функции. Эти callback-функции можно также назвать функциями обработки событий, связанных с элементами управления. Это значит, что в представленном m-файле MyFig1 .m необходимо добавить свойства: 'Callback', 'MyGO' - для hBut1 и 'Callback', 'MyCLEAR' - для hBut2.

Внеся указанные выше изменения в файл MyFig1.m, перейдем к написанию функций MyGO и MyCLEAR. Начнем с более простой функции MyCLEAR. Вот текст этой функции:

function MyClear

global hAxes

axes (hAxes);

cla

Ранее мы уже в файле MyFig1 .m обозначили все описатели созданных нами графических объектов как глобальные. Таким образом, у нас имеется очень простой к ним доступ. Все функции, где эти описатели маркированы как глобальные, имеют к ним непосредственный доступ. Именно так все и происходит в функции MyCLEAR (ее мы записываем в файл MyCLEAR.m).

Поясним текст этой функции. Здесь команда cla стирает содержимое текущего (активного) объекта axes. Предварительный вызов функции axes (hAxes) позволяет гарантировать, что будет стерто содержимое именно объекта axes с описателем hAxes, так как такой вызов делает активным объект, описатель которого указан в качестве параметра. У нас в данной конфигурации имеется всего один объект axes, так что этот вызов не является обязательным, но он иллюстрирует решение проблемы в случае нескольких объектов axes.

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

function MyGo

global hAxes hEd1 hEd2 hEd3 hEd4 hList

str1=get(hEd1,'String');

str2=get(hEd2,'String');

str3=get(hEd3,'String');

str4=get(hEd4,'String');

x0 = str2num(str1);

y0 = str2num(str2);

xf = str2num(str3);

N = str2num(str4);

index = get(hList,'Value');

cellArr = get(hList,'String');

funName = cellArr{index};

dx= (xf-x0)/N;

X=x0;

Y=y0;

for k=1:N

Y=[Y,Y(end)+feval(funName,X(end))*dx];

X=[X,X(end)+dx];

end

axes(hAxes);

plot(X,Y);

С помощью функций get, указывая им в качестве первого аргумента описатели редактирующих полей, мы получаем значения свойства 'String', то есть текстовое содержимое этих полей, которое вводится пользователем с клавиатуры. Это и есть входные данные. Только пока что числовые данные представлены в строковой (текстовой) форме.

Для получения имени функции, представляющей правую часть дифференциального уравнения, сначала у списка запрашиваем через свойство 'Value' индекс подсвеченной (выбранной) строки, затем читаем значение свойства 'String'. Последнее для списка является массивом ячеек, откуда и выбираем имя функции с помощью операции индексирования массива ячеек фигурными скобками.

С помощью функций str2num преобразовываем текстовые величины в числа. Затем устанавливаем начальные значения для массива значений аргумента и массива значений функции на сетке. Шаг сетки (шаг интегрирования) вычисляется делением длины отрезка интегрирования на количество шагов , которое пользователь вводит с клавиатуры. Наращивание массивов осуществляется операцией конкатенации, в то время как очередное значение функции вычисляется согласно алгоритму метода Эйлера и равно

feval( funName, X(end) ) * dx.

Ключевое слово системы Matlab end в индексирующих выражениях означает ссылку на последний элемент массива, что нам и требуется в данном случае. С помощью функции feval осуществляется вызов функции с именем funName, прочитанным из списка, и ей передается аргумент X(end) (текущее значение независимого аргумента).

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

Проверим работу приложения для пользовательских функций MyFunct1, MyFunct2 и f3. Текст функции MyFunct1:

function ret = MyFunct1(x)

ret=x^3;

записываем в файл MyFunct1 .m.

Текст функции MyFunct2

function ret = MyFunct2(x)

ret=1/sin(x);

записываем в файл MyFunct2 .m. Текст функции f3

function ret = f3( х ) ret = sin( x * x ) ;

записываем в файл f3.m.

Получим решение для данных функций.

sin

cos

MyFunct1

MyFunct2

f3

Рисунок 3.4 Пример численного решения дифференциального уравнения

4. Решение прикладных задач средствами Matlab

4. 1 Создание интерфейса

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

. (4.1)

Начальное распределение температуры есть

. (4.2)

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

В результате получается графическое окно, изображенное на рисунке 4.1:

Рисунок 4.1 Оформления интерфейса

Кнопка «Решение» предназначена для начала расчетов и показа результата в виде графиков на поверхности объектов axes. Кнопка «Очистить графики» позволяет очистить объекты axes от содержимого. Наши командные кнопки при помощи callback-функций связаны с m-функциями `function_pdepe' и `MyClear'.

Текст m-функции `MyClear':

function MyClear

global hAxes1 hAxes2

axes(hAxes1);

cla

axes(hAxes2);

cla

Текст m-функции `function_pdepe' см. приложение Е.

4.2 Зависимость решений уравнения теплопроводности от заданных параметров

...