여러 경계 조건이 있는 BVP 풀기

이 예제에서는 다중 점 경계값 문제를 푸는 방법을 보여줍니다. 살펴보려는 해는 적분 구간 내에서 조건을 충족합니다.

$\left[0,\lambda \right]$에 있는 $\mathit{x}$에 대한 다음 방정식을 가정해 보겠습니다.

${\mathit{v}}^{\prime }=\frac{\mathit{C}-1}{\mathit{n}}$,

${\mathit{C}}^{\prime }=\frac{\mathrm{vC}-\min \left(\mathit{x},1\right)}{\eta }$.

문제의 알려진 파라미터는 $\mathit{n}$, $\kappa$, $\lambda >1$ 및 $\eta =\frac{{\lambda }^2}{\mathit{n}\cdot {\kappa }^2}$입니다.

${\mathit{C}}^{\prime }\left(\mathit{x}\right)$에 대한 방정식에 항 $\min \left(\mathit{x},1\right)$은 $\mathit{x}=1$에서 매끄럽지 않으므로 문제를 직접 풀 수 없습니다. 대신 문제를 둘로 나누어 하나는 $\left[0,1\right]$ 구간에서 풀고 다른 하나는 $\left[1,\lambda \right]$ 구간에서 풀 수 있습니다. 두 영역 간의 연결은 $\mathit{x}=1$에서 해가 연속적이어야 한다는 점을 기반으로 합니다. 또한 해는 경계 조건도 충족해야 합니다.

$\mathit{v}\left(0\right)=0$,

$\mathit{C}\left(\lambda \right)=1$.

각 영역에 대한 방정식은 다음과 같습니다.

영역 1: $0\le \mathit{x}\le 1$

${\mathit{v}}^{\prime }=\frac{\mathit{C}-1}{\mathit{n}}$,

${\mathit{C}}^{\prime }=\frac{\mathrm{vC}-\mathit{x}}{\eta }$.

영역 2: $1\le \mathit{x}\le \lambda$

${\mathit{v}}^{\prime }=\frac{\mathit{C}-1}{\mathit{n}}$,

${\mathit{C}}^{\prime }=\frac{\mathrm{vC}-1}{\eta }$.

접점 $\mathit{x}=1$은 두 영역에 모두 포함됩니다. 이 점에서 솔버는 왼쪽 해와 오른쪽 해를 모두 생성하며 해가 연속적이려면 두 해가 같아야 합니다.

MATLAB®에서 이 연립방정식을 풀려면 경계값 문제 솔버 bvp5c를 호출하기 전에 방정식, 경계 조건 및 초기 추측값을 코딩해야 합니다. 필요한 함수를 이 예제와 같이 파일 끝에 로컬 함수로 포함시킬 수도 있고, MATLAB 경로에 있는 디렉터리에 이름이 지정된 별도의 파일로 저장할 수도 있습니다.

방정식 코딩하기

${\mathit{v}}^{\prime }\left(\mathit{x}\right)$와 ${\mathit{C}}^{\prime }\left(\mathit{x}\right)$에 대한 방정식은 풀려는 영역에 따라 다릅니다. 다중 점 경계값 문제의 경우 도함수는 도함수를 계산할 영역을 식별하는 데 사용되는 세 번째 입력 인수 region을 받아야 합니다. 솔버는 왼쪽 영역에서 오른쪽 영역으로 1부터 번호를 지정합니다.

방정식을 코딩하기 위해 시그니처 dydx = f(x,y,region,p)가 있는 함수를 작성합니다. 여기서,

• x는 독립 변수입니다.

• y는 종속 변수입니다.

• dydx(1)은 ${\mathit{v}}^{\prime }\left(\mathit{x}\right)$에 대한 방정식을 제공하고 dydx(2)는 ${\mathit{C}}^{\prime }\left(\mathit{x}\right)$에 대한 방정식을 제공합니다.

• region은 도함수가 계산될 영역의 번호입니다. 이 두 영역 문제에서 region은 1 또는 2입니다.

• p는 상수 파라미터 $\left[\mathit{n},\kappa ,\lambda ,\eta \right]$의 값을 포함하는 벡터입니다.

풀려는 영역에 따라 다른 방정식이 반환되도록 switch 문을 사용합니다. 함수는 다음과 같습니다.

function dydx = f(x,y,region,p) % equations being solved
n = p(1);
eta = p(4);

dydx = zeros(2,1);
dydx(1) = (y(2) - 1)/n;

switch region
    case 1    % x in [0 1]
        dydx(2) = (y(1)*y(2) - x)/eta;
    case 2    % x in [1 lambda]
        dydx(2) = (y(1)*y(2) - 1)/eta;
end
end

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

경계 조건 코딩하기

두 영역에서 2개의 1계 미분 방정식을 푸는 경우 4개의 경계 조건이 필요합니다. 이러한 조건 중 다음의 두 조건은 원래 문제에서 가져옵니다.

$\mathit{v}\left(0\right)=0$,

$\mathit{C}\left(\lambda \right)-1=0$.

다음의 두 조건은 접점 $\mathit{x}=1$에서 왼쪽 해와 오른쪽 해가 연속되도록 합니다.

