연속법(Continuation)을 사용하여 BVP 풀기

이 예제에서는 문제를 보다 간단한 일련의 여러 문제로 효과적으로 분해하는 연속법(Continuation)을 사용하여 수치적으로 어려운 경계값 문제를 푸는 방법을 보여줍니다.

$0<\mathit{e}\ll 1$에 대해 다음과 같은 미분 방정식이 있다고 가정해 보겠습니다.

$\mathit{e}{{\mathit{y}}^{\prime }}^{\prime }+{\mathrm{xy}}^{\prime }=-\mathit{e}{\pi }^2\cos \left(\pi \mathit{x}\right)-\pi \mathrm{xsin}\left(\pi \mathit{x}\right)$.

이 문제는 구간 $\left[-1,1\right]$에 있으며 다음 경계 조건이 적용됩니다.

$\mathit{y}\left(-1\right)=-2$,

$\mathit{y}\left(1\right)=0$.

$\mathit{e}={10}^{-4}$일 때 방정식의 해가 $\mathit{x}=0$ 근처에서 급격히 천이하므로 수치적으로 풀기가 어렵습니다. 대신 이 예제에서는 연속법을 사용하여 $\mathit{e}$의 여러 값을 $\mathit{e}={10}^{-4}$이 될 때까지 반복해서 적용합니다. 중간 해는 각각 다음 문제에 대한 초기 추측값으로 사용됩니다.

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

방정식 코딩하기

대입 ${\mathit{y}}_1=\mathit{y}$ 및 ${\mathit{y}}_2={\mathit{y}}^{\prime }$를 사용하여 방정식을 다음과 같은 1계 연립방정식으로 다시 작성할 수 있습니다.

${{\mathit{y}}_1}^{\prime }={\mathit{y}}_2$,

${{\mathit{y}}_2}^{\prime }=-\frac{\mathit{x}}{\mathit{e}}{\mathit{y}}^{\prime }-{\pi }^2\cos \left(\pi \mathit{x}\right)-\frac{\pi \mathit{x}}{\mathit{e}}\sin \left(\pi \mathit{x}\right)$.

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

• x는 독립 변수입니다.

• y는 종속 변수입니다.

• dydx(1)은 ${{\mathit{y}}_1}^{\prime }$에 대한 방정식을 제공하고 dydx(2)는 ${{\mathit{y}}_2}^{\prime }$에 대한 방정식을 제공합니다.

shockode([x1 x2 ...],[y1 y2 ...])가 [shockode(x1,y1) shockode(x2,y2) ...]를 반환하도록 함수를 벡터화합니다. 이 방식은 솔버 성능을 향상시킵니다.

대응하는 함수는 다음과 같습니다.

function dydx = shockode(x,y)
pix = pi*x;
dydx = [y(2,:)
       -x/e.*y(2,:) - pi^2*cos(pix) - pix/e.*sin(pix)];
end

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

경계 조건 코딩하기

BVP 솔버를 사용하려면 경계 조건이 $\mathit{g}\left(\mathit{y}\left(\mathit{a}\right),\mathit{y}\left(\mathit{b}\right)\right)=0$ 형식이어야 합니다. 이 형식에서 경계 조건은 다음과 같습니다.

$\mathit{y}\left(-1\right)+2=0$,

$\mathit{y}\left(1\right)=0$.

경계 조건을 코딩하기 위해 시그니처 res = shockbc(ya,yb)가 있는 함수를 작성합니다. 여기서,

• ya는 구간 $\left[\mathit{a},\mathit{b}\right]$에서 시작 부분의 경계 조건 값입니다.

• yb는 구간 $\left[\mathit{a},\mathit{b}\right]$에서 끝부분의 경계 조건 값입니다.

대응하는 함수는 다음과 같습니다.

function res = shockbc(ya,yb) % boundary conditions
res = [ya(1)+2
       yb(1)];
end

야코비 행렬 코딩하기

ODE 함수 및 경계 조건에 대한 해석적 야코비 행렬은 이 문제에서 쉽게 계산할 수 있습니다. 야코비 행렬을 제공하면 더 이상 솔버가 유한 차분으로 근사치를 계산할 필요가 없으므로 솔버의 효율성이 높아집니다.

ODE 함수의 경우 야코비 행렬은 다음과 같습니다.

