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

Для тестирования программы рассмотрим несколько решений с различными параметрами.

Решение уравнения теплопроводности в различные моменты времени. Результаты приведены в рисунках 4.2 4.4.



Рисунок 4.2 - Решение при tmax=0.2

Рисунок 4.3 - Решение при tmax=0.5

Рисунок 4.4 - Решение при tmax=1

Из рисунков видно, что с увеличением tmax уменьшается температура в центральной точке стержня.

Рассмотрим, как влияет наличие источника на решение. Результаты приведены в рисунках 4.5 4.7.

Рисунок 4.5 - Решение при отсутствии источника

Рисунок 4.6 - Решение при параметре источника А=1



Рисунок 4.7 - Решение при параметре источника А=2

Из рисунков видно, что с увеличением параметра А увеличивается температура в центральной точке стержня.

Рассмотрим, как влияют начальные условия на решение. Результаты приведены в рисунках 4.8 4.9.

Рисунок 4.8 - Решение при отсутствии начальных условий

Рисунок 4.9 - Решение при параметрах C=1 и D=0 в начальных условиях

Рассмотрим, как влияет тип уравнения на решение. Результаты приведены в рисунках 4.10 4.12.

Рисунок 4.10 - Решение при m=0



Рисунок 4.11 - Решение при m=1

Рисунок 4.12 - Решение при m=2



Заключение

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

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

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

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

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

Результаты работы представлены в сборник докладов «Творчество молодых-2014».



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

1 Несис, Е.И. Методы математической физики: учеб. пособие / Е.И. Несис. - М., 1977. - 199 с.

2 Рындин, Е.А. Решение задач математической физики в системе Matlab: учебное пособие / Е.А. Рындин, И.Е. Лысенко. - Таганрог: Изд-во ТРТУ, 2005. 62 с.

3 Говорухин, В. Ф. Компьютер в математическом исследовании / В.Ф. Говорухин, Б.П. Цибулин. // - М., 2004. - 619 с.

4 Мартынов, Н.Н. Matlab 5.x. Вычисления, визуализация, программирование / Н.Н. Мартынов, А.П. Иванов. // - М.:Кудиц-образ, 2000. 336 с.



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

% Очистка ранее сохраненных результатов вычислений

% из оперативной памяти

clear all

% Закрытие раннее выведенных графических окон

close all

% Очистка экрана

clc

% Ввод исходных данных с клавиатуры

input ('Начальная координата области решения: ' );

x0= ans;

input ('Конечная координата области решения: ' );

xn=ans;

input ('Число точек координатной сетки: ' );

n=ans;

input ('Функция правой части уравнения в одинарных кавычках: ' );

f=ans;

input ('Вид ГУ на левой границе ( 1-Дирихле, 2- Неймана): ' );

v1=ans;

input ('Значение ГУ на левой границе: ');

g1=ans;

input ('Вид ГУ на правой границе ( 1-Дирихле, 2- Неймана): ' );

v2=ans;

input ('Значение ГУ на правой границе: ');

g2=ans;

%Задание равномерной координатной сетки с шагом dx

x=x0:(xn-x0)/(n-1):xn;

dx=x(2)-x(1);

% Вычисление значений функций, заданных символьно,

% в узлах координатной сетки

F=inline(f,'x');

FF=F(x);

% Задание матрицы коэффициентов СЛАУ размерностью n x n,

% все элементы которой равны 0

a=zeros(n,n);

% Задание матрицы-строки свободных членов СЛАУ размерностью 1 x n,

% все элементы которой равны 0

b=zeros(1,n);

% Определение коэффициентов и свободных членов СЛАУ,

% соответствующих граничным условиям и проверка корректности

% значений параметров v1, v2

b(1)=g1;

if v1==1

a(1,1)=1;

elseif v1==2

a(1,1)=-1/dx;

a(1,2)=1/dx;

else

error ('Parameter v1 have incorrect value');

end

b(n)=g2;

if v2++1

if v2==1

a(n,n)=1;

elseif v2==2

a(n,n)=1/dx;

a(n,n-1)=-1/dx;

else

error ('Parameter v2 have incorrect value');

end

% Определение коэффициентов и свободных членов СЛАУ,

% соответствующих внутренним точкам области

for i=2:n-1

a(i,i)=-2/dx^2;

a(i,i+1)=1/dx^2;

a(i,i-1)=1/dx^2;

b(i)=FF(i);

end

% Решение СЛАУ

u=b/a';

% Построение графика функции правой части f(x)

