bvp4c

경계값 문제 풀기 — 4계 방법

구문

sol = bvp4c(odefun,bcfun,solinit)

sol = bvp4c(odefun,bcfun,solinit,options)

설명

sol = bvp4c(odefun,bcfun,solinit)는 odefun으로 지정된 y′ = f(x,y) 형식의 연립미분방정식을 적분합니다. bcfun은 경계 조건이며 solinit는 초기 해 추측값입니다. bvpinit 함수를 사용하여 초기 추측값 solinit를 만듭니다. 초기 추측값은 bcfun의 경계 조건이 적용되는 점도 정의합니다.

sol = bvp4c(odefun,bcfun,solinit,options)도 options에 의해 정의된 적분 설정을 사용합니다. 적분 설정은 bvpset 함수를 사용하여 만든 인수입니다. 예를 들어, 절대 허용오차와 상대 허용오차를 지정하려면 AbsTol 옵션과 RelTol 옵션을 사용하고 odefun의 해석적 편도함수를 제공하려면 FJacobian 옵션을 사용하십시오.

예제

2계 BVP 풀기

MATLAB®에서 함수를 사용하여 2계 BVP를 풉니다. 이 예제에서는 다음과 같은 2계 방정식을 사용합니다.

${{\mathit{y}}^{\prime }}^{\prime }+\mathit{y}=0$.

이 방정식은 구간 $\left[0,\pi /2\right]$에 대해 정의되며 다음 경계 조건이 적용됩니다.

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

$\mathit{y}\left(\pi /2\right)=2$.

MATLAB에서 이 방정식을 풀려면 이 방정식을 1계 연립방정식으로 표현하는 함수, 경계 조건에 대한 함수, 그리고 초기 추측값에 대한 함수를 작성해야 합니다. 그러면 BVP 솔버가 3개의 입력값을 사용하여 방정식을 풉니다.

방정식 코딩하기

방정식을 코딩하는 함수를 작성합니다. 대입 ${\mathit{y}}_1=\mathit{y}$와 ${\mathit{y}}_2={\mathit{y}}^{\prime }$을 사용하여 방정식을 1계 연립방정식으로 바꿉니다.

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

${{\mathit{y}}_2}^{\prime }=-{\mathit{y}}_1$.

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

function dydx = bvpfcn(x,y)
dydx = zeros(2,1);
dydx = [y(2)
       -y(1)];
end

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

경계 조건 코딩하기

경계 조건을 $\mathit{g}\left(\mathit{y}\left(\mathit{a}\right),\mathit{y}\left(\mathit{b}\right)\right)=0$ 형식으로 코딩하는 함수를 작성합니다. 이 형식에서 경계 조건은 다음과 같습니다.

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

$\mathit{y}\left(\pi /2\right)-2=0$.

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

function res = bcfcn(ya,yb)
res = [ya(1)
       yb(1)-2];
end

초기 추측값 만들기

bvpinit 함수를 사용하여 방정식 해의 초기 추측값을 만듭니다. 방정식은 ${{\mathit{y}}^{\prime }}^{\prime }$을 $\mathit{y}$에 연결시키므로 해가 삼각 함수와 관련이 있다고 합리적으로 추측할 수 있습니다. 적분 구간에서 5개 점으로 구성된 메시를 사용합니다. 솔버는 메시의 첫 번째 값과 마지막 값에 경계 조건을 적용합니다.

초기 추측값에 대한 함수는 $\mathit{x}$를 입력값으로 받고 ${\mathit{y}}_1$과 ${\mathit{y}}_2$의 값에 대한 추측값을 반환합니다. 함수는 다음과 같습니다.

function g = guess(x)
g = [sin(x)
     cos(x)];
end

xmesh = linspace(0,pi/2,5);
solinit = bvpinit(xmesh, @guess);

방정식 풀기

bvp4c를 도함수, 경계 조건 함수 및 초기 추측값과 함께 사용하여 문제를 풉니다.

sol = bvp4c(@bvpfcn, @bcfcn, solinit);

해 플로팅하기

plot(sol.x, sol.y, '-o')

로컬 함수(Local Function)

다음은 bvp4c가 방정식을 푸는 데 사용하는 로컬 함수입니다.

function dydx = bvpfcn(x,y) % equation to solve
dydx = zeros(2,1);
dydx = [y(2)
       -y(1)];
end
%--------------------------------
function res = bcfcn(ya,yb) % boundary conditions
res = [ya(1)
       yb(1)-2];
