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

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

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

$π^{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 соответствуйте коэффициентам в стандартной форме уравнения УЧП, ожидаемой 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, уравнение УЧП, начальные условия, граничные условия и сетки для 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')

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

$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')

Теперь сравните числовые и аналитические решения в $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)')

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

Перечисленный здесь локальные функции помощника, которые решатель УЧП 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

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