plot (x,FF,x,ones(size(x)).*(sum(FF)/length(FF)),'r.','LineWidth',1.5)

xlabel('x','FontSize',13)

ylabel('f(x)','FontSize',13)

grid on

% Построение графика искомой функции u(x)

figure

plot (x,u,x,ones(size(x)).*(sum(u)/length(u)),'r.','LineWidth',1.5)

xlabel('x','FontSize',13)

ylabel('U(x)','FontSize',13)

grid on

end



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

% Функция решения волнового уравнения

% d2U/dt2=d2U/dx2+d2U/dy2

% на прямоугольной области с граничными условиями

% Дирихле и/или Неймана

function[x,y,t,U]=f_wave2d(t0,ts,s,x0,xn,n,y0,ym,m,vt1,gt1,vt2,gt2,v1,g1,v2,g2,v3,g3,v4,g4)

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

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

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

% x0 - начальная координата области решения по оси х;

% xn - конечная координата области решения по оси х;

% y0 - начальная координата области решения по оси y;

% ym - конечная координата области решения по оси y;

% n - число точек координатной сетки вдоль оси х;

% m - число точек координатной сетки вдоль оси y;

% s - число точек сетки вдоль оси времени t;

% vt1- параметр, значение которого определяет

% тип начального условия в момент времени t(1)

% (1 - Дирихле, 2 - Неймана);

% gt1 - функция в правой части начального условия

% в момент времени t(1),

% задаваемая строкой символов, заключенных

% в одинарные кавычки;

% vt2 - параметр, значение которого определяет

% тип начального условия в момент времени t(s)

% (1 - Дирихле, 2 - Неймана);

% gt2 - функция в правой части начального условия

% в момент времени t(s),

% задаваемая строкой символов, заключенных

% в одинарные кавычки;

% v1 - параметр, значение которого определяет

% тип граничного условия (ГУ) на первой границе

% области х = х(1) (1 - ГУ Дирихле, 2 - ГУ Неймана);

% g1 - функция в правой части граничного условия на первой границе,

% задаваемая строкой символов, заключенных

% в одинарные кавычки;

% v2 - параметр, значение которого определяет

% тип граничного условия на второй границе

% области х = х(n) (1 - ГУ Дирихле, 2 - ГУ Неймана);

% g2 - функция в правой части граничного условия на второй границе,

% задаваемая строкой символов, заключенных

% в одинарные кавычки;

% v3 - параметр, значение которого определяет

% тип граничного условия на третьей границе

% области y = y(1) (1 - ГУ Дирихле, 2 - ГУ Неймана);

% g3 - функция в правой части граничного условия на третьей границе,

% задаваемая строкой символов, заключенных

% в одинарные кавычки;

% v4 - параметр, значение которого определяет

% тип граничного условия на четвертой границе

% области y = y(m) (1 - ГУ Дирихле, 2 - ГУ Неймана);

% g4 - функция в правой части граничного условия на четвертой границе,

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

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

% х - вектор-строка координатной сетки по оси х размерности 1 х n;

% y - вектор-строка координатной сетки по оси y размерности 1 х m;

% t - вектор-строка сетки по оси времени размерностью 1 х s;

% U - матрица значений результирующей функции

% в узлах координатной сетки размерностью n х m x s.

% Функции и переменные по умолчанию

if exist('t0')==0

t0=0;

end

if exist('ts')==0

ts=0.2;

end

if exist('s')==0

s=6;

end

if exist('x0')==0

x0=-1;

end

if exist('xn')==0

xn=1;

end

if exist('n')==0

n=18;

end

if exist('y0')==0

y0=-1;

end

if exist('ym')==0

ym=1;

end

if exist('m')==0

m=18;

end

if exist('vt1')==0

vt1=1;

end

if exist('gt1')==0

gt1='sin(4*x)-sin(4*y)';

end

if exist('vt2')==0

vt2=1;

end

if exist('gt2')==0

gt2='sin(4*y)-sin(4*x)';

end

if exist('v1')==0

v1=2;

end

if exist('g1')==0

g1='0';

end

if exist('v2')==0

v2=2;

end

if exist('g2')==0

g2='0';

end

if exist('v3')==0

v3=2;

end

if exist('g3')==0

g3='0';

end

if exist('v4')==0

v4=2;

end

if exist('g4')==0

g4='0';

end

% Задание равномерной координатной сетки

x=x0:(xn-x0)/(n-1):xn; dx=x(2)-x(1);

