Решение одиночного УЧП

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

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

$π^{2} \frac{\partial u}{\partial t} = \frac{\partial^{2} u}{\partial x^{2}} .$

Уравнение определяется интервалом $0 \leq x \leq 1$ для раз $t \geq 0$ . В $t = 0$ , решение удовлетворяет начальному условию

$u (x, 0) = \sin (π x) .$

Кроме того, в $x = 0$ и $x = 1$ , решение удовлетворяет граничным условиям

$u (0, t) = 0,$

$π e^{- t} + \frac{\partial u}{\partial x} (1, t) = 0 .$

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

Кодовое уравнение

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

$c (x, t, u, \frac{\partial u}{\partial x}) \frac{\partial u}{\partial t} = x^{- m} \frac{\partial}{\partial x} (x^{m} f (x, t, u, \frac{\partial u}{\partial x})) + s (x, t, u, \frac{\partial u}{\partial x})$ .

Написанный в этой форме, УЧП становится

$π^{2} \frac{\partial u}{\partial t} = x^{0} \frac{\partial}{\partial x} (x^{0} \frac{\partial u}{\partial x}) + 0$ .

С помощью уравнения в надлежащей форме можно считать релевантные условия:

$m = 0$

$c (x, t, u, \frac{\partial u}{\partial x}) = π^{2}$

$f (x, t, u, \frac{\partial u}{\partial x}) = \frac{\partial u}{\partial x}$

$s (x, t, u, \frac{\partial u}{\partial x}) = 0$

Теперь можно создать функцию, чтобы кодировать уравнение. Функция должна иметь подпись [c,f,s] = pdex1pde(x,t,u,dudx):

x - независимая пространственная переменная.
t - независимая временная переменная.
u - зависимая переменная, дифференцируемая по отношению к x и t.
dudx - частная пространственная производная $\partial u / \partial x$ .
Выходные выходы c, f, и s соответствуют коэффициентам в стандартной форме уравнения УЧП, ожидаемым по pdepe. Эти коэффициенты кодируются в терминах входных переменных x, t, u, и dudx.

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

function [c,f,s] = pdex1pde(x,t,u,dudx)
c = pi^2;
f = dudx;
s = 0;
end

(Примечание. Все функции включены в качестве локальных функций в конце примера.)

Начальное условие кода

Затем напишите функцию, которая возвращает начальное условие. Начальное условие применяется при первом временном значении tspan(1). Функция должна иметь подпись u0 = pdex1ic(x).

Соответствующая функция является

function u0 = pdex1ic(x)
u0 = sin(pi*x);
end

Граничные условия кода

Теперь напишите функцию, которая оценивает граничные условия. Для задач, поставленных на интервале $a \leq x \leq b$ , граничные условия применяются для всех $t$ и либо $x = a$ или $x = b$ . Стандартная форма для граничных условий, ожидаемых решателем,

$p (x, t, u) + q (x, t) f (x, t, u, \frac{\partial u}{\partial x}) = 0 .$

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

Для $x = 0$ , уравнение

$u (0, t) = 0 \to u + 0 \cdot \frac{\partial u}{\partial x} = 0 .$

Коэффициенты:

$p (0, t, u) = u$
$q (0, t) = 0$

Для $x = 1$ , уравнение

$π e^{- t} + \frac{\partial u}{\partial x} (1, t) = 0 \to π e^{- t} + 1 \cdot \frac{\partial u}{\partial x} (1, t) = 0 .$

Коэффициенты:

$p (1, t, u) = π e^{- t}$
$q (1, t) = 1$

Поскольку функция граничного условия выражена в терминах $f (x, t, u, \frac{\partial u}{\partial x})$ , и этот термин уже определен в основной функции PDE, вам не нужно задавать эту часть уравнения в функции граничных условий. Вам нужно только задать значения $p (x, t, u)$ и $q (x, t)$ на каждом контуре.

Граничная функция должна использовать сигнатуру функции [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t):

