Решение Дифференциального уравнения в среде MATLAB

Введение

MATLAB как язык программирования был разработан Кливом Моулером (англ. Cleve Moler) в конце 1970-х годов, когда он был деканом факультета компьютерных наук в Университете Нью-Мексико. Целью разработки служила задача дать студентам факультета возможность использования программных библиотек Linpack и EISPACK без необходимости изучения Фортрана. Вскоре новый язык распространился среди других университетов и был с большим интересом встречен учёными, работающими в области прикладной математики. До сих пор в Интернете можно найти версию 1982 года, написанную на Фортране, распространяемую с открытым исходным кодом. Инженер Джон Литтл (англ. John N. (Jack) Little) познакомился с этим языком во время визита Клива Моулера в Стэнфордский университет в 1983 году. Поняв, что новый язык обладает большим коммерческим потенциалом, он объединился с Кливом Моулером и Стивом Бангертом (англ. Steve Bangert). Совместными усилиями они переписали MATLAB на C и основали в 1984 компанию The MathWorks для дальнейшего развития. Эти переписанные на С библиотеки долгое время были известны под именем JACKPAC.

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

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

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

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

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

Математическая статистика и анализ данных — статистические функции, статистическая регрессия, цифровая фильтрация, быстрое преобразование Фурье и другие.

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

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

Разреженные матрицы — специальный класс данных пакета MATLAB, использующийся в специализированных приложениях.

Целочисленная арифметика — выполнение операций целочисленной арифметики в среде MATLAB.

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

Цифровая обработка сигналов, изображений и данных: DSP Toolbox, Image Processing Toolbox, Wavelet Toolbox, Communication Toolbox, Filter Design Toolbox — наборы функций, позволяющих решать широкий спектр задач обработки сигналов, изображений, проектирования цифровых фильтров и систем связи.

Системы управления: Control Systems Toolbox, µ-Analysis and Synthesis Toolbox, Robust Control Toolbox, System Identification Toolbox, LMI Control Toolbox, Model Predictive Control Toolbox, Model-Based Calibration Toolbox — наборы функций, облегчающих анализ и синтез динамических систем, проектирование, моделирование и идентификацию систем управления, включая современные алгоритмы управления, такие как робастное управление, H∞-управление, ЛМН-синтез, µ-синтез и другие.

Финансовый анализ: GARCH Toolbox, Fixed-Income Toolbox, Financial Time Series Toolbox, Financial Derivatives Toolbox, Financial Toolbox, Datafeed Toolbox — наборы функций, позволяющие быстро и эффективно собирать, обрабатывать и передавать различную финансовую информацию.

Анализ и синтез географических карт, включая трёхмерные: Mapping Toolbox.

Сбор и анализ экспериментальных данных: Data Acquisition Toolbox, Image Acquisition Toolbox, Instrument Control Toolbox, Link for Code Composer Studio — наборы функций, позволяющих сохранять и обрабатывать данные, полученные в ходе экспериментов, в том числе в реальном времени. Поддерживается широкий спектр научного и инженерного измерительного оборудования.

Визуализация и представление данных: Virtual Reality Toolbox — позволяет создавать интерактивные миры и визуализировать научную информацию с помощью технологий виртуальной реальности и языка VRML.

Средства разработки: MATLAB Builder for COM, MATLAB Builder for Excel, MATLAB Builder for NET, MATLAB Compiler, Filter Design HDL Coder — наборы функций, позволяющих создавать независимые приложения из среды MATLAB.

Взаимодействие с внешними программными продуктами: MATLAB Report Generator, Excel Link, Database Toolbox, MATLAB Web Server, Link for ModelSim — наборы функций, позволяющие сохранять данные в различных видов таким образом, чтобы другие программы могли с ними работать.

Базы данных: Database Toolbox — инструменты работы с базами данных.

Научные и математические пакеты: Bioinformatics Toolbox, Curve Fitting Toolbox, Fixed-Point Toolbox, Fuzzy Logic Toolbox, Genetic Algorithm and Direct Search Toolbox, OPC Toolbox, Optimization Toolbox, Partial Differential Equation Toolbox, Spline Toolbox, Statistic Toolbox, RF Toolbox — наборы специализированных математических функций, позволяющие решать широкий спектр научных и инженерных задач, включая разработку генетических алгоритмов, решения задач в частных производных, целочисленные проблемы, оптимизацию систем и другие.

Нейронные сети: Neural Network Toolbox — инструменты для синтеза и анализ нейронных сетей.

