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

Скрипт Open 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, необходимо закодировать уравнение PDE, начальные условия и граничные условия, затем выбрать подходящую 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}) .$

В этой форме УЧП

$\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 дополнительный вход, содержащий физические константы.
Выходные параметры cF, и s соответствуйте коэффициентам в стандартной форме уравнения PDE, ожидаемой 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):

Входные параметры xl и ul соответствовать $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

Выберите Solution Mesh

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

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

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

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

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компонент th решения $u_{k}$ оцененный в t(i) и x(j). Для этой проблемы u имеет только один компонент, но в целом можно извлечь kкомпонент решения th с командой 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 object. The axes object 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 object. The axes object 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 object. The axes object with title Emitter discharge current I(t) contains 2 objects of type line. These objects represent From PDEPE + PDEVAL, From series.

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

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

Перечисленный здесь локальные функции помощника что решатель УЧП 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, Э.К. и Д.Л. То. Введение в дифференциальные уравнения с частными производными с приложениями. Дувр, Нью-Йорк, 1986.

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

pdepe

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