39076

MATLAB. Численное решение дифференциальных уравнений

Лабораторная работа

Информатика, кибернетика и программирование

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

Русский

2015-01-11

315.39 KB

426 чел.

Лабораторная работа №6

MATLAB. Численное решение дифференциальных уравнений.

Исходные данные: математический пакет MATLAB.

Задание:

  1.  Генерация случайных чисел
  2.  Решение дифференциального уравнения с начальными условиями (первого  и второго порядков)
  3.  Решение уравнений в частных производных

Теоретическая справка:

Генерация случайных чисел в MatLab осуществляет функция randn(n,m)

x = -4:0.1:4;

y = randn(10000,1);

% График  нормального распределения

Решение обыкновенных дифференциальных уравнений

Для решения обыкновенных дифференциальных уравнений (ODE) могут быть применены численные методы, которые в MATLAB реализованы в специальных функциях-решателях: ode45, ode23, ode113.

Функции ode23 и ode45 предназначены для численного интегрирования систем ОДУ. Они применимы как для решения простых дифференциальных уравнений, так и для моделирования сложных динамических систем.

Любая система нелинейных ОДУ может быть представлена как система дифференциальных уравнений 1-го порядка в явной форме Коши:

Синтаксис:

           [t, X] = ode23(‘<имя функции>‘, t0, tf, x0)

            [t, X] = ode23(‘<имя функции>‘, t0, tf, x0, tol, trace)

            [t, X] = ode45(‘<имя функции>‘, t0, tf, x0)

            [t, X] = ode45(‘<имя функции>‘, t0, tf, x0, tol, trace)

где x - вектор состояния;

      t - время;

      f - нелинейная вектор-функция от переменных x, t.

Функции

[t, X] = ode23(‘<имя функции>‘, t0, tf, x0, tol, trace) и

[t, X] = = ode45(‘<имя функции>‘, t0, tf, x0, tol, trace) интегрируют системы ОДУ, используя формулы Рунге - Кутта соответственно 2-го и 3-го или 4-го и 5-го порядка.

Эти функции имеют следующие параметры:

Входные параметры:

‘<имя функции>‘ - строковая переменная, являющаяся именем М-файла, в котором вычисляются правые части системы ОДУ;

 t0 - начальное значение времени; tfinal - конечное значение времени;

 x0 - вектор начальных условий;

 tol - задаваемая точность; по умолчанию для ode23 tol = 1.e-3, для ode45 tol = 1.e-6);

 trace - флаг, регулирующий вывод промежуточных результатов; по умолчанию равен нулю, что подавляет вывод промежуточных результатов;

Выходные параметры:

 t - текущее время;

 X - двумерный массив, где каждый столбец соответствует одной переменной.

Алгоритм:

Функции ode23 и ode45 реализуют методы Рунге - Кутты с автоматическим выбором шага, описанные в работе [2]. Такие алгоритмы используют тем большее количество шагов, чем медленнее изменяется функция. Поскольку функция ode45 использует формулы более высокого порядка, обычно требуется меньше шагов интегрирования и результат достигается быстрее. Для сглаживания полученных процессов можно использовать функцию interp1. 

Установление заданной относительной погрешности  RelTolODESET(odeset).

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

Пример Уравнение Ван-дер-Поля с заданной относительной погрешностью.

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

% Формат записи функции odeset включает строковый служебный символ RelTol и задаваемую %величину относительной погрешности (0.1 для d1 и 0.2 для d2).

Рассмотрим пример решения дифференциального уравнения с начальными условиями.

Решить задачу Коши

Точное решение имеет вид

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

 Funcrion dydx=F(x,y) 

Dydx=zeros(1,1);

Dydx(1)=2*(x^2+y(1));

Решение будет таким

>>[X Y]=ode45(@F,[01],[1]);

>>plot(X,Y);

Hold on; gtext(‘y(x)’)

>>[X Y]

 

  Формирование прямоугольной сетки на плоскости — meshgrid.

» [x,y]=meshgrid(-5:0.5:5,-5:0.5:5);

» plot(x,y),xlabel('X'),ylabel('Y')

%  Результатом действия функции meshgrid является формирование "основания" в плоскости XOY для построения над этим основанием  пространственной фигуры.

Построение пространственных сетчатых фигур — mesh.

Пример

[x,y]=meshgrid(-5:0.1:5,-5:0.1:5);