end
%--------------------------------
function g = guess(x) % initial guess for y and y'
g = [sin(x)
     cos(x)];
end
%--------------------------------

bvp4c와 bvp5c 솔버 비교하기

2개의 다른 솔버를 사용하여 엄격하지 않은 허용오차로 BVP를 풀고 결과를 비교합니다.

다음과 같은 2계 ODE가 있다고 가정하겠습니다.

${{\mathit{y}}^{\prime }}^{\prime }+\frac{2}{\mathit{x}}{\mathit{y}}^{\prime }+\frac{1}{{\mathit{x}}^4}\mathit{y}=0$.

이 방정식은 구간 $\left[\frac{1}{3\pi },1\right]$에 대해 정의되며 다음 경계 조건이 적용됩니다.

$\mathit{y}\left(\frac{1}{3\pi }\right)=0$,

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

MATLAB®에서 이 방정식을 풀려면 이 방정식을 1계 연립방정식으로 표현하는 함수와 경계 조건에 대한 함수를 작성하고 옵션 값을 설정하고 초기 추측값을 만들어야 합니다. 그러면 BVP 솔버가 4개의 입력값을 사용하여 방정식을 풉니다.

방정식 코딩하기

대입 ${\mathit{y}}_1=\mathit{y}$와 ${\mathit{y}}_2={\mathit{y}}^{\prime }$을 사용하여 ODE를 1계 연립방정식으로 바꿀 수 있습니다.

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

${{\mathit{y}}_2}^{\prime }=-\frac{2}{\mathit{x}}{\mathit{y}}_2-\frac{1}{{\mathit{x}}^4}{\mathit{y}}_1$.

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

function dydx = bvpfcn(x,y)
dydx = [y(2)
       -2*y(2)/x - y(1)/x^4];
end

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

경계 조건 코딩하기

경계 조건 함수에서는 경계 조건이 $\mathit{g}\left(\mathit{y}\left(\mathit{a}\right),\mathit{y}\left(\mathit{b}\right)\right)=0$과 같은 형식을 가져야 합니다. 이 형식에서 경계 조건은 다음과 같습니다.

$\mathit{y}\left(\frac{1}{3\pi }\right)=0$,

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

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

function res = bcfcn(ya,yb)
res = [ya(1)
       yb(1)-sin(1)];
end

옵션 설정하기

두 솔버 사이의 오차 제어 차이점을 보기 위해 bvpset를 사용하여 솔버 통계량 표시를 설정하고 엄격하지 않은 허용오차를 지정합니다. 효율을 높이기 위해 다음과 같은 해석적 야코비 행렬을 지정합니다.

$\mathit{J}=\frac{\partial {\mathit{f}}_{\mathit{i}}}{\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\\ -\frac{1}{{\mathit{x}}^4}& -\frac{2}{\mathit{x}}\end{array}\right]$.

야코비 행렬의 값을 반환하는 함수는 다음과 같습니다.

function dfdy = jac(x,y)
dfdy = [0      1
       -1/x^4 -2/x];
end

opts = bvpset('FJacobian',@jac,'RelTol',0.1,'AbsTol',0.1,'Stats','on');

초기 추측값 만들기

bvpinit를 사용하여 해의 초기 추측값을 만듭니다. 구간 $\left[1/3\pi ,1\right]$에 있는 10개 점으로 구성된 초기 메시를 사용하여 상수 함수를 초기 추측값으로 지정합니다.

xmesh = linspace(1/(3*pi), 1, 10);
solinit = bvpinit(xmesh, [1; 1]);

방정식 풀기

bvp4c와 bvp5c를 모두 사용하여 방정식을 풉니다.

sol4c = bvp4c(@bvpfcn, @bcfcn, solinit, opts);

The solution was obtained on a mesh of 9 points.
The maximum residual is  9.794e-02. 
There were 157 calls to the ODE function. 
There were 28 calls to the BC function. 

sol5c = bvp5c(@bvpfcn, @bcfcn, solinit, opts);

The solution was obtained on a mesh of 11 points.
The maximum error is  6.742e-02. 
There were 244 calls to the ODE function. 
There were 29 calls to the BC function. 

결과 플로팅하기

비교를 위해 ${\mathit{y}}_1$에 대한 두 가지 계산의 결과를 해석적 해와 함께 플로팅합니다. 해석적 해는 다음과 같습니다.

${\mathit{y}}_1=\sin \left(\frac{1}{\mathit{x}}\right)$,

