В. Пример 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;