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

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

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

Рассмотрим PDE

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

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

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

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

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

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

Для фиксированного $x$ решение уравнения $u (x,$ t) описывает свертывание избыточного заряда $как$ t→∞. Это коллапс производит ток, называемый разрядным током эмиттера, который имеет другую ${постоянную}_{}$ Ip:

$I (t) {\frac{_{}}{} \frac{}{} =[IpDK∂∂xu (x,}_{t)}]$ x = 0.

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

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

Определение физических констант

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

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

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

Прежде чем кодировать уравнение, необходимо убедиться, что оно в форме, 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{}{}^{∂u∂t=x0} \frac{}{∂∂x} (^{} x0 \frac{}{} D\partialu\partialx) \frac{-}{} \frac{}{}$ DηL∂u∂x.

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

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

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

x является независимой пространственной переменной.
t - независимая переменная времени.
u является зависимой переменной, дифференцируемой по отношению к x и t.
dudx - частная пространственная производная $∂u/∂x$ .
C является дополнительным вводом, содержащим физические константы.
Продукция c, f, и 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,_{}$ t0) для любого значения x. Используйте сигнатуру функцииu0 = transistorIC(x,C) для записи функции.

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

$u (x, 0) \frac{=K}{} \frac{LD (^{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≤x≤b, граничные условия применяются для всех t и $либо x$ = a, либо x = b. Стандартная форма для граничных условий, ожидаемых решателем:

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

Записанные в этой форме граничные условия для этой задачи

- Для $x =$ 0 уравнение равно $\frac{}{} u+0\cdotd\partialu\partialx=0$ . Коэффициенты:

$_{pL} (x, t, u)$ = u,
$_{qL} (x, t)$ = 0.

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

$_{pR} (x, t, u)$ = u,
$_{qR} (x, t)$ = 0.

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

Входные данные xl и ul соответствуют $x$ и $u$ для левой границы.
Исходные данные xr и ur соответствуют $x$ и $u$ для правой границы.
t - независимая переменная времени.
Продукция pl и ql соответствуют $_{pL} (x, t,$ u) $_{и} qL ($ x, t) для левой $границы$ (x = 0 для этой задачи).
Продукция pr и qr соответствуют $_{pR} (x, t,$ u) $_{и} qR ($ x, t) для правой $границы$ (x = 1 для этой задачи).

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

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

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

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

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

Уравнение решения

Наконец, решите уравнение с помощью симметрии $m$ , уравнения PDE, начального условия, граничных условий и сеток для $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 возвращает решение в массиве 3-D sol, где sol(i,j,k) аппроксимирует k-й компонент раствора, $_{}$ оцененный при t(i) и x(j). Для этой проблемы u имеет только один компонент, но в целом можно извлечь kКомпонент решения th с командой u = sol(:,:,k).

u = sol(:,:,1);

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

Создайте график поверхности решения $u$ , выводимого на печать в выбранных точках сетки для $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) =^{}_{} \frac{2π2Ip^{(}}{}_{1-e-, (1-e))}^{} \frac{^{}}{^{}^{}^{}}^{\frac{}{^{}}^{}^{∑n=1∞n2n2π2+η2/4e-dtL2} (^{} n2xeon2}$ +, 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{_{}}{} \frac{}{} =[IpDK∂∂xu (x,}_{t)}]$ x = 0.

Вычислите аналитические и числовые решения для $I (t$ ) и постройте график результатов. Использоватьpdeval для вычисления значения $∂u/∂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 используя более тонкую сетку раствора.

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

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

См. также

pdepe

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