Решение УЧП и вычисление частичных производных

Открыть Live Script

Этот пример показывает, как решить транзистор с дифференциальным уравнением с частными производными (PDE) и использовать результаты, чтобы получить производные с частными производными, которые являются частью решения большей задачи.

Рассмотрите УЧП

$\frac{\partial u}{\partial t} = D \frac{\partial^{2} u}{\partial x^{2}} - \frac{D η}{L} \frac{\partial u}{\partial x} .$

Это уравнение возникает в теории транзисторов [1], и $u (x, t)$ - функция, описывающая концентрацию носителей избыточного заряда (или отверстий) в основе PNP транзистора. $D$ и $η$ являются физическими константами. Уравнение удерживается на интервале $0 \leq x \leq L$ для раз $t \geq 0$ .

Начальное условие включает константу $K$ и дается

$u (x, 0) = \frac{K L}{D} (\frac{1 - e^{- η (1 - x / L)}}{η}) .$

Задача имеет граничные условия, заданные как

$u (0, t) = u (L, t) = 0 .$

Для фиксированного $x$ , решение уравнения $u (x, t)$ описывает сверновение избыточного заряда как $t \to \infty$ . Это сверновение производит ток, называемый током разряда излучателя, который имеет другую постоянную $I_{p}$ :

$I (t) = {[\frac{I_{p} D}{K} \frac{\partial}{\partial x} u (x, t)]}_{x = 0} .$

Уравнение верно для $t > 0$ из-за несоответствия граничных значений в $x = 0$ для $t = 0$ и $t > 0$ . Поскольку УЧП имеет решение серии замкнутой формы для $u (x, t)$ можно вычислить ток разряда излучателя как аналитически, так и численно, и сравнить результаты.

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

Задайте физические константы

Чтобы отслеживать физические константы, создайте массив структур с полями для каждого. Когда вы позже задаете функции для уравнений, начального условия и граничных условий, можно пройти в этой структуре как дополнительный аргумент, чтобы функции имели доступ к константам.

C.L = 1;
C.D = 0.1;
C.eta = 10;
C.K = 1;
C.Ip = 1;

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

Прежде чем вы сможете кодировать уравнение, необходимо убедиться, что это в том виде, в котором 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}) .$

В этой форме PDE является

$\frac{\partial u}{\partial t} = x^{0} \frac{\partial}{\partial x} (x^{0} D \frac{\partial u}{\partial x}) - \frac{D η}{L} \frac{\partial u}{\partial x} .$

Таким образом, значения коэффициентов в уравнении

$m = 0$ (Декартовы координаты без угловой симметрии)
$c (x, t, u, \frac{\partial u}{\partial x}) = 1$
$f (x, t, u, \frac{\partial u}{\partial x}) = D \frac{\partial u}{\partial x}$
$s (x, t, u, \frac{\partial u}{\partial x}) = - \frac{D η}{L} \frac{\partial u}{\partial x}$

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

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

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

function [c,f,s] = transistorPDE(x,t,u,dudx,C)
D = C.D;
eta = C.eta;
L = C.L;

c = 1;
f = D*dudx;
s = -(D*eta/L)*dudx;
end

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

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

Затем напишите функцию, которая возвращает начальное условие. Начальное условие применяется в первый раз и предоставляет значение $u (x, t_{0})$ для любого значения $x$ . Используйте сигнатуру функции u0 = transistorIC(x,C) для записи функции.

Начальное условие является

$u (x, 0) = \frac{K L}{D} (\frac{1 - e^{- η (1 - x / L)}}{η}) .$

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

function u0 = transistorIC(x,C)
K = C.K;
L = C.L;
D = C.D;
eta = C.eta;

u0 = (K*L/D)*(1 - exp(-eta*(1 - x/L)))/eta;
end

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

Теперь напишите функцию, которая оценивает граничные условия $u (0, t) = u (1, t) = 0$ . Для задач, поставленных на интервале $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 \cdot d \frac{\partial u}{\partial x} = 0 .$ Коэффициенты:

$p_{L} (x, t, u) = u,$
$q_{L} (x, t) = 0 .$

- Аналогично для $x = 1$ , уравнение $u + 0 \cdot d \frac{\partial u}{\partial x} = 0 .$ Коэффициенты:

$p_{R} (x, t, u) = u,$
$q_{R} (x, t) = 0 .$

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

Входы x l и u l соответствуют $x$ и $u$ для левого контура.
Входы xr и ur соответствуют $x$ и $u$ для правого контура.
t - независимая временная переменная.
Выходные выходы pl и ql соответствуют $p_{L} (x, t, u)$ и $q_{L} (x, t)$ для левого контура ( $x = 0$ для этой задачи).
Выходные выходы pr и qr соответствуют $p_{R} (x, t, u)$ и $q_{R} (x, t)$ для правого контура ( $x = 1$ для этой задачи).

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

function [pl,ql,pr,qr] = transistorBC(xl,ul,xr,ur,t)
pl = ul;
ql = 0;
pr = ur;
qr = 0;
end

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

