PDE 풀기 및 편도함수 계산하기

이 예제에서는 트랜지스터 편미분 방정식(PDE)을 풀고 그 결과를 이용하여 더 큰 문제를 푸는 과정에 필요한 편도함수를 얻는 방법을 보여줍니다.

다음과 같은 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$ 와 $η$ 는 물리 상수입니다. 방정식은 시간 $t \geq 0$ 에 대해 구간 $0 \leq x \leq L$ 에서 성립합니다.

초기 조건은 상수 $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$ 및 $t > 0$ 에 대해 $x = 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 (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에서 요구하는 표준 PDE 방정식 형식의 계수에 대응합니다.

그 결과, 이 예제의 방정식은 다음과 같은 함수로 표현할 수 있습니다.

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

(참고: 모든 함수는 예제 끝에 로컬 함수로 포함되어 있습니다.)

초기 조건 코딩하기

이제 초기 조건을 반환하는 함수를 작성합니다. 초기 조건은 첫 번째 시간 값에서 적용되고 모든 $x$ 값에 대한 $u (x, t_{0})$ 값을 제공합니다. 함수 시그니처 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은 왼쪽 경계(이 문제의 경우 $x = 0$ )의 $p_{L} (x, t, u)$ 와 $q_{L} (x, t)$ 에 대응합니다.
출력값 pr과 qr은 오른쪽 경계(이 문제의 경우 $x = 1$ )의 $p_{R} (x, t, u)$ 와 $q_{R} (x, t)$ 에 대응합니다.

이 예제의 경계 조건은 다음과 같은 함수로 표현됩니다.

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

해 메시 선택하기

해 메시는 솔버가 해를 계산하는 지점인 $x$ 와 $t$ 의 값을 정의합니다. 이 문제에 대한 해는 급격하게 변하므로 구간 $0 \leq x \leq L$ 의 공간 점 50개와 구간 $0 \leq t \leq 1$ 의 시간 점 50개로 구성된, 상대적으로 세밀한 메시를 사용하십시오.

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

방정식 풀기

마지막으로 대칭 $m$ , PDE 방정식, 초기 조건, 경계 조건, $x$ 와 $t$ 에 대한 메시를 사용하여 방정식을 풉니다. pdepe에서 PDE 함수는 4개의 입력값을 사용하고 초기 조건 함수는 1개의 입력값을 사용해야 하므로, 물리 상수의 구조체를 추가 입력값으로 전달하는 함수 핸들을 만드십시오.

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차원 배열 sol로 해를 반환합니다. 여기서 sol(i,j,k)는 t(i)와 x(j)에서 구한 해의 k번째 성분 $u_{k}$ 에 대한 근삿값을 계산합니다. 이 문제의 경우 u는 하나의 성분만 가지지만, 일반적으로 u = sol(:,:,k)를 사용하여 k번째 해 성분을 추출할 수 있습니다.

u = sol(:,:,1);

해 플로팅하기

$x$ 와 $t$ 에 대해 선택한 메시 점에서 플로팅된 해 $u$ 의 곡면 플롯을 만듭니다.

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), xlabel Distance x, ylabel Time t 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, xlabel Distance x, ylabel Solution u(x,t) 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)} .$

급수에 40개 항을 사용하는 $I (t)$ 에 대한 해석적 해를 계산하는 함수를 작성합니다. 유일한 변수는 시간이지만, 상수의 구조체에 함수에 대한 두 번째 입력값을 지정하십시오.

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

pdepe로 계산한 것처럼 $u (x, t)$ 에 대한 수치 해를 사용하면, 다음 방정식으로 $x = 0$ 에서 $I (t)$ 에 대한 수치 근삿값을 계산할 수도 있습니다.

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

$I (t)$ 에 대한 해석적 해와 수치 해를 계산하고 결과를 플로팅합니다. pdeval을 사용하여 $x = 0$ 에서 $\partial u / \partial x$ 값을 계산하십시오.

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), xlabel Time t contains 2 objects of type line. One or more of the lines displays its values using only markers These objects represent From PDEPE + PDEVAL, From series.

결과는 서로 상당히 일치합니다. 더 세밀한 해 메시를 사용하여 pdepe의 수치 결과를 보다 더 개선할 수 있습니다.

로컬 함수(Local Function)

여기 나열된 함수는 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
% ----------------------------------------------------