Нечёткая логика: Fuzzy Logic Toolbox — инструменты для построения и анализа нечётких множеств.

Символьные вычисления: Symbolic Math Toolbox — инструменты для символьных вычислений с возможностью взаимодействия с символьным процессором программы Maple.

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

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

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

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

Пакет MATLAB включает различные интерфейсы для получения доступа к внешним подпрограммам, написанным на других языках программирования, данным, клиентам и серверам, общающимся через технологии Component Object Model или Dynamic Data Exchange, а также периферийным устройствам, которые взаимодействуют напрямую с MATLAB. Многие из этих возможностей известны под названием MATLAB API.

Раздел 1

В данной курсовой работе попробуем решить вариационную задачу

(4.9)

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

clear all

format long

disp('Решаем пример 4')

syms x y z Dzx Dzy D2zx2 D2zxy D2zy2 % переменные

F = Dzx^2-2*Dzy^2+2*y*z*(sin(pi*x)+x/5); % подынтегральная функция

zc = x/10+y^2/50; % функция граничных условий

x1=0; % координаты ограничивающего прямоугольника

x2=1;

y1=0;

y2=2;

fprintf('Подынтегральная функция: F=%s\n',char(F))

fprintf('Граничное условие на контуре: z=%s\n',char(zc))

fprintf('Область: %d<=x<=%d; %d<=y<=%d\n',x1,x2,y1,y2)

Решаем пример 4

Подынтегральная функция: F=Dzx^2-2*Dzy^2+2*y*z*(sin(pi*x)+1/5*x)

Граничное условие на контуре: z=1/10*x+1/50*y^2

Область: 0<=x<=1; 0<=y<=2

Находим частные производные Fz, Fp и Fq.

dFdz = diff(F,z)

dFdp = diff(F,Dzx)

dFdq = diff(F,Dzy)

dFdz =

2*y*(sin(pi*x)+1/5*x)

dFdp =

2*Dzx

dFdq =

-4*Dzy

Сформируем полные частные производные  Fp/ x и  Fq/ y. При их формировании учитываем, что z=z(x,y), p=p(x,y), q=q(x,y). Используем формулу (1.7).

d_dFdp_dx = diff(dFdp,x);

d_dFdp_dz = diff(dFdp,z);

d_dFdp_dp = diff(dFdp,Dzx);

d_dFdp_dq = diff(dFdp,Dzy);

dFpdx = d_dFdp_dx + d_dFdp_dz*Dzx + d_dFdp_dp*D2zx2 + d_dFdp_dq*D2zxy

d_dFdq_dy = diff(dFdq,y);

d_dFdq_dz = diff(dFdq,z);

d_dFdq_dp = diff(dFdq,Dzx);

d_dFdq_dq = diff(dFdq,Dzy);

dFqdy = d_dFdq_dy + d_dFdq_dz*Dzy + d_dFdq_dp*D2zxy + d_dFdq_dq*D2zy2

dFpdx =

2*D2zx2

dFqdy =

-4*D2zy2

Формируем уравнение Эйлера-Остроградского.

Euler = simple(dFdz-dFpdx-dFqdy)

EuR = -subs(Euler,{z,D2zx2,D2zy2,D2zxy},{0,0,0,0}); % правая часть

EuL = Euler+EuR; % левая часть

deqEuler = [ char(EuL) '=' char(EuR) ]; % уравнение

fprintf('Уравнение Эйлера-Остроградского:\n%s\n',deqEuler)

Euler =

2*y*sin(pi*x)+2/5*y*x-2*D2zx2+4*D2zy2

Уравнение Эйлера-Остроградского:

-2*D2zx2+4*D2zy2=-2*y*sin(pi*x)-2/5*y*x

Для решения дифференциальных уравнений в частных производных в MATLAB'е есть специальный инструментарий – Partial Differential Equation Toolbox (PDE), в котором используется метод конечных элементов (FEM). Для применения PDE нужно привести дифференциальное уравнение к виду

(4.10)

и дополнить его граничными условиями типа Дирихле

(4.11)

или Неймана

. (4.12)

Здесь u(x,y) – искомая функция, C – матрица 2 x 2, элементы которой являются коэффициентами при
 2u/ x2,  2u/ x y,  2u/ y2; a – коэффициент при u; f – правая часть, h, r, q, g – заданные функции x, y. Величины C, a, f, h, r, q, g могут быть как постоянными, так и переменными. В последнем случае они должны вычисляться в центрах тяжести конечных элементов и задаваться как массивы.

