특이 항을 갖는 BVP 풀기

이 예제에서는 엠덴 방정식을 푸는 방법을 보여줍니다. 엠덴 방정식은 가스의 구면체를 모델링하는 과정에서 발생하는 특이 항을 갖는 경계값 문제입니다.

대칭을 사용하여 모델의 PDE를 줄이면 방정식이 다음과 같이 구간 $[0, 1]$ 에서 정의된 2계 ODE가 됩니다.

$y^{''} + \frac{2}{x} y^{'} + y^{5} = 0$ .

$x = 0$ 에서 $(2 / x)$ 항은 특이값을 갖지만 대칭이 경계 조건 $y^{'} (0) = 0$ 을 의미합니다. 이 경계 조건에서 $(2 / x) y^{'}$ 항은 $x \to 0$ 일 때 잘 정의되어 있습니다. 경계 조건 $y (1) = \sqrt{3} / 2$ 의 경우 BVP에는 MATLAB®에서 계산된 수치 해와 비교할 수 있는 다음 해석적 해가 있습니다.

$y (x) = {[\sqrt{1 + \frac{x^{2}}{3}}]}^{- 1}$ .

BVP 솔버 bvp4c는 다음 형식을 가진 특이 BVP를 풀 수 있습니다.

$y^{'} = S \frac{y}{x} + f (x, y)$ .

행렬 $S$ 는 상수여야 하며 $x = 0$ 에서 경계 조건은 필요 조건 $S \cdot y (0) = 0$ 과 일관되어야 합니다. bvpset의 'SingularTerm' 옵션을 사용하여 S 행렬을 솔버로 전달합니다.

$y_{1} = y$ 와 $y_{2} = y^{'}$ 를 사용하여 다음과 같이 1계 연립방정식으로 엠덴 방정식을 다시 작성할 수 있습니다.

${y_{1}}^{'} = y_{2}$ ,

${y_{2}}^{'} = - \frac{2}{x} y_{2} - y_{1}^{5}$ .

경계 조건이 다음과 같이 됩니다.

$y_{2} (0) = 0$ ,

$y_{1} (1) = \sqrt{3} / 2$ .

필요한 행렬 형식에서 시스템은 다음과 같습니다.

$[\begin{matrix} {y_{1}}^{'} \\ {y_{2}}^{'} \end{matrix}] = \frac{1}{x} [\begin{matrix} 0 & 0 \\ 0 & - 2 \end{matrix}] [\begin{matrix} y_{1} \\ y_{2} \end{matrix}] + [\begin{matrix} y_{2} \\ - y_{1}^{5} \end{matrix}]$ .

행렬 형식에서 다음은 분명합니다. $S = [\begin{matrix} 0 & 0 \\ 0 & - 2 \end{matrix}]$ 및 $f (x, y) = [\begin{matrix} y_{2} \\ - y_{1}^{5} \end{matrix}]$

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

방정식 코딩하기

$f (x, y)$ 에 대한 방정식을 코딩하는 함수를 만듭니다. 이 함수는 시그니처 dydx = emdenode(x,y)를 사용해야 합니다. 여기서,

x는 독립 변수입니다.
y는 종속 변수입니다.
dydx(1)은 ${y_{1}}^{'}$ 에 대한 방정식을 제공하고 dydx(2)는 ${y_{2}}^{'}$ 에 대한 방정식을 제공합니다.

이 입력값은 솔버에 의해 자동으로 함수로 전달되지만 변수 이름에 따라 방정식 코딩 방식이 달라집니다. 이 경우:

function dydx = emdenode(x,y)
dydx = [y(2) 
       -y(1)^5];
end

S를 포함한 항이 옵션을 사용하여 솔버로 전달되므로 함수에 해당 항은 포함되지 않습니다.

경계 조건 코딩하기

이제 경계점에서 경계 조건의 잔차 값을 반환하는 함수를 작성합니다. 이 함수는 시그니처 res = emdenbc(ya,yb)를 사용해야 합니다. 여기서,

ya는 구간의 시작 부분에서의 경계 조건 값입니다.
yb는 구간의 끝부분에서의 경계 조건 값입니다.

이 문제의 경우 경계 조건 중 하나는 $y_{1}$ 에 대한 것이고 다른 하나는 $y_{2}$ 에 대한 것입니다. 잔차 값을 계산하려면 경계 조건을 $g (x, y) = 0$ 형식으로 넣어야 합니다.

이 형식에서 경계 조건은 다음과 같습니다.

$y_{2} (0) = 0$ ,

$y_{1} (1) - \sqrt{3} / 2 = 0$ .

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

function res = emdenbc(ya,yb)
res = [ya(2)
       yb(1) - sqrt(3)/2];
end

초기 추측값 만들기

마지막으로, 해의 초기 추측값을 만듭니다. 이 문제에서는 경계 조건을 충족하는 상수 초기 추측값과 0~1 사이의 5개 점으로 구성된 간단한 메시를 사용합니다. BVP 솔버는 풀이 과정에서 이 점들을 변경해 활용하기 때문에 많은 메시 점을 사용할 필요가 없습니다.

$y_{1} = \sqrt{3} / 2$ ,

$y_{2} = 0$ .

guess = [sqrt(3)/2; 0];
xmesh = linspace(0,1,5);
solinit = bvpinit(xmesh, guess);

방정식 풀기

S에 대한 행렬을 만들어 'SingularTerm' 옵션의 값으로 bvpset에 전달합니다. 마지막으로 ODE 함수, 경계 조건 함수, 초기 추측값 및 options 구조체를 사용하여 bvp4c를 호출합니다.

S = [0 0; 0 -2];
options = bvpset('SingularTerm',S);
sol = bvp4c(@emdenode, @emdenbc, solinit, options);

해 플로팅하기

$[0, 1]$ 에서 해석적 해의 값을 계산합니다.

x = linspace(0,1);
truy = 1 ./ sqrt(1 + (x.^2)/3);

비교를 위해 bvp4c에서 계산된 해와 해석적 해를 플로팅합니다.

plot(x,truy,sol.x,sol.y(1,:),'ro');
title('Emden Problem -- BVP with Singular Term.')
legend('Analytical','Computed');
xlabel('x');
ylabel('solution y');

Figure contains an axes object. The axes object with title Emden Problem -- BVP with Singular Term., xlabel x, ylabel solution y contains 2 objects of type line. One or more of the lines displays its values using only markers These objects represent Analytical, Computed.

로컬 함수(Local Function)

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

function dydx = emdenode(x,y) % equation being solved 
dydx = [y(2) 
       -y(1)^5];
end
%-------------------------------------------
function res = emdenbc(ya,yb) % boundary conditions
res = [ya(2)
       yb(1) - sqrt(3)/2];
end
%-------------------------------------------

참고 항목

bvp4c | bvp5c | bvpinit | bvpset