y=y0:(ym-y0)/(m-1):ym; dy=y(2)-y(1);

t=t0:(ts-t0)/(s-1):ts; dt=t(2)-t(1);

% Вычисление значений функций, заданных символьно,

% в узлах координатной сетки

GT1=inline(gt1,'x','y');

GT2=inline(gt2,'x','y');

G1=inline(g1,'y');

G2=inline(g2,'y');

G3=inline(g3,'x');

G4=inline(g4,'x');

% Определение размерности СЛАУ

N=s*n*m;

% Задание матрицы коэффициентов СЛАУ размерностью N x N,

% все элементы которой равны 0

a=zeros(N,N);

% Задание матрицы-строки свободных членов СЛАУ размерностью 1 x N, matlab интерфейс теплопроводность уравнение

% все элементы которой равны 0

b=zeros(1,N);

% Определение коэффициентов и свободных членов СЛАУ,

% соответствующих начальным и граничным условиям, и проверка

% корректности значений параметров vt1, vt2, v1, v2, v3, v4

for i=1:n

for j=1:m

b(m*(i-1)+j)=GT1(x(i),y(j));

if vt1==1

a(m*(i-1)+j,m*(i-1)+j)=1;

elseif vt1==2

a(m*(i-1)+j,m*(i-1)+j)=-1/dt;

a(m*(i-1)+j,n*m+m*(i-1)+j)=1/dt;

else

error('Parameter vt1 have incorrect value');

end

b(n*m*(s-1)+m*(i-1)+j)=GT2(x(i),y(j));

if vt2==1

a(n*m*(s-1)+m*(i-1)+j,n*m*(s-1)+m*(i-1)+j)=1;

elseif vt2==2

a(n*m*(s-1)+m*(i-1)+j,n*m*(s-1)+m*(i-1)+j)=1/dt;

a(n*m*(s-1)+m*(i-1)+j,n*m*(s-2)+m*(i-1)+j)=-1/dt;

else

error('Parameter vt2 have incorrect value');

end

end

end

for l=1:s

for j=1:m

b(n*m*(l-1)+j)=G1(y(j));

if v1==1

a(n*m*(l-1)+j,n*m*(l-1)+j)=1;

elseif v1==2

a(n*m*(l-1)+j,n*m*(l-1)+j)=-1/dx;

a(n*m*(l-1)+j,n*m*(l-1)+m+j)=1/dx;

else

error('Parameter v1 have incorrect value');

end

b(n*m*(l-1)+m*(n-1)+j)=G2(y(j));

if v2==1

a(n*m*(l-1)+m*(n-1)+j,n*m*(l-1)+m*(n-1)+j)=1;

elseif v2==2

a(n*m*(l-1)+m*(n-1)+j,n*m*(l-1)+m*(n-1)+j)=1/dx;

a(n*m*(l-1)+m*(n-1)+j,n*m*(l-1)+m*(n-2)+j)=-1/dx;

else

error('Parameter v2 have incorrect value');

end

end

for i=2:n-1

b(n*m*(l-1)+m*(i-1)+1)=G3(x(i));

if v3==1

a(n*m*(l-1)+m*(i-1)+1,n*m*(l-1)+m*(i-1)+1)=1;

elseif v3==2

a(n*m*(l-1)+m*(i-1)+1,n*m*(l-1)+m*(i-1)+1)=-1/dy;

a(n*m*(l-1)+m*(i-1)+1,n*m*(l-1)+m*(i-1)+2)=1/dy;

else

error('Parameter v3 have incorrect value');

end

b(n*m*(l-1)+m*(i-1)+m)=G4(x(i));

if v4==1

a(n*m*(l-1)+m*(i-1)+m,n*m*(l-1)+m*(i-1)+m)=1;

elseif v4==2

a(n*m*(l-1)+m*(i-1)+m,n*m*(l-1)+m*(i-1)+m)=1/dy;

a(n*m*(l-1)+m*(i-1)+m,n*m*(l-1)+m*(i-1)+m-1)=-1/dy;

else

error('Parameter v4 have incorrect value');

end

end

end

% Определение коэффициентов и свободных членов СЛАУ,

% соответствующих внутренним точкам области

for l=2:s-1

for i=2:n-1

for j=2:m-1

a(n*m*(l-1)+m*(i-1)+j,n*m*(l-1)+m*(i-1)+j)=-2/dt^2+2/dx^2+2/dy^2;

a(n*m*(l-1)+m*(i-1)+j,n*m*(l-1)+m*i+j)=-1/dx^2;