Mesh решения определяет значения $x$ и $t$ где решение оценивается решателем. Поскольку решение этой задачи изменяется быстро, используйте относительно мелкий mesh из 50 пространственных точек в интервале $0 \leq x \leq L$ и 50 временных точек в интервале $0 \leq t \leq 1$ .

x = linspace(0,C.L,50);
t = linspace(0,1,50);

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

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

m = 0;
eqn = @(x,t,u,dudx) transistorPDE(x,t,u,dudx,C);
ic = @(x) transistorIC(x,C);
sol = pdepe(m,eqn,ic,@transistorBC,x,t);

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

u = sol(:,:,1);

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

Создайте объемную поверхностную диаграмму решения $u$ нанесенный на график в выбранных точках mesh для $x$ и $t$ .

surf(x,t,u)
title('Numerical Solution (50 mesh points)')
xlabel('Distance x')
ylabel('Time t')
zlabel('Solution u(x,t)')

Figure contains an axes. The axes with title Numerical Solution (50 mesh points) contains an object of type surface.

Теперь постройте график просто $x$ и $u$ для получения вида сбоку контуров на объемной поверхностной диаграмме.

plot(x,u)
xlabel('Distance x')
ylabel('Solution u(x,t)')
title('Solution profiles at several times')

Figure contains an axes. The axes with title Solution profiles at several times contains 50 objects of type line.

Вычисление тока разряда излучателя

Использование последовательного решения для $u (x, t)$ ток разрядки излучателя может быть выражен как бесконечный ряд [1]:

$I (t) = 2 π^{2} I_{p} (\frac{1 - e^{- η}}{η}) \sum_{n = 1}^{\infty} \frac{n^{2}}{n^{2} π^{2} + η^{2} / 4} e^{- \frac{dt}{L^{2}} (n^{2} π^{2} + η^{2} / 4)} .$

Напишите функцию, чтобы вычислить аналитическое решение для $I (t)$ использование 40 терминов в серии. Единственная переменная - это время, но задайте второй вход в функцию для структуры констант.

function It = serex3(t,C) % Approximate I(t) by series expansion.
Ip = C.Ip;
eta = C.eta;
D = C.D;
L = C.L;

It = 0;
for n = 1:40 % Use 40 terms
  m = (n*pi)^2 + 0.25*eta^2;
  It = It + ((n*pi)^2 / m)* exp(-(D/L^2)*m*t);
end
It = 2*Ip*((1 - exp(-eta))/eta)*It;
end

Использование числового решения для $u (x, t)$ как вычислено pdepe, можно также вычислить числовое приближение для $I (t)$ в $x = 0$ с

$I (t) = {[\frac{I_{p} D}{K} \frac{\partial}{\partial x} u (x, t)]}_{x = 0} .$

Вычислим аналитические и числовые решения для $I (t)$ и постройте график результатов. Использование pdeval для вычисления значения $\partial u / \partial x$ в $x = 0$ .

nt = length(t);
I = zeros(1,nt);
seriesI = zeros(1,nt);
iok = 2:nt;
for j = iok
   % At time t(j), compute du/dx at x = 0.
   [~,I(j)] = pdeval(m,x,u(j,:),0);
   seriesI(j) = serex3(t(j),C);
end
% Numeric solution has form I(t) = (I_p*D/K)*du(0,t)/dx
I = (C.Ip*C.D/C.K)*I;

plot(t(iok),I(iok),'o',t(iok),seriesI(iok))
legend('From PDEPE + PDEVAL','From series')
title('Emitter discharge current I(t)')
xlabel('Time t')

Figure contains an axes. The axes with title Emitter discharge current I(t) contains 2 objects of type line. These objects represent From PDEPE + PDEVAL, From series.

Результаты соответствуют достаточно хорошо. Можно дополнительно улучшить численный результат от pdepe при использовании более мелкого mesh решения.

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

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

function [c,f,s] = transistorPDE(x,t,u,dudx,C) % Equation to solve
D = C.D;
eta = C.eta;
L = C.L;

c = 1;
f = D*dudx;
s = -(D*eta/L)*dudx;
end
% ----------------------------------------------------
function u0 = transistorIC(x,C) % Initial condition
K = C.K;
L = C.L;
D = C.D;
eta = C.eta;

u0 = (K*L/D)*(1 - exp(-eta*(1 - x/L)))/eta;
end
% ----------------------------------------------------
function [pl,ql,pr,qr] = transistorBC(xl,ul,xr,ur,t) % Boundary conditions
pl = ul;
ql = 0;
pr = ur;
qr = 0;
end
% ----------------------------------------------------
function It = serex3(t,C) % Approximate I(t) by series expansion.
Ip = C.Ip;
eta = C.eta;
D = C.D;
L = C.L;

It = 0;
for n = 1:40 % Use 40 terms
  m = (n*pi)^2 + 0.25*eta^2;
  It = It + ((n*pi)^2 / m)* exp(-(D/L^2)*m*t);
end
It = 2*Ip*((1 - exp(-eta))/eta)*It;
end
% ----------------------------------------------------

Ссылки

[1] Zachmanoglou, E.C. and D.L. Thoe. Введение в дифференциальные уравнения с частными производными с применением. Дувр, Нью-Йорк, 1986.

См. также

pdepe

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