${\mathit{y}}_2=-\frac{1}{{\mathit{x}}^2}\cos \left(\frac{1}{\mathit{x}}\right)$.

xplot = linspace(1/(3*pi),1,200);
yplot = [sin(1./xplot); -cos(1./xplot)./xplot.^2];

plot(xplot,yplot(1,:),'k',sol4c.x,sol4c.y(1,:),'r*',sol5c.x,sol5c.y(1,:),'bo')
title('Comparison of BVP Solvers with Crude Error Tolerance')
legend('True','BVP4C','BVP5C')
xlabel('x')
ylabel('solution y')

플롯을 통해 bvp5c는 계산에서 정확한 오차를 직접 제어하지만 bvp4c는 정확한 오차를 간접적으로만 제어한다는 사실을 확인할 수 있습니다. 허용오차가 보다 엄격해지면 두 솔버 사이의 차이점이 이렇게까지 뚜렷하지는 않아집니다.

로컬 함수(Local Function)

다음은 BVP 솔버가 문제를 푸는 데 사용하는 로컬 함수입니다.

function dydx = bvpfcn(x,y) % equation to solve
dydx = [y(2)
       -2*y(2)/x - y(1)/x^4];
end
%---------------------------------
function res = bcfcn(ya,yb) % boundary conditions
res = [ya(1)
       yb(1)-sin(1)];
end
%---------------------------------
function dfdy = jac(x,y) % analytical jacobian for f
dfdy = [0      1
       -1/x^4 -2/x];
end
%---------------------------------

입력 인수

odefun — 풀려는 함수

함수 핸들

풀려는 함수로, 적분할 함수를 정의하는 함수 핸들로 지정됩니다. odefun과 bcfun은 동일한 개수의 입력 인수를 받아야 합니다.

odefun을 코딩하려면 스칼라 x와 열 벡터 y에 대해 함수 시그니처 dydx = odefun(x,y)를 사용하십시오. 반환 값 dydt는 f(x,y)에 대응하는 single이나 double 데이터형의 열 벡터입니다. odefun은 입력 인수 x와 y 모두를 받아야만 하는데, 이는 이 중 하나가 함수에 사용되지 않더라도 마찬가지입니다.