${\mathit{J}}_{\mathrm{ODE}}=\frac{\partial \mathit{f}}{\partial \mathit{y}}=\left[\begin{array}{cc}\frac{\partial {\mathit{f}}_1}{\partial {\mathit{y}}_1}& \frac{\partial {\mathit{f}}_1}{\partial {\mathit{y}}_2}\\ \frac{\partial {\mathit{f}}_2}{\partial {\mathit{y}}_1}& \frac{\partial {\mathit{f}}_2}{\partial {\mathit{y}}_2}\end{array}\right]=\left[\begin{array}{cc}0& 1\\ 0& -\frac{\mathit{x}}{\mathit{e}}\end{array}\right]$.

대응하는 함수는 다음과 같습니다.

function jac = shockjac(x,y,e)
jac = [0   1
       0  -x/e];
end

마찬가지로, 경계 조건의 경우 야코비 행렬은 다음과 같습니다.

${\mathit{J}}_{\mathit{y}\left(\mathit{a}\right)}=\left[\begin{array}{cc}1& 0\\ 0& 0\end{array}\right]$ , ${\mathit{J}}_{\mathit{y}\left(\mathit{b}\right)}=\left[\begin{array}{cc}0& 0\\ 1& 0\end{array}\right]$.

대응하는 함수는 다음과 같습니다.

function [dBCdya,dBCdyb] = shockbcjac(ya,yb)
dBCdya = [1 0; 0 0];
dBCdyb = [0 0; 1 0];
end

초기 추측값 생성하기

상수 추측값을 사용하여 $\left[-1,1\right]$의 5개 점으로 구성된 메시에 대한 해를 구합니다.

sol = bvpinit([-1 -0.5 0 0.5 1],[1 0]);

방정식 풀기

$\mathit{e}={10}^{-4}$를 사용하여 방정식을 직접 풀려는 경우 솔버가 $\mathit{x}=0$ 전환점 근처에서 문제의 좋지 않은 조건을 극복할 수 없습니다. 대신 $\mathit{e}={10}^{-4}$에 대한 해를 구하기 위해 이 예제에서는 연속법을 사용하여 ${10}^{-2}$, ${10}^{-3}$ 및 ${10}^{-4}$에 대한 일련의 문제를 풉니다. 각 반복에서 솔버의 출력값은 다음 반복에서 해에 대한 추측값으로 작용합니다(따라서 bvpinit에서 초기 추측값에 대한 변수가 sol이며 솔버의 출력값도 sol이라는 이름을 갖게 됩니다).

야코비 행렬의 값은 $\mathit{e}$의 값에 따라 다르므로 야코비 행렬에 대한 shockjac 함수와 shockbcjac 함수를 지정하는 루프에 옵션을 설정합니다. 또한 shockode는 값 벡터를 처리하도록 코딩되어 있으므로 벡터화를 활성화합니다.

e = 0.1;
for i = 2:4
   e = e/10;
   options = bvpset('FJacobian',@(x,y) shockjac(x,y,e),...
                    'BCJacobian',@shockbcjac,...
                    'Vectorized','on');
   sol = bvp4c(@(x,y) shockode(x,y,e),@shockbc, sol, options);
end

해 플로팅하기

메시 $\mathit{x}$와 해 $\mathit{y}\left(\mathit{x}\right)$에 대한 bvp4c의 결과를 플로팅합니다. 솔버는 연속법을 사용하여 $\mathit{x}=0$에서의 불연속을 처리할 수 있습니다.

plot(sol.x,sol.y(1,:),'-o');
axis([-1 1 -2.2 2.2]);
title(['There Is a Shock at x = 0 When e =' sprintf('%.e',e) '.']);
xlabel('x');
ylabel('solution y');

로컬 함수(Local Function)

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

function dydx = shockode(x,y,e) % equation to solve
pix = pi*x;
dydx = [y(2,:)
       -x/e.*y(2,:) - pi^2*cos(pix) - pix/e.*sin(pix)];
end
%-------------------------------------------
function res = shockbc(ya,yb) % boundary conditions
res = [ya(1)+2
       yb(1)];
end
%-------------------------------------------
function jac = shockjac(x,y,e) % jacobian of shockode
jac = [0   1
       0  -x/e];
end
%-------------------------------------------
function [dBCdya,dBCdyb] = shockbcjac(ya,yb) % jacobian of shockbc
dBCdya = [1 0; 0 0];
dBCdyb = [0 0; 1 0];
end
%-------------------------------------------

참고 항목

bvp4c | bvpinit | bvpset

관련 항목

• 경계값 문제 풀기
• 연속법(Continuation)을 사용하여 BVP 일관성 확인하기