Z = 1*x.* exp(-x.^2 - y.^2);% Коэффициент 1 заменить: 2, 5, 10, 20

 mesh(Z),xlabel('X'),ylabel('Y'),zlabel('Z'),title('Z-поверхность')

% Для сравнения применить  plot3(x,y,Z),grid вместо  mesh(Z).

Пример

[x,y]=meshgrid(-15:0.2:15,-15:0.2:15);

 R=sqrt(x.^2+y.^2)+0.001; %Коэфф. 0.001 введен для исключения деления на ноль

Z=1*sin(R)./R;

mesh(Z)

% Коэффициент 1 заменить: 5, 10, 20, 0.5, 0.1

% Для сравнения применить  plot3(x,y,Z),grid

Сетчатая поверхность с проекциями линий постоянного уровня meshс.

Пример

 [x,y]=meshgrid(-15:0.2:15,-15:0.2:15);

R=sqrt(x.^2+y.^2)+0.001;

Z=1*sin(R)./R;

meshc(Z)

% Коэффициент 1 заменить: 5, 10, 20, 0.5, 0.1

Пример 4.2.

[x,y]=meshgrid(-5:0.1:5,-5:0.1:5);

Z = 1*x.* exp(-x.^2 - y.^2);

meshc(Z),xlabel('X'),ylabel('Y'),zlabel('Z'),title('Z-поверхность')

% Коэффициент 1 заменить: 5, 10, 20, 0.5, 0.1

Сетчатая поверхность с пьедесталом плоскости отсчета на нулевом уровне— meshz.

Пример 5.1.

[x,y]=meshgrid(-15:0.2:15,-15:0.2:15);

R=sqrt(x.^2+y.^2)+0.001;

Z=1*sin(R)./R;

meshz(Z)

% Коэффициент 1 заменить: 5, 10, 20, 0.5, 0.1

Пример

[x,y]=meshgrid(-5:0.1:5,-5:0.1:5);

Z = 1*x.* exp(-x.^2 - y.^2);

 meshz(Z),xlabel('X'),ylabel('Y'),zlabel('Z'),title('Z-поверхность')

% Коэффициент 1 заменить: 5, 10, 20, 0.5, 0.1

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

Эта модель Вольтерра-Лотка с логистической поправкой описывается системой уравнений

С условиями заданной численности "жертв" и "хищников" в начальный момент t=0.

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

 function f=VolterraLog(t,x)

 a=4;

 b=2.5;

 c=2;

 d=1;

 alpha=0.1;

 f(1)= (a-b*x(2))*x(1)-alpha*x(1)^2;

 f(2)= (-c+d*x(1))*x(2)-alpha*x(2)^2;

 f=f';

>>opt=odeset('OutputSel',[1 2], 'OutputFcn', 'odephas2');

>> [T,X]=ode45('VolterraLog', [0 10],[3 1],opt );

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

=x2x3,

=-x1x3,

=-0.51Чx1x2

выступает в виде:

 function f=Eiler(t,x)

  f(1)=  x(2)*x(3);

  f(2)=  -x(1)*x(3) ;

  f(3)=  -0.51*x(1)*x(2) ;

  f=f';>> opt=odeset('OutputSel',[1 2 3], 'OutputFcn', 'odephas3');

>> [T,X]=ode45('Eiler', [0 7.25], [0 0 1], opt ); %рис.8.10

>> [T,X]=ode45('Eiler', [0:0.25:7.25],[0 1 1]);

>> plot(T,X)

Еще один пример применения функций:

. function dn = population(t,n)

dn = zeros(2,1);  

a = 0.1; f = 0.45; r = 4; q = 1.5;

dn(1) = r * n(1) - a * n(2) * n(1);

dn(2) = f * a * n(2) * n(1) - q * n(2);

x = (-5:0.05:5);

Y = pdf('normal',x,0,1);

plot(x,Y)

pause

n0 = 90; c0 = 5;

opt=odeset('OutputSel',[1 2], 'OutputFcn','odephas2');

[T,Y] = ode45(@lab4_f,[0 10],[n0 c0],opt);

pause

plot(T,Y(:,1),T,Y(:,2))

Initmesh

Синтаксис:

[p,e,t]=initmesh(g,p1,v1,p2,v2,p3,v3,p4,v4,p5,v5,p6,v6,p7,v7,p8,v8)