Раздел 2.

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

                    задание геометрии (области решения) и построение FEM-сетки;

                    задание граничных условий (функций h, r, q, g);

                    задание функций, входящих в дифференциальное уравнение (c, a, f);

                    решение дифференциального уравнения;

                    отображение результатов в виде графика.

Задать геометрию в виде прямоугольника проще всего с помощью команды poimesh. Она формирует сетку на квадрате [–1,1] x [–1,1]. Эта команда возвращает 3 выходных параметра: p, e и t. В переменной p возвращаются координаты узлов сформированной сетки: 1-я строка – x-е, 2-я – y-е. Массив p имеет размеры 2 x np, где np – число узлов. В переменной t возвращаются данные по треугольникам. Это массив размером 4 x nel, где nel – число элементов. Первые 3 строки содержат номера узлов для каждого элемента в порядке обхода против часовой стрелки. 4-е число в нашем случае – всегда 1. В переменной e возвращаются данные по граничным точкам сетки. Размер массива e 7 x nb, где nb – число участков границы в FEM-сетке. В каждом столбце массива e содержатся данные по одной граничной линии. 1-е и 2-е числа – это номера узлов (в том порядке, в котором они перечислены в массиве p).

Сформируем треугольную FEM – сетку и нарисуем её. Напечатаем количество узлов и элементов.

del=min(x2-x1,y2-y1)/20 % размер элемента

nx=round((x2-x1)/del)

ny=round((y2-y1)/del)

[p,e,t]=poimesh('squareg',nx,ny); % формируем FEM-сетку

p(1,:)=(p(1,:)+1)/2*(x2-x1)+x1;

p(2,:)=(p(2,:)+1)/2*(y2-y1)+y1; % приводим к прямоугольной

np=size(p,2); % число узлов

nel=size(t,2); % число элементов

fprintf('Число узлов сетки разбиения np=%d\n',np)

fprintf('Число конечных элементов nel=%d\n',nel)

pdemesh(p,e,t) % рисуем сетку

da=daspect; % текущие масштабы осей

da(1:2)=min(da(1:2)); % выравняли масштабы

daspect(da); % установили одинаковые масштабы

title('\bfFEM Grid') % заголовок

xlabel('x')

ylabel('y') % метки осей

del =

0.05000000000000

nx =

20

ny =

40

Число узлов сетки разбиения np=861

Число конечных элементов nel=1600

Рис. 4.1. Разбивка области на конечные элементы

Следующий этап – это задание граничных условий. Граничные условия можно задать или в виде матрицы граничных условий (boundary condition matrix), или в виде *.m-файла граничных условий (boundary M-file). Второй вариант проще. Файл граничных условий должен иметь такую структуру:

[q,g,h,r]=pdebound(p,e,u,time)

Входные параметры: p,e – данные по сетке разбиения, u – решение, time – время. Выходные параметры – матрицы граничных условий Дирихле (4.11) или Неймана (4.12). PDE позволянт решать параболические и гиперболические уравнения (зависящие от времени), а также нелинейные задачи, поэтому граничные условия и параметры дифференциального уравнения могут зависеть также от времени и решения.

Для граничных условий Неймана матрицы q и g должны содержать значения параметров q и g в средних точках границ. Размер матрицы q: N2 x nb, где N – число уравнений системы, а nb – число сторон конечных элементов, выходящих на границу. Размер матрицы g: N x nb. В каждом из столбцов этих матриц должны возвращаться коэффициенты q и g в средних точках соответствующей границы. Для нескольких уравнений элементы q должны следовать по столбцам. Так, например, для 2-х уравнений в каждом столбце матрицы q нужно возвратить
q11, q21, q12, q22. В случае граничных условий Дирихле в этих матрицах должны возвращаться нулевые значения.

Для граничных условий Дирихле должны формироваться матрицы h и r. Матрица h имеет размеры N2 x (2nb), и содержит значения h сперва во всех начальных точках каждой границы, и сразу за ними – в конечных точках границ. Размер матрицы r: N x (2nb), этот масси заполняется аналогично.

Сформируем файл для вычисления граничных условий Дирихле. Число уравнений у нас N=1, число граничных точек
nb находим из массива e. Формируем строки для записи в файл и записываем их. Файл должен быть размещён в каком-либо каталоге, доступном системе MATLAB.

s{1}='function [q,g,h,r]=boundmem(p,e,u,time)'; % заголовок

s{2}='nb=size(e,2);'; % определили размерности

s{3}='q=zeros(1,nb); g=zeros(1,nb);'; % задали массивы условий Неймана