연립 PDE 풀기
이 예제에서는 2개의 편미분 방정식으로 구성된 연립방정식의 해를 정식화, 계산 및 플로팅하는 방법을 보여줍니다.
연립 PDE가 있다고 가정해 보겠습니다.
(함수 는 약식으로 사용됩니다.)
방정식은 시간 에 대해 구간 에서 성립합니다. 초기 조건은 다음과 같습니다.
경계 조건은 다음과 같습니다.
MATLAB®에서 이 방정식을 풀려면 방정식, 초기 조건 및 경계 조건을 코딩하고 적합한 해 메시를 선택한 후에 솔버 pdepe
를 호출해야 합니다. 필요한 함수를 이 예제와 같이 파일 끝에 로컬 함수로 포함시킬 수도 있고, MATLAB 경로에 있는 디렉터리에 이름이 지정된 별도의 파일로 저장할 수도 있습니다.
방정식 코딩하기
방정식을 코딩하려면 먼저 방정식이 pdepe
솔버에서 요구하는 형식인지 확인해야 합니다.
이 형식에서 PDE 계수는 행렬 값을 가지며 방정식은 다음과 같습니다.
방정식의 계수 값은 다음과 같습니다.
(대각 값만 해당)
이제 방정식을 코딩하는 함수를 만들 수 있습니다. 이 함수에는 시그니처 [c,f,s] = pdefun(x,t,u,dudx)
가 있어야 합니다.
x
는 독립적 공간 변수입니다.t
는 독립적 시간 변수입니다.u
는x
와t
에 따라 미분되는 종속 변수입니다. 요소를 2개 가진 벡터로, 여기서u(1)
은 이고u(2)
는 입니다.dudx
는 편미분 공간 도함수 입니다. 요소를 2개 가진 벡터로, 여기서dudx(1)
은 이고dudx(2)
는 입니다.출력값
c
,f
및s
는pdepe
에서 요구하는 표준 PDE 방정식 형식의 계수에 대응합니다.
그 결과, 이 예제의 방정식은 다음과 같은 함수로 표현할 수 있습니다.
function [c,f,s] = pdefun(x,t,u,dudx) c = [1; 1]; f = [0.024; 0.17] .* dudx; y = u(1) - u(2); F = exp(5.73*y)-exp(-11.47*y); s = [-F; F]; end
(참고: 모든 함수는 예제 끝에 로컬 함수로 포함되어 있습니다.)
초기 조건 코딩하기
이제 초기 조건을 반환하는 함수를 작성합니다. 초기 조건은 첫 번째 시간 값에서 적용되고 모든 x 값에 대한 값을 제공합니다. 초기 조건 개수는 방정식 개수와 같아야 하기 때문에 이 문제의 경우 초기 조건은 2개입니다. 함수 시그니처 u0 = pdeic(x)
를 사용하여 함수를 작성하십시오.
초기 조건은 다음과 같습니다.
대응하는 함수는 다음과 같습니다.
function u0 = pdeic(x) u0 = [1; 0]; end
경계 조건 코딩하기
이제 경계 조건을 계산하는 함수를 작성합니다.
구간 에 있는 문제의 경우, 경계 조건은 모든 , 그리고 또는 에 적용됩니다. 솔버에서 요구하는 경계 조건의 표준 형식은 다음과 같습니다.
이 형식으로 작성된 경우 의 편도함수에 대한 경계 조건은 플럭스 에 대해 표현해야 합니다. 따라서 이 문제의 경계 조건은 다음과 같습니다.
인 경우, 방정식은 다음과 같습니다.
계수는 다음과 같습니다.
마찬가지로 인 경우 방정식은 다음과 같습니다.
계수는 다음과 같습니다.
경계 함수는 함수 시그니처 [pl,ql,pr,qr] = pdebc(xl,ul,xr,ur,t)
를 사용해야 합니다.
입력값
xl
과ul
은 왼쪽 경계의 와 에 대응합니다.입력값
xr
과ur
은 오른쪽 경계의 와 에 대응합니다.t
는 독립적 시간 변수입니다.출력값
pl
과ql
은 왼쪽 경계(이 문제의 경우 )의 와 에 대응합니다.출력값
pr
과qr
은 오른쪽 경계(이 문제의 경우 )의 와 에 대응합니다.
이 예제의 경계 조건은 다음과 같은 함수로 표현됩니다.
function [pl,ql,pr,qr] = pdebc(xl,ul,xr,ur,t) pl = [0; ul(2)]; ql = [1; 0]; pr = [ur(1)-1; 0]; qr = [0; 1]; end
해 메시 선택하기
가 작을 때 이 문제의 해는 급격하게 변합니다. pdepe
가 급격한 변화에 대응하기에 적합한 시간 스텝을 선택하더라도 출력 플롯에서 거동을 확인하려면 적합한 출력 시간을 선택해야 합니다. 공간 메시의 경우, 이 해에는 의 양쪽 끝에 경계층이 있으므로 이러한 급격한 변화에 대응하기 위해 메시 점을 지정해야 합니다.
x = [0 0.005 0.01 0.05 0.1 0.2 0.5 0.7 0.9 0.95 0.99 0.995 1]; t = [0 0.005 0.01 0.05 0.1 0.5 1 1.5 2];
방정식 풀기
마지막으로 대칭 , PDE 방정식, 초기 조건, 경계 조건, 와 에 대한 메시를 사용하여 방정식을 풉니다.
m = 0; sol = pdepe(m,@pdefun,@pdeic,@pdebc,x,t);
pdepe
는 3차원 배열 sol
로 해를 반환합니다. 여기서 sol(i,j,k)
는 t(i)
와 x(j)
에서 구한 해의 k
번째 성분 에 대한 근삿값을 계산합니다. 각 해 성분을 별도의 변수로 추출하십시오.
u1 = sol(:,:,1); u2 = sol(:,:,2);
해 플로팅하기
와 에 대해 선택한 메시 점에서 플로팅된 과 에 대한 해의 곡면 플롯을 만듭니다.
surf(x,t,u1) title('u_1(x,t)') xlabel('Distance x') ylabel('Time t')
surf(x,t,u2) title('u_2(x,t)') xlabel('Distance x') ylabel('Time t')
로컬 함수(Local Function)
여기 나열된 함수는 PDE 솔버 pdepe
가 해를 계산하기 위해 호출하는 로컬 헬퍼 함수입니다. 또는 이러한 함수를 MATLAB 경로에 있는 디렉터리에 고유의 파일로 저장할 수도 있습니다.
function [c,f,s] = pdefun(x,t,u,dudx) % Equation to solve c = [1; 1]; f = [0.024; 0.17] .* dudx; y = u(1) - u(2); F = exp(5.73*y)-exp(-11.47*y); s = [-F; F]; end % --------------------------------------------- function u0 = pdeic(x) % Initial Conditions u0 = [1; 0]; end % --------------------------------------------- function [pl,ql,pr,qr] = pdebc(xl,ul,xr,ur,t) % Boundary Conditions pl = [0; ul(2)]; ql = [1; 0]; pr = [ur(1)-1; 0]; qr = [0; 1]; end % ---------------------------------------------