Описание:

Возвращает треугольную конечноэлементную сетку, построенную в расчётной области, геометрия которой описана в m-функции g.

Обязательным является только первый входной параметр g.

Входные параметры.

p1,v1,… - список необязательных ключевых («именных») параметров функции initmesh:

 pk - строки символов - имена указываемых параметров;

 vk - значения указываемых параметров.

Имена ключевых параметров, их назначение и допустимые значения представлены в таблице:pk Значение vk / {По умолчанию} Описание

Hmax Числовое значение {estimate} Максимальный размер треугольников

Hgrad Числовое значение {1.3} Показатель нерегулярности треугольников

Box on|{off} Генерация сетки в пределах ограниченного прямоугольника

Init on|{off} Выполнить только начальную триангуляцию границ

Jiggle off|{mean}|min Вызов функции jigglemesh

JiggleIter Числовое значение {10} Максимальное число итераций при регуляризации сетки

Параметры Box и Init связаны с работой алгоритма построения сетки. Их полезно устанавливать в ‘on’ для изучения алгоритма построения

сетки. Параметры Jiggle и JiggleIter используются для управления уровнем регуляризации конечноэлементной сетки (подробнее см. jigglemesh).

Выходные параметры.

p - массив узлов конечноэлементной сетки (столбцам соответствуют узлы):

- первая строка - горизонтальные координаты узлов,

- вторая строка -вертикальные координаты узлов;

e - матрица граничных элементов на границах раздела зон (см. pdegeom):

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

первые две строки - номера номера начальных и конечных узлов граничных элементов;

строки 3, 4 - длина «дуги» от начала граничного сегмента до начального и конечного узла граничного элемента, отнесённая к длине «дуги» граничного сегмента;

строка 5 - номера граничных сегментов, которым принадлежат граничные элементы;

строки 6, 7 - номера зон, примыкающих слева и справа к граничным элементам;

t - матрица треугольных конечных элементов (столбцам соответствуют треугольники):

- t(1:3,ie) - глобальные номера узлов треугольника с номером ie,

- t(4,ie) - номер зоны, которой принадлежит треугольник с номером ie.

Пакет Partial Differential Equations Toolbox 

Специалистов в области численных методов решения дифференциальных уравнений в частных производных несомненно заинтересует обширный пакет Partial Differential Equations Toolbox (PDETB). Хотя этот пакет является самостоятельным приложением и в ядро MATLAB не входит, мы приведем краткое описание некоторых его возможностей с парой примеров. Вы можете вызвать пакет с его графическим интерфейсом командой pdetool.

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

[p,e,t]=initmesh('Ishapeg');  

u=assempde('Ishapeb'p.e.t.1.0,1):  

pdeplottp.e.t.'xydata',u.'zdata'.u.'mesh'.'off'):

В состав пакета входит ряд полезных демонстрационных примеров с именами от pdedemol до pdedemo8. Их можно запустить как из командной строки (путем указания имени), так из окна демонстрационных примеров Demos. Рекомендуется, однако, сделать это из командной строки, так как примеры ориентированы на управление в командном режиме — в основном оно сводится к нажатию клавиши Enter для прерывания пауз, введенных для просмотра промежуточных результатов вычислений.

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

 -div( l/sqrt(l+grad|u|  ^  2) * grad(u) ) = 0

при граничном условии u=х*2. Ниже представлен текст файла pdedemo3.m с убранными англоязычными комментариями:

%PDEDEM03 Решение проблемы минимизации параболической поверхности

 echo on 

 clc 

 g='circleg':  %  Единичная окружность

 b='circleb2':  %  х  ^  2 - граничное условие

 c='l./sqrt(l+ux.*2+uy.*2)':

а=0;  

 f=0:

 rto1=le-3; % Погрешность вычислений решателем

 pause  %  Пауза в вычислениях

 clc 

%  Генерация поверхности

 [p,e,t]=initmesh(g); [p,e.t]-refinemesh(g.p.e.t);

%  Решение нелинейной проблемы

u=pdenonlin(b.p.e,t.c.a.f,'tol'.rtol);

pdesurf(p.t.u);

pause  %  Пауза в вычислениях

 echo off 

Весьма интересны и поучительны примеры с анимацией:

pdedemo2 — появление и распространение волн,