a(n*m*(l-1)+m*(i-1)+j,n*m*(l-1)+m*(i-2)+j)=-1/dx^2;

a(n*m*(l-1)+m*(i-1)+j,n*m*(l-1)+m*(i-1)+j+1)=-1/dy^2;

a(n*m*(l-1)+m*(i-1)+j,n*m*(l-1)+m*(i-1)+j-1)=-1/dy^2;

a(n*m*(l-1)+m*(i-1)+j,n*m*l+m*(i-1)+j)=1/dt^2;

a(n*m*(l-1)+m*(i-1)+j,n*m*(l-2)+m*(i-1)+j)=1/dt^2;

end

end

end

% Решение СЛАУ

u=b/a';

% Преобразование вектора-строки значений искомой функции

% в узлах координатной сетки в матрицу размерности n x m x s,

% удобную для представления результатов в графическом виде

for l=1:s

for i=1:n

for j=1:m

U(i,j,l)=u(n*m*(l-1)+m*(i-1)+j);

end

end

end

% Построение графика искомой функции U(x,y,t)

for l=1:s

figure

surf(y,x,U(:,:,l))

xlabel('y','FontSize',13)

ylabel('x','FontSize',13)

zlabel('U','FontSize',13)

grid on

colormap('cool')

axis([min(y) max(y) min(x) max(x) min(min(min(U))) max(max(max(U)))])

pause(0.1)

M(l)=getframe;

end

for l=s+1:2*s-2

figure

surf(y,x,U(:,:,2*s-l))

xlabel('y','FontSize',13)

ylabel('x','FontSize',13)

zlabel('U','FontSize',13)

grid on

colormap('cool')

axis([min(y) max(y) min(x) max(x) min(min(min(U))) max(max(max(U)))])

pause(0.1)

M(l)=getframe;

end

% Отображение волнового процесса в динамическом режиме

figure

ans=1;

while ans==1

movie(M,10,10)

ans=menu('Повторить просмотр результатов?','ДА','НЕТ');

end



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

% Решение выводится в трехмерный массив SOL,

% причем первый индекс отвечает временному слою,

% второй - определяет номер узла

% третий - дает номер компоненты вектора решения

% SOL(i,j,k) даёт значение решения U^k(x^j,t^i)

function pde_demo

% m - параметр позволяющий решать одномерные задачи для радиальной

% координаты в декартовой ( m = 0), в цилиндрической или полярной( m = 1),

% а также в сферической ( m = 2) системе координат

m=0;

% a,b - концы отрезка изменения переменной x

a=0; b=3*pi/2;

% p - параметр

p=1;

% n - число ячеек на отрезке

n=16;

% xx - узлы сетки

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

% начальным и конечным элементами которого являются точки x1 и x2.

xx=linspace(a,b,n);

% tt - задает временные слои

tt=0:.2:3;

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

% ic - функция задающая начальные условия

% bc - функция задающая краевые условия

sol = pdepe (m,@pde,@ic,@bc,xx,tt,[],p);

u = sol(:,:,1);

% subplot(mnp), где mnp - 3 цифры,

%производит разбивку графического окна на несколько подокон,

% создавая при этом новые объекты axes;

% значение m указывает, на сколько частей разбивается окно по горизонтали,

% n - по вертикали, а p - номер подокна,

% куда будет выводиться очередной график.

subplot(121),

% mesh(X, Y, Z) выводит на экран сетчатую поверхность

% для значений массива Z,

% определенных на множестве значений массивов X и Y.

mesh(xx,tt,u),

% Команда axis([xmin xmax ymin ymax zmin zmax]) устанавливает масштаб

% по осям x, y, z для активного графического окна.

axis ([a b tt(1) tt(end) 0 1.5]), xlabel('x'), ylabel ('t')

subplot(122),

% plot - построение графика

plot (xx, u(end,:),'ok',xx, exp(-tt(end))*(p+sin(xx))),

% xlabel - задание нозваний координатных осей

xlabel('x'), ylabel(strcat('u(x',num2str(tt(end)),')')),

% axis tight - не позволяет изменится размеру кадра

axis tight

% Вспомогательная функция вычисления правой части

function [C,F,S]= pde(x,t,u,DuDx,p)

% Выходные данные представлены диагональной матрицей C и векторами F,C

C = 1;

F = DuDx;

S = -p*exp(-t);

% Вспомогательная функция задания начальных условий

function u0 = ic(x,p)