Входы xl и ul соответствуют $u$ и $x$ для левого контура.
Входы xr и ur соответствуют $u$ и $x$ для правого контура.
t - независимая временная переменная.
Выходные выходы pl и ql соответствуют $p (x, t, u)$ и $q (x, t)$ для левого контура ( $x = 0$ для этой задачи).
Выходные выходы pr и qr соответствуют $p (x, t, u)$ и $q (x, t)$ для правого контура ( $x = 1$ для этой задачи).

Граничные условия в этом примере представлены функцией:

function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)
pl = ul;
ql = 0;
pr = pi * exp(-t);
qr = 1;
end

Выберите Mesh решения

Перед решением уравнения необходимо задать точки mesh $(t, x)$ при котором вы хотите pdepe для оценки решения. Задайте точки как векторы t и x. Векторы t и x играть различные роли в решателе. В частности, стоимость и точность решения сильно зависят от длины вектора x. Однако расчет намного менее чувствительно к значениям в векторе t.

Для этой задачи используйте mesh с 20 равномерно разнесенными точками в пространственном интервале [0,1] и пятью значениями t из временного интервала [0,2].

x = linspace(0,1,20);
t = linspace(0,2,5);

Решение уравнения

Наконец, решить уравнение используя симметрию m, уравнение УЧП, начальные условия, граничные условия и сетки для x и t.

m = 0;
sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t);

pdepe возвращает решение в трехмерный массив sol, где sol(i,j,k) аппроксимирует kпервый компонент решения $u_{k}$ оценивается в t(i) и x(j). Размер sol является length(t)-by- length(x)-by- length(u0), с u0 задает начальное условие для каждого компонента решения. Для этой задачи u имеет только один компонент, так что sol является матрицей 5 на 20, но в целом можно извлечь kкомпонент решения с командой u = sol(:,:,k).

Извлеките первый компонент решения из sol.

u = sol(:,:,1);

Решение для построения графика

Создать объемную поверхностную диаграмму решения.

surf(x,t,u)
title('Numerical solution computed with 20 mesh points')
xlabel('Distance x')
ylabel('Time t')

Figure contains an axes. The axes with title Numerical solution computed with 20 mesh points contains an object of type surface.

Начальное условие и краевые условия этой задачи были выбраны так, чтобы было аналитическое решение, заданное как

$u (x, t) = e^{- t} \sin (π x) .$

Постройте график аналитического решения с теми же точками mesh.

surf(x,t,exp(-t)'*sin(pi*x))
title('True solution plotted with 20 mesh points')
xlabel('Distance x')
ylabel('Time t')

Figure contains an axes. The axes with title True solution plotted with 20 mesh points contains an object of type surface.

Теперь сравните числовые и аналитические решения в $t_{f}$ , окончательное значение $t$ . В этом примере $t_{f} = 2$ .

plot(x,u(end,:),'o',x,exp(-t(end))*sin(pi*x))
title('Solution at t = 2')
legend('Numerical, 20 mesh points','Analytical','Location','South')
xlabel('Distance x')
ylabel('u(x,2)')

Figure contains an axes. The axes with title Solution at t = 2 contains 2 objects of type line. These objects represent Numerical, 20 mesh points, Analytical.

Локальные функции

Здесь перечислены локальные вспомогательные функции, которые УЧП решателя pdepe вызывает для вычисления решения. Также можно сохранить эти функции как собственные файлы в директории по пути MATLAB.

function [c,f,s] = pdex1pde(x,t,u,dudx) % Equation to solve
c = pi^2;
f = dudx;
s = 0;
end
%----------------------------------------------
function u0 = pdex1ic(x) % Initial conditions
u0 = sin(pi*x);
end
%----------------------------------------------
function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t) % Boundary conditions
pl = ul;
ql = 0;
pr = pi * exp(-t);
qr = 1;
end
%----------------------------------------------

См. также

pdepe

Документация