bvp4c
경계값 문제 풀기 — 4계 방법
설명
예제
2계 BVP 풀기
MATLAB®에서 함수를 사용하여 2계 BVP를 풉니다. 이 예제에서는 다음과 같은 2계 방정식을 사용합니다.
.
이 방정식은 구간 에 대해 정의되며 다음 경계 조건이 적용됩니다.
,
.
MATLAB에서 이 방정식을 풀려면 이 방정식을 1계 연립방정식으로 표현하는 함수, 경계 조건에 대한 함수, 그리고 초기 추측값에 대한 함수를 작성해야 합니다. 그러면 BVP 솔버가 3개의 입력값을 사용하여 방정식을 풉니다.
방정식 코딩하기
방정식을 코딩하는 함수를 작성합니다. 대입 와 을 사용하여 방정식을 1계 연립방정식으로 바꿉니다.
,
.
대응하는 함수는 다음과 같습니다.
function dydx = bvpfcn(x,y) dydx = zeros(2,1); dydx = [y(2) -y(1)]; end
참고: 모든 함수는 예제 끝에 로컬 함수로 포함되어 있습니다.
경계 조건 코딩하기
경계 조건을 형식으로 코딩하는 함수를 작성합니다. 이 형식에서 경계 조건은 다음과 같습니다.
,
.
대응하는 함수는 다음과 같습니다.
function res = bcfcn(ya,yb) res = [ya(1) yb(1)-2]; end
초기 추측값 만들기
bvpinit
함수를 사용하여 방정식 해의 초기 추측값을 만듭니다. 방정식은 을 에 연결시키므로 해가 삼각 함수와 관련이 있다고 합리적으로 추측할 수 있습니다. 적분 구간에서 5개 점으로 구성된 메시를 사용합니다. 솔버는 메시의 첫 번째 값과 마지막 값에 경계 조건을 적용합니다.
초기 추측값에 대한 함수는 를 입력값으로 받고 과 의 값에 대한 추측값을 반환합니다. 함수는 다음과 같습니다.
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가 있다고 가정하겠습니다.
.
이 방정식은 구간 에 대해 정의되며 다음 경계 조건이 적용됩니다.
,
.
MATLAB®에서 이 방정식을 풀려면 이 방정식을 1계 연립방정식으로 표현하는 함수와 경계 조건에 대한 함수를 작성하고 옵션 값을 설정하고 초기 추측값을 만들어야 합니다. 그러면 BVP 솔버가 4개의 입력값을 사용하여 방정식을 풉니다.
방정식 코딩하기
대입 와 을 사용하여 ODE를 1계 연립방정식으로 바꿀 수 있습니다.
,
.
대응하는 함수는 다음과 같습니다.
function dydx = bvpfcn(x,y) dydx = [y(2) -2*y(2)/x - y(1)/x^4]; end
참고: 모든 함수는 예제 끝에 로컬 함수로 포함되어 있습니다.
경계 조건 코딩하기
경계 조건 함수에서는 경계 조건이 과 같은 형식을 가져야 합니다. 이 형식에서 경계 조건은 다음과 같습니다.
,
.
대응하는 함수는 다음과 같습니다.
function res = bcfcn(ya,yb) res = [ya(1) yb(1)-sin(1)]; end
옵션 설정하기
두 솔버 사이의 오차 제어 차이점을 보기 위해 bvpset
를 사용하여 솔버 통계량 표시를 설정하고 엄격하지 않은 허용오차를 지정합니다. 효율을 높이기 위해 다음과 같은 해석적 야코비 행렬을 지정합니다.
.
야코비 행렬의 값을 반환하는 함수는 다음과 같습니다.
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
를 사용하여 해의 초기 추측값을 만듭니다. 구간 에 있는 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.
결과 플로팅하기
비교를 위해 에 대한 두 가지 계산의 결과를 해석적 해와 함께 플로팅합니다. 해석적 해는 다음과 같습니다.
,
.
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)
를 사용하십시오. 반환 값 dydx
는 f(x,y)에 대응하는 single
이나 double
데이터형의 열 벡터입니다. odefun
은 입력 인수 x
와 y
모두를 받아야만 하는데, 이는 이 중 하나가 함수에 사용되지 않더라도 마찬가지입니다.
예를 들어, 을 풀려면 다음 함수를 사용하십시오.
function dydx = odefun(x,y) dydx = 5*y-3; end
연립방정식에 대한 odefun
의 출력값은 벡터입니다. 벡터의 각 요소는 하나의 방정식의 해입니다. 예를 들어 다음을 풀려면
다음 함수를 사용하십시오.
function dydx = odefun(x,y) dydx = zeros(2,1); dydx(1) = y(1)+2*y(2); dydx(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
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
은 다음 필드를 가집니다.
필드 | 설명 |
---|---|
|
|
|
|
|
|
|
|
|
|
| 해와 관련 있는 계산 비용 통계량, 즉 메시 점의 개수, 잔차 오차, 그리고 |
세부 정보
다중 점 경계값 문제(Multipoint Boundary Value Problem)
다중 점 경계값 문제의 경우, 경계 조건이 적분 구간의 여러 점에 적용됩니다.
bvp4c
는 a = a0 < a1 < a2 < ...< an = b가 구간 [a,b]의 경계점인, 다중 점 경계값 문제를 풀 수 있습니다. 점 a1,a2,...,an−1은 [a,b]를 여러 영역으로 나누는 접점을 나타냅니다. bvp4c
는 1부터 시작하는 인덱스를 사용하여 왼쪽에서 오른쪽으로(a에서 b로) 영역을 나열합니다. 영역 k, [ak−1,ak]에서 bvp4c
는 다음과 같이 도함수를 계산합니다.
yp = odefun(x,y,k)
경계 조건 함수 bcfun(yleft,yright)
에서 yleft(:,k)
는 [ak−1,ak]의 왼쪽 경계에서의 해입니다. 마찬가지로 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
가 있는 다음 형식의 문제를 비롯하여 특이 경계값 문제의 클래스를 풉니다.
구간은 [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차 정확도를 가지며 C1 연속해를 제공하는 선점 다항식입니다. 메시 선택과 오차 제어는 연속해(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
을 사용해 코드 실행 속도를 높일 수 있습니다.
이 함수는 스레드 기반 환경을 완전히 지원합니다. 자세한 내용은 스레드 기반 환경에서 MATLAB 함수 실행하기 항목을 참조하십시오.
버전 내역
R2006a 이전에 개발됨
MATLAB 명령
다음 MATLAB 명령에 해당하는 링크를 클릭했습니다.
명령을 실행하려면 MATLAB 명령 창에 입력하십시오. 웹 브라우저는 MATLAB 명령을 지원하지 않습니다.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list:
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)