Решите один УЧП

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

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

$π^{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$ .
Выходные параметры cF, и s соответствуйте коэффициентам в стандартной форме уравнения PDE, ожидаемой pdepe. Эти коэффициенты закодированы в терминах входных переменных xTU, и 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})$ , и этот термин уже задан в основной функции УЧП, вы не должны задавать эту часть уравнения в функции граничного условия. Вы должны только задать значения $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

Выберите Solution 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, уравнение PDE, начальные условия, граничные условия и сетки для x и t.

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

pdepe возвращает решение в трехмерном массиве sol, где sol(i,j,k) аппроксимирует kкомпонент th решения $u_{k}$ оцененный в t(i) и x(j). Размер sol length(t)- length(x)- length(u0), начиная с u0 задает начальное условие для каждого компонента решения. Для этой проблемы, u имеет только один компонент, таким образом, sol 5 20 матрица, но в целом можно извлечь kкомпонент решения th с командой 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 object. The axes object 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 object. The axes object 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 object. The axes object 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

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