예를 들어, $$y\text{'}=5y-3$$을 풀려면 다음 함수를 사용하십시오.

function dydt = odefun(t,y)
dydt = 5*y-3;
end

연립방정식에 대한 odefun의 출력값은 벡터입니다. 벡터의 각 요소는 하나의 방정식의 해입니다. 예를 들어 다음을 풀려면

$$\begin{array}{l}y{\text{'}}_1=y_1+2y_2\\ y{\text{'}}_2=3y_1+2y_2\end{array}$$

다음 함수를 사용하십시오.

function dydt = odefun(t,y)
dydt = zeros(2,1);
dydt(1) = y(1)+2*y(2);
dydt(2) = 3*y(1)+2*y(2);
end

bvp4c는 해에 특이점이 있거나 다중 점 경계 조건을 갖는 문제도 풀 수 있습니다.

예: sol = bvp4c(@odefun, @bcfun, solinit)

알 수 없는 파라미터

풀려는 BVP에 알 수 없는 파라미터가 포함된 경우 함수 시그니처 dydx = odefun(x,y,p)를 대신 사용할 수 있습니다. 여기서 p는 파라미터 값으로 구성된 벡터입니다. 이 함수 시그니처를 사용하면 BVP 솔버가 solinit에 제공된 파라미터 값에 대한 초기 추측값부터 시작하여 알 수 없는 파라미터의 값을 계산합니다.

데이터형: function_handle

bcfun — 경계 조건

함수 핸들

경계 조건으로, 경계 조건의 잔차 오차를 계산하는 함수 핸들로 지정됩니다. odefun과 bcfun은 동일한 개수의 입력 인수를 받아야 합니다.

bcfun을 코딩하려면 열 벡터 ya와 yb에 대해 함수 시그니처 res = bcfun(ya,yb)를 사용하십시오. 반환 값 res는 경계 점에서의 경계 조건의 잔차 값에 해당하는 single형이나 double형의 열 벡터입니다.

예를 들어, y(a) = 1이고 y(b) = 0인 경우 경계 조건 함수는 다음과 같습니다.

function res = bcfun(ya,yb)
res = [ya(1)-1
       yb(1)];
end

y(a) = 1이므로 점 x = a에서 ya(1)-1의 잔차 값은 0이 되어야 합니다. 마찬가지로, y(b) = 0이므로 점 x = b에서 yb(1)의 잔차 값은 0이 되어야 합니다.

경계 조건이 적용되는 경계 점 x = a와 x = b는 초기 추측값 구조체 solinit에서 정의됩니다. 2점 경계값 문제의 경우 a = solinit.x(1)이고 b = solinit.x(end)가 됩니다.

예: sol = bvp4c(@odefun, @bcfun, solinit)

알 수 없는 파라미터

풀려는 BVP에 알 수 없는 파라미터가 포함된 경우 함수 시그니처 res = bcfun(ya,yb,p)를 대신 사용할 수 있습니다. 여기서 p는 파라미터 값으로 구성된 벡터입니다. 이 함수 시그니처를 사용하면 BVP 솔버가 solinit에 제공된 파라미터 값에 대한 초기 추측값부터 시작하여 알 수 없는 파라미터의 값을 계산합니다.

데이터형: function_handle

solinit — 해의 초기 추측값

구조체

해의 초기 추측값으로, 구조체로 지정됩니다. bvpinit를 사용하여 solinit를 만듭니다.

초기값 문제와 달리, 경계값 문제는 해를 가지지 않거나, 유한개의 해를 가지거나, 무한개의 많은 해를 가질 수 있습니다. 필요한 해에 대한 추측값을 제공하는 것은 BVP를 푸는 과정의 중요한 부분입니다. 이 추측값의 품질은 솔버 성능은 물론, 계산 성공 여부에도 중요한 요소로 작용할 수 있습니다. 양호한 초기 추측값을 만드는 방법에 대한 지침은 해의 초기 추측값 항목을 참조하십시오.

예: sol = bvp4c(@odefun, @bcfun, solinit)

데이터형: struct

options — options 구조체

구조체

options 구조체입니다. bvpset 함수를 사용하여 options 구조체를 만들거나 수정할 수 있습니다.

예: options = bvpset('RelTol',1e-5,'Stats','on')은 1e-5의 상대 허용오차를 지정하고 솔버 통계량 표시를 켭니다.

데이터형: struct

출력 인수

sol — 해 구조체

구조체

해 구조체입니다. sol.field1과 같이 점 인덱싱을 사용하여 sol의 필드에 액세스할 수 있습니다. 해 (sol.x,sol.y)는 초기 메시 solinit.x에서 정의된 적분 구간에서 연속하고 거기에서 연속 1계 도함수를 갖습니다. sol을 deval 함수와 함께 사용하여 구간 내에 있는 다른 점에서 해를 계산할 수 있습니다.

구조체 sol은 다음 필드를 가집니다.

필드	설명
x	bvp4c에서 선택된 메시입니다. 이 메시는 일반적으로 초기 메시 solinit.x와 다른 점을 포함합니다.
y	sol.x의 메시 점에서의 y(x)에 대한 근삿값입니다.
yp	sol.x의 메시 점에서의 y′(x)에 대한 근삿값입니다.
parameters	solinit.parameters에 지정된 알 수 없는 파라미터에 대한 최종 값입니다.
solver	'bvp4c'
stats	해와 관련 있는 계산 비용 통계량, 즉 메시 점의 개수, 잔차 오차, 그리고 odefun과 bcfun에 대한 호출의 개수입니다.

세부 정보

다중 점 경계값 문제(Multipoint Boundary Value Problem)

다중 점 경계값 문제의 경우, 경계 조건이 적분 구간의 여러 점에 적용됩니다.

bvp4c는 a = a₀ < a₁ < a₂ < ...< a_n = b가 구간 [a,b]의 경계 점인, 다중 점 경계값 문제를 풀 수 있습니다. 점 a₁,a₂,...,a_n−1은 [a,b]를 여러 영역으로 나누는 접점을 나타냅니다. bvp4c는 1부터 시작하는 인덱스를 사용하여 왼쪽에서 오른쪽으로(a에서 b로) 영역을 나열합니다. 영역 k, [a_k−1,a_k]에서 bvp4c는 다음과 같이 도함수를 계산합니다.

yp = odefun(x,y,k) 

경계 조건 함수 bcfun(yleft,yright)에서 yleft(:,k)는 [a_k−1,a_k]의 왼쪽 경계에서의 해입니다. 마찬가지로 yright(:,k)는 k 영역의 오른쪽 경계에서의 해입니다. 특히, yleft(:,1) = y(a)이고 yright(:,end) = y(b)입니다.

bvpinit를 사용하여 초기 추측값을 만들 경우 각 접점의 xinit에 2개의 요소를 사용해야 합니다. 자세한 내용은 bvpinit에 대한 함수 도움말 페이지를 참조하십시오.

yinit가 함수인 경우 bvpinit는 y = yinit(x,k)를 호출하여 영역 k의 x에서 해의 초기 추측값을 가져옵니다. bpv4c에서 반환되는 해 구조체 sol에서 sol.x는 각 접점에 대해 2개의 항목을 가집니다. sol.y의 대응하는 열에는 접점의 왼쪽 해와 오른쪽 해가 각각 들어 있습니다.

3점 경계값 문제를 푸는 예제를 보려면 여러 경계 조건이 있는 BVP 풀기 항목을 참조하십시오.

특이 경계값 문제(Singular Boundary Value Problem)

bvp4c는 알 수 없는 파라미터 p가 있는 다음 형식의 문제를 비롯하여 특이 경계값 문제의 클래스를 풉니다.

$$\begin{array}{l}y\text{'}=S\frac{y}{x}+f\left(x,y,p\right),\\ 0=bc\left(y\left(0\right),y\left(b\right),p\right).\end{array}$$

구간은 [0, b]이어야 하고 b > 0이어야 합니다. 원기둥 또는 구의 대칭으로 인해 PDE(편미분 방정식)의 결과로 산출된 ODE의 매끄러운 해를 구할 때 이러한 문제가 자주 나타납니다. 특이 문제의 경우 (상수) 행렬 S를 bvpset의 'SingularTerm' 옵션 값으로 지정하면, odefun이 f(x,y,p)만 실행하게 됩니다. 경계 조건과 초기 추측값은 매끄러움을 위한 필요 조건 S·y(0) = 0과 모순되지 않아야 합니다.

특이 경계값 문제를 푸는 예제를 보려면 특이 항을 갖는 BVP 풀기 항목을 참조하십시오.

알고리즘

bvp4c는 3단계 로바토 IIIa(Three-stage Lobatto IIIa) 식을 구현하는 유한 미분 코드입니다 [1], [2]. 이는 선점(Collocation) 식으로, 적분 구간에서 일정한 4차 정확도를 가지며 C¹ 연속해를 제공하는 선점 다항식입니다. 메시 선택과 오차 제어는 연속해(Continuous Solution)의 잔차(Residual)를 기준으로 합니다.

선점 기법은 점들로 구성된 메시를 사용하여 적분 구간을 하위 구간으로 나눕니다. 솔버는 경계 조건에서 얻는 대수 방정식의 전역 시스템(Global System)을 풀어 수치 해를 구하고, 모든 하위 구간에 적용되는 선점 조건을 결정합니다. 그런 다음, 각 하위 구간에 대한 수치 해의 오차를 추정합니다. 해가 허용오차 조건을 충족하지 않을 경우 솔버는 이에 맞게 메시를 조정하고 프로세스를 반복합니다. 초기 메시의 점뿐만 아니라 메시 점에서 해의 초기 근삿값도 반드시 제공해야 합니다.

참고 문헌

[1] Shampine, L.F., and J. Kierzenka. "A BVP Solver based on residual control and the MATLAB PSE." ACM Trans. Math. Softw. Vol. 27, Number 3, 2001, pp. 299–316.

[2] Shampine, L.F., M.W. Reichelt, and J. Kierzenka. "Solving Boundary Value Problems for Ordinary Differential Equations in MATLAB with bvp4c." MATLAB File Exchange, 2004.

확장 기능

스레드 기반 환경

MATLAB®의 backgroundPool을 사용해 백그라운드에서 코드를 실행하거나 Parallel Computing Toolbox™의 ThreadPool을 사용해 코드 실행 속도를 높일 수 있습니다.

이 함수는 스레드 기반 환경을 완전히 지원합니다. 자세한 내용은 Run MATLAB Functions in Thread-Based Environment 항목을 참조하십시오.

버전 내역

R2006a 이전에 개발됨

참고 항목

bvp5c | bvpget | bvpinit | bvpset | bvpxtend | deval

도움말 항목

• 경계값 문제 풀기

×

예제 열기

이 예제의 수정된 버전이 있습니다. 사용자가 편집한 내용을 반영하여 이 예제를 여시겠습니까?

아니요, 수정된 버전을 덮어씁니다. 예