pdedemoS — вздутие поверхности (пузырек газа),

pdedemo6 — колебания плоскости (гиперболическая проблема) и т. д.

Переопределение (сгущение) треугольной сетки с помощью REFINEMESH

Синтаксис.

[p1,e1,t1]=refinemesh(g,p,e,t)

Описание.

Возвращает переопределённую версию треугольной конечноэлементной сетки, представленной геометрией g, матрицей узлов p, матрицей граничных элементов e и матрицей треугольников t.

Входные параметры.

g описывает геометрию PDE задачи. g может матрицей “расчленённой” геометрии или именем m-файла описания геометрии. См также decsg, pdegeom.

p, e, t - матрицы описания конечноэлементной сетки (см initmesh).

Синтаксис.

[p1,e1,t1,u1]=refinemesh(g,p,e,t,u)

Описание.

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

Входные параметры.

u - узловое распределение искомой функции (величины) на непереопределённой сетке. Строкам u и столбцам p соответствуют узлы. Строкам u1 и столбцам p1 соответствуют узлы переопределённой сетки. Каждый столбец u интерполируется отдельно.

Синтаксис.

[p1,e1,t1]=refinemesh(g,p,e,t,u,it) или [p1,e1,t1,u1]=refinemesh(g,p,e,t,u,it)

Входные параметры.

it - список зон для переопределения сетки, если it - матрица-строка, или список треугольников, если it - матрица-столбец.

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

Синтаксис.

[p1,e1,t1]=refinemesh(g,p,e,t,u,it,mode) или [p1,e1,t1,u1]=refinemesh(g,p,e,t,u,it,mode)

Входные параметры.

mode - метод переопределения сетки (по умолчанию - “регулярное переопределение”, когда каждый переопределяемый треугольник разделяется на четыре треугольника).

Этот параметр может принимать одно из следующих значений:

- longest - треугольники делятся на два так, что большая сторона делится пополам;

- regular - “регулярное переопределение”.

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

 Формирование линейного массива равноотстоящих узлов LINSPACE

Синтаксис:

            x = linspace(x1, x2)

             x = linspace(x1, x2, n)

Описание:

Функция x = linspace(x1, x2) формирует линейный массив размера 1 х 100, начальным и конечным элементами которого являются точки x1 и x2.

Функция x = linspace(x1, x2, n) формирует линейный массив размера 1 х n, начальным и конечным элементами которого являются точки x1 и x2.

 Решение параболической PDE задачи с помощью PARABOLIC

Синтаксис.

u1=parabolic(u0,tlist,b,p,e,t,c,a,f,d)

Возможны также следующие варианты вызова данной функции:

 u1=parabolic(u0,tlist,b,p,e,t,c,a,f,d,rtol),

u1=parabolic(u0,tlist,b,p,e,t,c,a,f,d,rtol,atol).

 Здесь rtol, atol - относительная и абсолютная погрешность решателя ODE.

Описание:

Производит решение скалярной PDE задачи, основанной на уравнении вида

d*du/dt-div(c*grad(u))+a*u=f  (1)

Входные параметры.

u0 - узловое распределение искомой величины в начальный момент времени t=0;

tlist – список моментов времени, для которых нужно вычислить решение уравнения (1);

b - имя пользовательской m-функции, вычисляющей матрицы описания граничных условий (см. pdebound);

p - массив узлов конечноэлементной сетки (столбцам соответствуют узлы):

- первая строка - горизонтальные координаты узлов,

- вторая строка - вертикальные координаты узлов;

e - матрица граничных элементов на границах раздела зон (см. initmesh, pdegeom);

t - матрица треугольных конечных элементов (столбцам соответствуют треугольники):

- t(1:3,ie) - глобальные номера узлов треугольника с номером ie,

- t(4,ie) - номер зоны, которой принадлежит треугольник с номером ie.

с - массив, описывающий распределение коэффициента c в расчётной области (см. уравнение (1));

a - массив, описывающий распределение коэффициента a в расчётной области (см. уравнение (1));

f - массив, описывающий распределение правой части PDE f в расчётной области (см. уравнение (1));

d - массив, описывающий распределение коэффициента d в расчётной области (см. уравнение (1)).

Кодирование входных параметров c, a, f, d более подробно описано в assempde. 

Выходной параметр.