u0 = p+sin(x);

% Вспомогательная функция задания краевых условий

% выходные параметры представлены векторами на левом (Pa,Qa)

% и правом(Pb,Qb) концах интервала

function [Pa, Qa, Pb,Qb] = bc(a,ua,b,ub,t,p)

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

Qa = 0;

Pb = 0;

Qb = 1;



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

function MyFig1

global hFig1 hAxes

global hEd1 hEd2 hEd3 hEd4

globalhList

global hBut1 hBut2

global hTxt1 hTxt2 hTxt3 hTxt4 hTxt5 hTxt6 hTxt7

hFig1 = figure

hAxes=axes('Parent',hFig1,'Color',[1 1 1],...

'Units', 'points',...

'Position', [20 60 230 230],...

'FontSize', 10);

hEd1 = uicontrol(hFig1,'Style','edit',...

'BackgroundColor',[1 1 1 ],...

'Position', [435 400 190 30],...

'HorizontalAlignment','left');

hEd2 = uicontrol(hFig1,'Style','edit',...

'BackgroundColor',[1 1 1 ],...

'Position', [435 320 190 30],...

'HorizontalAlignment','left');

hEd3 = uicontrol(hFig1,'Style','edit',...

'BackgroundColor',[1 1 1 ],...

'Position', [435 230 190 30],...

'HorizontalAlignment','left');

hEd4 = uicontrol(hFig1,'Style','edit',...

'BackgroundColor',[1 1 1 ],...

'Position', [435 150 190 30],...

'HorizontalAlignment','left');

hList = uicontrol(hFig1,'Style','listbox',...

'Position', [485 15 140 90],...

'String', {'sin','cos','MyFunct1','MyFunct2','f3'},...

'HorizontalAlignment','left');

hBut1 = uicontrol(hFig1,'Style','pushbutton',...

'Position', [60 25 140 40],...

'String','Go',...

'Callback','MyGo');

hBut2 = uicontrol(hFig1,'Style','pushbutton',...

'Position', [240 25 140 40],...

'String','Clear',...

'Callback','MyClear');

hTxt1 = uicontrol(hFig1,'Style','Text',...

'BackgroundColor',[0.7 0.7 0.7],...

'Position', [90 475 240 20],...

'String','Euler solution');

hTxt2 = uicontrol(hFig1,'Style','Text',...

'BackgroundColor',[0.7 0.7 0.7],...

'Position', [435 440 60 20],...

'String','X initial');

hTxt3 = uicontrol(hFig1,'Style','Text',...

'BackgroundColor',[0.7 0.7 0.7],...

'Position', [435 360 60 20],...

'String','Y initial');

hTxt4 = uicontrol(hFig1,'Style','Text',...

'BackgroundColor',[0.7 0.7 0.7],...

'Position', [435 260 60 20],...

'String','X final');

hTxt5 = uicontrol(hFig1,'Style','Text',...

'BackgroundColor',[0.7 0.7 0.7],...

'Position', [435 180 60 20],...

'String','N');

