PDE 풀기 및 편도함수 계산하기
이 예제에서는 트랜지스터 편미분 방정식(PDE)을 풀고 그 결과를 이용하여 더 큰 문제를 푸는 과정에 필요한 편도함수를 얻는 방법을 보여줍니다.
다음과 같은 PDE가 있다고 가정해 보겠습니다.
이 방정식은 트랜지스터 이론 [1]에서 기인하며 는 PNP 트랜지스터의 베이스에서 과잉 전하 반송자(또는 정공)의 농도를 설명하는 함수입니다. 와 는 물리 상수입니다. 방정식은 시간 에 대해 구간 에서 성립합니다.
초기 조건은 상수 를 포함하며 다음과 같이 주어집니다.
이 문제는 다음과 같이 주어지는 경계 조건을 가집니다.
가 고정된 경우 방정식 의 해는 과잉 전하 붕괴를 로 설명합니다. 이러한 붕괴로 인해 다른 상수 를 가진 이미터 방전 전류가 발생합니다.
및 에 대해 에서의 경계값이 일관되지 않아 방정식은 에서 유효합니다. PDE는 에 대해 닫힌 형식의 급수해를 갖기 때문에 수치적으로뿐 아니라 해석적으로도 이미터 방전 전류를 계산하고 그 결과를 비교할 수 있습니다.
MATLAB®에서 이 문제를 풀려면 PDE 방정식, 초기 조건 및 경계 조건을 코딩하고 적합한 해 메시를 선택한 후에 솔버 pdepe
를 호출해야 합니다. 필요한 함수를 이 예제와 같이 파일 끝에 로컬 함수로 포함시킬 수도 있고, MATLAB 경로에 있는 디렉터리에 이름이 지정된 별도의 파일로 저장할 수도 있습니다.
물리 상수 정의하기
물리 상수를 추적하려면 각 상수에 대한 필드가 포함된 구조체형 배열을 만드십시오. 나중에 방정식, 초기 조건 및 경계 조건에 대한 함수를 정의할 때 함수가 상수에 액세스할 수 있도록 이 구조체를 추가 인수로 전달할 수 있습니다.
C.L = 1; C.D = 0.1; C.eta = 10; C.K = 1; C.Ip = 1;
방정식 코딩하기
방정식을 코딩하려면 먼저 방정식이 pdepe
솔버에서 요구하는 형식인지 확인해야 합니다.
이 형식으로 작성된 PDE는 다음과 같습니다.
방정식의 계수 값은 다음과 같습니다.
(각 대칭성이 없는 카테시안 좌표)
이제 방정식을 코딩하는 함수를 만들 수 있습니다. 이 함수에는 시그니처 [c,f,s] = transistorPDE(x,t,u,dudx,C)
가 있어야 합니다.
x
는 독립적 공간 변수입니다.t
는 독립적 시간 변수입니다.u
는x
와t
에 따라 미분되는 종속 변수입니다.dudx
는 편미분 공간 도함수 입니다.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
(참고: 모든 함수는 예제 끝에 로컬 함수로 포함되어 있습니다.)
초기 조건 코딩하기
이제 초기 조건을 반환하는 함수를 작성합니다. 초기 조건은 첫 번째 시간 값에서 적용되고 모든 값에 대한 값을 제공합니다. 함수 시그니처 u0 = transistorIC(x,C)
를 사용하여 함수를 작성하십시오.
초기 조건은 다음과 같습니다.
대응하는 함수는 다음과 같습니다.
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
경계 조건 코딩하기
이제 경계 조건 을 계산하는 함수를 작성합니다. 구간 에 있는 문제의 경우, 경계 조건은 모든 , 그리고 또는 에 적용됩니다. 솔버에서 요구하는 경계 조건의 표준 형식은 다음과 같습니다.
이 문제의 경계 조건을 표준 형식으로 작성하면 다음과 같습니다.
- 인 경우, 방정식은 입니다. 계수는 다음과 같습니다.
- 마찬가지로 인 경우, 방정식은 입니다. 계수는 다음과 같습니다.
경계 함수는 함수 시그니처 [pl,ql,pr,qr] = transistorBC(xl,ul,xr,ur,t)
를 사용해야 합니다.
입력값 x
l
과 ul
은 왼쪽 경계의 와 에 대응합니다.입력값
xr
과ur
은 오른쪽 경계의 와 에 대응합니다.t
는 독립적 시간 변수입니다.출력값
pl
과ql
은 왼쪽 경계(이 문제의 경우 )의 와 에 대응합니다.출력값
pr
과qr
은 오른쪽 경계(이 문제의 경우 )의 와 에 대응합니다.
이 예제의 경계 조건은 다음과 같은 함수로 표현됩니다.
function [pl,ql,pr,qr] = transistorBC(xl,ul,xr,ur,t) pl = ul; ql = 0; pr = ur; qr = 0; end
해 메시 선택하기
해 메시는 솔버가 해를 계산하는 지점인 와 의 값을 정의합니다. 이 문제에 대한 해는 급격하게 변하므로 구간 의 공간 점 50개와 구간 의 시간 점 50개로 구성된, 상대적으로 세밀한 메시를 사용하십시오.
x = linspace(0,C.L,50); t = linspace(0,1,50);
방정식 풀기
마지막으로 대칭 , PDE 방정식, 초기 조건, 경계 조건, 와 에 대한 메시를 사용하여 방정식을 풉니다. 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
는 하나의 성분만 가지지만, 일반적으로 u = sol(:,:,k)
를 사용하여 k
번째 해 성분을 추출할 수 있습니다.
u = sol(:,:,1);
해 플로팅하기
와 에 대해 선택한 메시 점에서 플로팅된 해 의 곡면 플롯을 만듭니다.
surf(x,t,u) title('Numerical Solution (50 mesh points)') xlabel('Distance x') ylabel('Time t') zlabel('Solution u(x,t)')
이제 곡면 플롯의 등고선을 측면에서 볼 수 있도록 와 만 플로팅합니다.
plot(x,u) xlabel('Distance x') ylabel('Solution u(x,t)') title('Solution profiles at several times')
이미터 방전 전류 계산하기
에 대한 급수해를 사용하면 이미터 방전 전류를 다음과 같이 무한 급수 [1]로 표현할 수 있습니다.
급수에 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
pdepe
로 계산한 것처럼 에 대한 수치 해를 사용하면, 다음 방정식으로 에서 에 대한 수치 근삿값을 계산할 수도 있습니다.
에 대한 해석적 해와 수치 해를 계산하고 결과를 플로팅합니다. pdeval
을 사용하여 에서 값을 계산하십시오.
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')
결과는 서로 상당히 일치합니다. 더 세밀한 해 메시를 사용하여 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 % ----------------------------------------------------
참고 문헌
[1] Zachmanoglou, E.C.와 D.L. Thoe. Introduction to Partial Differential Equations with Applications. Dover, New York, 1986.