Решение одного PDE

Открыть сценарий в реальном времени

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

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

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

Уравнение определяется на интервале, $0≤x≤1$ для времени $t≥0$ . При $t =$ 0 решение удовлетворяет начальному условию

$u (x, 0) = sin$ (øx).

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

$u (0, t)$ = 0,

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

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

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

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

$c (\frac{}{} x,t,u,\partialu\partialx \frac{)}{}^{} \frac{}{} {∂u∂t=x-m∂∂x}^{} (xmf (\frac{}{} x,t,u,∂u∂x)) + s \frac{(}{}$ x,t,u,∂u∂x).

Написанный в этой форме PDE становится

$^{} \frac{}{}^{} \frac{}{π2∂u∂t=x0∂∂x} (^{} \frac{}{} x0\partialu\partialx)$ + 0.

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

$m =$ 0

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

$f (\frac{}{} x,t,u,\partialu\partialx) \frac{}{}$ =∂u∂x

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

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

x является независимой пространственной переменной.
t - независимая переменная времени.
u является зависимой переменной, дифференцируемой по отношению к x и t.
dudx - частная пространственная производная $∂u/∂x$ .
Продукция c, f, и s соответствуют коэффициентам в стандартной форме уравнения PDE, ожидаемым 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≤x≤b$ , граничные условия применяются для всех $t$ и $либо x =$ a, $либо x$ = b. Стандартная форма для граничных условий, ожидаемых решателем:

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

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

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

$u (0, t) = 0 → \frac{}{}$ u+0⋅∂u∂x=0.

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

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

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

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

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

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

Поскольку функция граничного условия выражается в терминах $f (\frac{}{} x,t,u,\partialu\partialx$ ), и этот термин уже определен в основной функции 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

Выбрать сетку решения

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

Для этой задачи используйте сетку с 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 возвращает решение в массиве 3-D sol, где sol(i,j,k) аппроксимирует kтретий компонент раствора $_{uk}$ оценен на 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. The axes with title Numerical solution computed with 20 mesh points contains an object of type surface.

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

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

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

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.

Теперь сравните численные и аналитические решения при $_{tf}$ , конечное значение t. В этом примере $_{tf} =$ 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 решателя PDE для вычисления решения. Кроме того, эти функции можно сохранить в виде собственных файлов в каталоге по пути 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

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