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

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

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

π2ut=2ux2.

Уравнение определено на интервале 0x1 в течение многих времен t0. В t=0, решение удовлетворяет начальному условию

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

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

u(0,t)=0,

πe-t+ux(1,t)=0.

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

Уравнение кода

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

c(x,t,u,ux)ut=x-mx(xmf(x,t,u,ux))+s(x,t,u,ux).

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

π2ut=x0x(x0ux)+0.

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

m=0

c(x,t,u,ux)=π2

f(x,t,u,ux)=ux

s(x,t,u,ux)=0

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

  • x независимая пространственная переменная.

  • t независимая переменная времени.

  • u зависимая переменная, дифференцируемая относительно x и t.

  • dudx частичная пространственная производная u/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

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

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

p(x,t,u)+q(x,t)f(x,t,u,ux)=0.

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

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

u(0,t)=0u+0ux=0.

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

  • p(0,t,u)=u

  • q(0,t)=0

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

πe-t+ux(1,t)=0πe-t+1ux(1,t)=0.

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

  • p(1,t,u)=πe-t

  • q(1,t)=1

Поскольку функция граничного условия описывается в терминах f(x,t,u,ux), и этот термин уже задан в основной функции УЧП, вы не должны задавать эту часть уравнения в функции граничного условия. Вы должны только задать значения 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 решения 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')

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

u(x,t)=e-tsin(πx).

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

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

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

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

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

Смотрите также

Похожие темы

Для просмотра документации необходимо авторизоваться на сайте