hTxt6 = uicontrol(hFig1,'Style','Text',...

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

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

hTxt7 = uicontrol(hFig1,'Style','Text',...

'BackgroundColor',[0.7 0.7 0.7],...

'Position', [435 140 240 20],...

'String','Choose function');



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

function Interface_pdepe

global hFig hAxes1 hAxes2

global hEdit1 hEdit3 hEdit4 hEdit6 hEdit9 hEdit12 hEdit13 hEdit15 hEdit17

global hBut1 hBut2

global hText0 hText1 hText2 hText3 hText4 hText5 hText6 hText7 hText8 hText9 hText10 hText11 hText12 hText13 hText14 hText15 hText16 hText17

hFig = figure('Color',[0.8 0.8 0.8],...

'Name','Решение уравнения теплопроводности',...

'Position',[40 36 900 653], ...

'NumberTitle','off',...

'MenuBar','none',...

'PaperPosition',[18 180 576 432],...

'PaperUnits','points',...

'ToolBar','none');

% Создание окна для графика

hAxes1 = axes('Parent',hFig,...

'Color',[1 1 1], ...

'Units','pixels', ...

'Position',[52 296 353 321],...

'CameraUpVector',[0 1 0], ...

'CameraUpVectorMode','manual');

hAxes2 = axes('Parent',hFig,...

'Color',[1 1 1], ...

'Units','pixels', ...

'Position',[500 296 353 321],...

'CameraUpVector',[0 1 0], ...

'CameraUpVectorMode','manual');

% Создание кнопки решение

hBut1 = uicontrol('Parent',hFig,'Style','pushbutton',...

'String','Решение',...

'Position', [39 160 100 30],...

'Units','points', ...

'BackgroundColor',[0.925490196078431 0.913725490196078 0.847058823529412], ...

'ListboxTop',0, ...

'Callback','function_pdepe');

% Создание кнопки очистить графики

hBut2 = uicontrol('Parent',hFig,'Style','pushbutton',...

'String','Очистить графики',...

'Position', [540 160 100 30],...

'Units','points', ...

'BackgroundColor',[0.925490196078431 0.913725490196078 0.847058823529412], ...

'ListboxTop',0, ...

'Callback','MyClear');

% Для текста 1 <= x <=3*pi/2

hText0 = uicontrol('Parent',hFig,'Style','text', ...

'Position',[39 115 165 25], ...

'String','1 <= x <= 3*pi/2', ...

'FontSize', 15,...

'Units','points', ...

'BackgroundColor',[0.8 0.8 0.8], ...

'HorizontalAlignment','left', ...

'ListboxTop',0);

% Для текста Количество разбиений N =

hText1 = uicontrol('Parent',hFig,'Style','text', ...

'Position',[39 85 185 35], ...

'String','Количество разбиений N =', ...

'FontSize', 14,...

'Units','points', ...

'BackgroundColor',[0.8 0.8 0.8], ...

'HorizontalAlignment','left', ...

'ListboxTop',0);

% Для ввода значения N

hEdit1 = uicontrol('Parent',hFig,'Style','edit', ...

'Position',[227 103 35 15], ...

'Units','points', ...

'BackgroundColor',[1 1 1], ...

'HorizontalAlignment','left', ...

'ListboxTop',0);

% Для текста 0 <= t <= tmax

hText2 = uicontrol('Parent',hFig,'Style','text', ...

'Position',[39 80 165 19], ...

'String','0 <= t <= tmax', ...

'FontSize', 15,...

'Units','points', ...

'BackgroundColor',[0.8 0.8 0.8], ...

'HorizontalAlignment','left', ...

'ListboxTop',0);

% Для текста tmax =

hText3 = uicontrol('Parent',hFig,'Style','text', ...

'Position',[39 55 105 19], ...

'String','tmax = ', ...

'FontSize', 15,...

'Units','points', ...

'BackgroundColor',[0.8 0.8 0.8], ...

'HorizontalAlignment','left', ...

'ListboxTop',0);

% Для ввода значения tmax

hEdit3 = uicontrol('Parent',hFig,'Style','edit', ...

'Position',[89 57 35 15], ...

'Units','points', ...

'BackgroundColor',[1 1 1], ...

'HorizontalAlignment','left', ...

'ListboxTop',0);

% Для текста dt

hText4 = uicontrol('Parent',hFig,'Style','text', ...

'Position',[149 55 165 19], ...

'String','dt = ', ...

'FontSize', 15,...

'Units','points', ...

'BackgroundColor',[0.8 0.8 0.8], ...

'HorizontalAlignment','left', ...

'ListboxTop',0);

% Для ввода значения dt

hEdit4 = uicontrol('Parent',hFig,'Style','edit', ...

'Position',[184 57 35 15], ...

'Units','points', ...

'BackgroundColor',[1 1 1], ...

'HorizontalAlignment','left', ...

'ListboxTop',0);

% Для текста Тип уравнения

hText5 = uicontrol('Parent',hFig,'Style','text', ...

'Position',[39 25 105 25], ...

'String','Тип уравнения', ...

'FontSize', 15,...

'Units','points', ...

'BackgroundColor',[0.8 0.8 0.8], ...

'HorizontalAlignment','left', ...

'ListboxTop',0, ...

'Tag','StaticText2');

% Для текста m =

hText6 = uicontrol('Parent',hFig,'Style','text', ...

'Position',[149 35 105 15], ...

'String','m =', ...

'FontSize', 15,...

'Units','points', ...

'BackgroundColor',[0.8 0.8 0.8], ...

'HorizontalAlignment','left', ...

'ListboxTop',0, ...

'Tag','StaticText2');

% Для ввода значения m

hEdit6 = uicontrol('Parent',hFig,'Style','edit', ...

'Position',[184 33 35 15], ...

'Units','points', ...

'BackgroundColor',[1 1 1], ...

'HorizontalAlignment','left', ...

'ListboxTop',0, ...

'Tag','Edit_m',...

'Callback','pde_demo(''press_m'');');

% Для текста Функция источника

hText7 = uicontrol('Parent',hFig,'Style','text', ...

'Position',[300 120 145 25], ...

'String','Функция источника', ...

'FontSize', 15,...

'Units','points', ...

'BackgroundColor',[0.8 0.8 0.8], ...

'HorizontalAlignment','left', ...

'ListboxTop',0);

% Для текста F(x,t)= A*exp(-t)

hText8 = uicontrol('Parent',hFig,'Style','text', ...

'Position',[300 100 165 25], ...

'String','F(x,t)= A*exp(-t)', ...

'FontSize', 15,...

'Units','points', ...

'BackgroundColor',[0.8 0.8 0.8], ...

'HorizontalAlignment','left', ...

'ListboxTop',0);

% Для текста A =

hText9 = uicontrol('Parent',hFig,'Style','text', ...

'Position',[300 85 165 15], ...

'String','A = ', ...

'FontSize', 15,...

'Units','points', ...

'BackgroundColor',[0.8 0.8 0.8], ...

'HorizontalAlignment','left', ...

'ListboxTop',0);

% Для ввода значения A

hEdit9 = uicontrol('Parent',hFig,'Style','edit', ...

'Position',[335 83 35 15], ...

'Units','points', ...

'BackgroundColor',[1 1 1], ...

'HorizontalAlignment','left', ...

'ListboxTop',0);

% Для текста Начальные условия

hText10 = uicontrol('Parent',hFig,'Style','text', ...

'Position',[300 45 145 35], ...

'String','Начальные условия', ...

'FontSize', 15,...

'Units','points', ...

'BackgroundColor',[0.8 0.8 0.8], ...

'HorizontalAlignment','left', ...

'ListboxTop',0);

% Для текста U(x,0) = C+D*sin(x)

hText11 = uicontrol('Parent',hFig,'Style','text', ...

'Position',[300 35 165 25], ...

'String','U(x,0) = C+D*sin(x)', ...

'FontSize', 15,...

'Units','points', ...

'BackgroundColor',[0.8 0.8 0.8], ...

'HorizontalAlignment','left', ...

'ListboxTop',0);

% Для текста C =

hText12 = uicontrol('Parent',hFig,'Style','text', ...

'Position',[300 25 165 15], ...

'String','C = ', ...

'FontSize', 15,...

'Units','points', ...

'BackgroundColor',[0.8 0.8 0.8], ...

'HorizontalAlignment','left', ...

'ListboxTop',0);

% Для ввода значения C

hEdit12 = uicontrol('Parent',hFig,'Style','edit', ...

'Position',[330 23 35 15], ...

'Units','points', ...

'BackgroundColor',[1 1 1], ...

'HorizontalAlignment','left', ...

'ListboxTop',0);

% Для текста D =

hText13 = uicontrol('Parent',hFig,'Style','text', ...

'Position',[380 25 165 15], ...

'String','D = ', ...

'FontSize', 15,...

'Units','points', ...

'BackgroundColor',[0.8 0.8 0.8], ...

'HorizontalAlignment','left', ...

'ListboxTop',0);

% Для ввода значения D

hEdit13 = uicontrol('Parent',hFig,'Style','edit', ...

'Position',[410 23 35 15], ...

'Units','points', ...

'BackgroundColor',[1 1 1], ...

'HorizontalAlignment','left', ...

'ListboxTop',0);

% Для текста Температура на левой границе Pa=Ua+U1

hText14 = uicontrol('Parent',hFig,'Style','text', ...

'Position',[480 103 165 45], ...

'String','Температура на левой границе Pa=Ua+U1', ...

'FontSize', 15,...

'Units','points', ...

'BackgroundColor',[0.8 0.8 0.8], ...

'HorizontalAlignment','left', ...

'ListboxTop',0);

% Для текста U1 =

hText15 = uicontrol('Parent',hFig,'Style','text', ...

'Position',[480 85 165 15], ...

'String','U1 = ', ...

'FontSize', 15,...

'Units','points', ...

'BackgroundColor',[0.8 0.8 0.8], ...

'HorizontalAlignment','left', ...

'ListboxTop',0);

% Для ввода значения U1

hEdit15 = uicontrol('Parent',hFig,'Style','edit', ...

'Position',[515 83 35 15], ...

'Units','points', ...

'BackgroundColor',[1 1 1], ...

'HorizontalAlignment','left', ...

'ListboxTop',0)

% Для текста Поток на правой границе Qb=Ub+Q2

hText16 = uicontrol('Parent',hFig,'Style','text', ...

'Position',[480 35 165 45], ...

'String','Поток на правой границе Qb=Ub+Q2', ...

'FontSize', 15,...

'Units','points', ...

'BackgroundColor',[0.8 0.8 0.8], ...

'HorizontalAlignment','left', ...

'ListboxTop',0);

% Для текста Q2 =

hText17 = uicontrol('Parent',hFig,'Style','text', ...

'Position',[480 25 165 15], ...

'String','Q2 = ', ...

'FontSize', 15,...

'Units','points', ...

'BackgroundColor',[0.8 0.8 0.8], ...

'HorizontalAlignment','left', ...

'ListboxTop',0);

% Для ввода значения Q2

hEdit17 = uicontrol('Parent',hFig,'Style','edit', ...

'Position',[515 23 35 15], ...

'Units','points', ...

'BackgroundColor',[1 1 1], ...

'HorizontalAlignment','left', ...

'ListboxTop',0);



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

% Решение выводится в трехмерный массив SOL,

% причем первый индекс отвечает временному слою,

% второй - определяет номер узла

% третий - дает номер компоненты вектора решения

% SOL(i,j,k) даёт значение решения U^k(x^j,t^i)

function function_pdepe

global hAxes1 hAxes2

global hEdit1 hEdit3 hEdit4 hEdit6 hEdit9 hEdit12 hEdit13 hEdit15 hEdit17

global A1 C1 D1 U1 Q2

str1=get(hEdit1,'String');

str3=get(hEdit3,'String');

str4=get(hEdit4,'String');

str6=get(hEdit6,'String');

str9=get(hEdit9,'String');

str12=get(hEdit12,'String');

str13=get(hEdit13,'String');

str15=get(hEdit15,'String');

str17=get(hEdit17,'String');

A1 = str2num(str9)

C1 = str2num(str12)

D1 = str2num(str13)

% m - параметр позволяющий решать одномерные задачи для радиальной

% координаты в декартовой ( m = 0), в цилиндрической или полярной( m = 1),

% а также в сферической ( m = 2) системе координат

m = str2num(str6);

% a,b - концы отрезка изменения переменной x

a = 1;

b = 3*pi/2;

n = str2num(str1)

% p - параметры

%p = str2num(str20)

%A1 = str2num(str8)

%C1 = str2num(str11)

%D1 = str2num(str12)

U1 = str2num(str15)

Q2 = str2num(str17)

% n - число ячеек на отрезке

% xx - узлы сетки

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

% начальным и конечным элементами которого являются точки x1 и x2.

xx=linspace(a,b,n);

% tt - задает временные слои

tmax = str2num(str3)

dt = str2num(str4)

tt =0:dt:tmax

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

% ic - функция задающая начальные условия

% bc - функция задающая краевые условия

sol = pdepe(m,@pde,@ic,@bc,xx,tt);

u = sol(:,:,1);

% mesh(X, Y, Z) выводит на экран сетчатую поверхность для значений

% массива Z, определенных на множестве значений массивов X и Y.

axes(hAxes1);

mesh(xx,tt,u),

% Команда axis([xmin xmax ymin ymax zmin zmax]) устанавливает масштаб

% по осям x, y, z для активного графического окна.

axis ([a b tt(1) tt(end) 0 1.5]), xlabel('x'), ylabel ('t')

% plot - построение графика

axes(hAxes2);

plot (xx, u(end,:)),

% xlabel - задание нозваний координатных осей

xlabel('x'), ylabel(strcat('u(x',num2str(tt(end)),')')),

% axis tight - не позволяет изменится размеру кадра

axis tight

% Вспомогательная функция вычисления правой части

function [C,F,S]= pde(x,t,u,DuDx)

global A1

% Выходные данные представлены диагональной матрицей C и векторами F,C

C = 1;

F = DuDx;

S = -A1*exp(-t);

% Вспомогательная функция задания начальных условий

function u0 = ic(x)

global C1 D1

u0 =C1+D1*sin(x);

% Вспомогательная функция задания краевых условий

% выходные параметры представлены векторами на левом (Pa,Qa)

% и правом(Pb,Qb) концах интервала

function [Pa, Qa, Pb,Qb] = bc(a,ua,b,ub,t)

global U1 Q2

Pa = ua+U1;

Qa = 0;

Pb = 0;

Qb = ub+Q2;

Размещено на Allbest.ru

...