${\mathit{v}}_{\mathit{L}}\left(1\right)-{\mathit{v}}_{\mathit{R}}\left(1\right)=0$,

${\mathit{C}}_{\mathit{L}}\left(1\right)-{\mathit{C}}_{\mathit{R}}\left(1\right)=0$.

다중 점 BVP의 경우, 경계 조건 함수의 인수 YL과 YR은 행렬이 됩니다. 특히, k번째 열 YL(:,k)는 k번째 영역의 왼쪽 경계에서의 해입니다. 마찬가지로, YR(:,k)는 k번째 영역의 오른쪽 경계에서의 해입니다.

이 문제에서 $\mathit{y}\left(0\right)$의 근삿값은 YL(:,1)로 계산되고 $\mathit{y}\left(\lambda \right)$의 근삿값은 YR(:,end)로 계산됩니다. $\mathit{x}=1$에서 해가 연속되려면 YR(:,1) = YL(:,2)여야 합니다.

4개의 경계 조건의 잔차 값을 인코딩하는 함수는 다음과 같습니다.

function res = bc(YL,YR)
res = [YL(1,1)               % v(0) = 0
       YR(1,1) - YL(1,2)     % Continuity of v(x) at x=1
       YR(2,1) - YL(2,2)     % Continuity of C(x) at x=1
       YR(2,end) - 1];       % C(lambda) = 1
end

초기 추측값 생성하기

다중 점 BVP의 경우 경계 조건은 적분 구간의 시작 부분과 끝부분에 자동으로 적용됩니다. 그러나 나머지 접점에 대해서는 xmesh에서 두 개의 항목을 지정해야 합니다. 경계 조건을 충족하는 간단한 추측값은 상수 추측값 y = [1; 1]입니다.

xc = 1;
xmesh = [0 0.25 0.5 0.75 xc xc 1.25 1.5 1.75 2];
yinit = [1; 1];
sol = bvpinit(xmesh,yinit);

방정식 풀기

상수 파라미터의 값을 정의하여 벡터 p에 넣습니다. bvp5c에 대한 함수를 제공할 때 구문 @(x,y,r) f(x,y,r,p)를 사용하여 파라미터로 구성된 벡터를 제공합니다.

각 해를 다음의 초기 추측값으로 사용하여 $\kappa$의 여러 값에 대한 해를 계산합니다. 각 $\kappa$ 값에 대해 오스몰 농도 $\mathrm{Os}=\frac{1}{\mathit{v}\left(\lambda \right)}$의 값을 계산합니다. 루프의 매 반복에서 계산된 값을 해석적 해의 근삿값과 비교합니다.

lambda = 2;
n = 5e-2;
for kappa = 2:5
   eta = lambda^2/(n*kappa^2);
   p = [n kappa lambda eta];
   sol = bvp5c(@(x,y,r) f(x,y,r,p), @bc, sol);
   K2 = lambda*sinh(kappa/lambda)/(kappa*cosh(kappa));
   approx = 1/(1 - K2);
   computed = 1/sol.y(1,end);
   fprintf('  %2i    %10.3f     %10.3f \n',kappa,computed,approx);
end

   2         1.462          1.454 
   3         1.172          1.164 
   4         1.078          1.071 
   5         1.039          1.034 

해 플로팅하기

$\mathit{v}\left(\mathit{x}\right)$와 $\mathit{C}\left(\mathit{x}\right)$에 대한 해 성분과 접점 $\mathit{x}=1$에서의 세로선을 플로팅합니다. $\kappa =5$에 대해 표시된 해는 루프의 최종 반복에서 나온 결과입니다.

plot(sol.x,sol.y(1,:),'--o',sol.x,sol.y(2,:),'--o')
line([1 1], [0 2], 'Color', 'k')
legend('v(x)','C(x)')
title('A Three-Point BVP Solved with bvp5c')
xlabel({'x', '\lambda = 2, \kappa = 5'})
ylabel('v(x) and C(x)')

로컬 함수(Local Function)

BVP 솔버 bvp5c가 해를 계산하기 위해 호출하는 로컬 헬퍼 함수는 다음과 같습니다. 또는 이러한 함수를 MATLAB 경로에 있는 디렉터리에 고유의 파일로 저장할 수도 있습니다.

function dydx = f(x,y,region,p) % equations being solved
n = p(1);
eta = p(4);

dydx = zeros(2,1);
dydx(1) = (y(2) - 1)/n;

switch region
    case 1    % x in [0 1]
        dydx(2) = (y(1)*y(2) - x)/eta;
    case 2    % x in [1 lambda]
        dydx(2) = (y(1)*y(2) - 1)/eta;
end
end
%-------------------------------------------
function res = bc(YL,YR) % boundary conditions
res = [YL(1,1)               % v(0) = 0
       YR(1,1) - YL(1,2)     % Continuity of v(x) at x=1
       YR(2,1) - YL(2,2)     % Continuity of C(x) at x=1
       YR(2,end) - 1];       % C(lambda) = 1
end
%-------------------------------------------

참고 항목

bvp4c | bvp5c | bvpinit

관련 항목

• 경계값 문제 풀기