В случае скалярной PDE задачи u1 - матрица размера (NP,length(tlist)), где NP - число узлов конечноэлементной сетки. Каждый столбец матрицы u1 представляет собой узловое распределение искомой величины u в соответствующий момент времени. В случае системы PDE из N уравнений u1 - матрица размера (NP*N,length(tlist)). Каждый столбец матрицы u1 состоит из N подстолбцов, каждый из которых представляет собой узловое распределение соответствующей искомой переменной в соответствующий момент времени.

Пример решения параболического уравнения.

%        du/dt-div(grad(u))=0

[p,e,t]=initmesh('squareg');

[p,e,t]=refinemesh('squareg',p,e,t);

u0=zeros(size(p,2),1);

ix=find(sqrt(p(1,:).^2+p(2,:).^2)<0.4);

u0(ix)=ones(size(ix));

tlist=linspace(0,0.1,20);

u1=parabolic(u0,tlist,'squareb1',p,e,t,1,0,1,1);


 

А также другие работы, которые могут Вас заинтересовать

7188. ОПРЕДЕЛЕНИЕ МОМЕНТА ИНЕРЦИИ КОЛЕСА МЕТОДОМ КОЛЕБАНИЙ 151 KB
  ОПРЕДЕЛЕНИЕ МОМЕНТА ИНЕРЦИИ КОЛЕСА МЕТОДОМ КОЛЕБАНИЙ Цель работы: Определение характеристик колебательного движения колеса, момента инерции колеса и сравнение его с теоретическим расчётом. Оборудование: экспериментальная установка,...
7189. ОПРЕДЕЛЕНИЕ ВЯЗКОСТИ ВОЗДУХА, СРЕДНЕЙ ДЛИНЫ СВОБОДНОГО ПРОБЕГА МОЛЕКУЛ И ИХ ЭФФЕКТИВНОГО ДИАМЕТРА 226.5 KB
  ОПРЕДЕЛЕНИЕ ВЯЗКОСТИ ВОЗДУХА, СРЕДНЕЙ ДЛИНЫ СВОБОДНОГО ПРОБЕГА МОЛЕКУЛ И ИХ ЭФФЕКТИВНОГО ДИАМЕТРА Цель работы: Определение вязкости воздуха, средней длины свободного пробега молекул и их эффективного диаметра с использованием легко измеряемых ...
7190. Поняття та особливості валютного контролю 559 KB
  Поняття та особливості валютного контролю Відповідальність за порушення валютного законодавства 1. Важливим видом фінансового контролю є валютний контроль - контроль за дотриманням валютного законодавства при здійсненні валютних операцій за участю р...
7191. Вивчення модифікаційної мінливості довжини та ширини листкової пластинки берези білої (Betula alba) 231 KB
  Тема: Вивчення модифікаційної мінливості довжини та ширини листкової пластинки берези білої (Betula alba) Розділ І. Теоретична частина. Модифікаційна мінливість. Модифікаційна мінливість - це мінливість. Яка виникає внаслідок впливу факторів зо...
7192. Генетика з біометрією. Робочий зошит для лабораторних робіт 2.76 MB
  Робочий зошит для лабораторних робіт з дисципліни Генетика з біометрією для студентів ІІ курсу біолого-технологічного факультету Зміст МОДУЛЬ № 1 Цитологічні основи спадковості ЗАНЯТТЯ №1. Вступ. Методи дослідження в генетиці ЗАНЯТТЯ №2. Будова кліт...
7193. Психолого-педагогічні основи діяльності працівників пожежної охорони 704 KB
  Вступ Мета, завдання та зміст курсу Психолого-педагогічні основи діяльності працівників пожежної охорони його місце у змісті пожежної освіти. Значення психолого-педагогічних знань для спеціаліста пожежної охорони. В умовах глобальних зрушень в...
7194. Генетика и наследственность. Конспект лекций 5.47 MB
  Введение в генетику -предмет генетики, -виды наследственности, -виды изменчивости, -история развития генетики, -место генетики в системе биологических наук, -задачи ветеринарной генетики, -методы исследовании используемые в генетике...
7195. Земельное право 247.5 KB
  Общие положения. 1. Дисциплина Земельное право включена в федеральный компонент цикла общепрофессиональных дисциплин ГОС ВПО по специальности 030105.65, 030500.68 Юриспруденция. 2. В современных условиях возрастает роль земельного права в пров...