연속법(Continuation)을 사용하여 BVP 일관성 확인하기
이 예제에서는 연속법(Continuation)을 사용하여 BVP 해를 더 큰 구간으로 조금씩 확장하는 방법을 보여줍니다.
포크너-스칸(Falkner-Skan) 경계값 문제 [1]은 평판 위를 흐르는 점성 있는 비압축성 층류의 상사해(Similarity Solution)를 구할 때 제기됩니다. 예제 수식은 다음과 같습니다.
.
문제는 로 무한 구간 에 있으며 다음 경계 조건이 적용됩니다.
,
,
.
무한 구간에서는 BVP를 풀 수 없으며 매우 큰 유한 구간에서 BVP를 푸는 것은 실용적이지 않습니다. 대신, 이 예제에서는 작은 구간 에 있는 일련의 여러 문제를 풀어 일 때 해가 일관된 동작을 하는지 확인합니다. 이렇게 문제를 더 간단한 문제로 나누어 각 문제의 해를 다음 문제의 초기 추측값으로 활용하는 방법을 연속법(Continuation)이라고 합니다.
MATLAB®에서 이 연립방정식을 풀려면 경계값 문제 솔버 bvp4c
를 호출하기 전에 방정식, 경계 조건 및 옵션을 코딩해야 합니다. 필요한 함수를 이 예제와 같이 파일 끝에 로컬 함수로 포함시킬 수도 있고, MATLAB 경로에 있는 디렉터리에 이름이 지정된 별도의 파일로 저장할 수도 있습니다.
방정식 코딩하기
방정식을 코딩하는 함수를 만듭니다. 이 함수는 시그니처 dfdx = fsode(x,f)
를 사용해야 합니다. 여기서,
x는 독립 변수입니다.
f는 종속 변수입니다.
대입 , 및 를 사용하여 3계 방정식을 1계 연립방정식으로 다시 작성할 수 있습니다. 다음과 같은 방정식이 됩니다.
,
,
.
대응하는 함수는 다음과 같습니다.
function dfdeta = fsode(x,f) b = 0.5; dfdeta = [ f(2) f(3) -f(1)*f(3) - b*(1 - f(2)^2) ]; end
참고: 모든 함수는 예제 끝에 로컬 함수로 포함되어 있습니다.
경계 조건 코딩하기
이제 경계점에서 경계 조건의 잔차 값을 반환하는 함수를 작성합니다. 이 함수는 시그니처 res = fsbc(f0,finf)
를 사용해야 합니다. 여기서,
f0
은 구간의 시작 부분에서의 경계 조건 값입니다.finf
는 구간의 끝부분에서의 경계 조건 값입니다.
잔차 값을 계산하려면 형식의 경계 조건을 넣어야 합니다. 이 형식에서 경계 조건은 다음과 같습니다.
,
,
.
대응하는 함수는 다음과 같습니다.
function res = fsbc(f0,finf) res = [f0(1) f0(2) finf(2) - 1]; end
초기 추측값 만들기
마지막으로 해에 대한 초기 추측값을 제공해야 합니다. 구간 에 대한 수렴을 구하기 위해서는 경계 조건을 충족하는 상수 추측값 하나와 5개 점으로 구성된 엉성한 메시면 충분합니다. 변수 infinity
는 적분 구간의 오른쪽 한계를 나타냅니다. 후속 반복에서 infinity
값이 3에서 최댓값 6으로 증가함에 따라 각각의 이전 반복의 해가 다음 반복에 대한 초기 추측값으로 적용됩니다.
infinity = 3; maxinfinity = 6; solinit = bvpinit(linspace(0,infinity,5),[0 0 1]);
방정식 풀기 및 해 플로팅하기
초기 구간 에서 문제를 풉니다. 해를 플로팅하고 의 값을 해석적 값 [1]과 비교합니다.
sol = bvp4c(@fsode,@fsbc,solinit); x = sol.x; f = sol.y; plot(x,f(2,:),x(end),f(2,end),'o'); axis([0 maxinfinity 0 1.4]); title('Falkner-Skan Equation, Positive Wall Shear, \beta = 0.5.') xlabel('x') ylabel('df/dx') hold on
fprintf('Cebeci & Keller report that f''''(0) = 0.92768.\n')
Cebeci & Keller report that f''(0) = 0.92768.
fprintf('Value computed using infinity = %g is %7.5f.\n', ... infinity,f(3,1))
Value computed using infinity = 3 is 0.92915.
이제 각 반복에 대한 infinity
의 값을 증가시켜 점점 큰 구간에서 문제를 풉니다. bvpinit
는 다음 차례 infinity
값의 초기 추측값으로 작용하도록 각 해를 새 구간에 외삽합니다. 각 반복은 계산된 값을 출력하고 이전 해 위에 해의 플롯을 겹쳐 놓습니다. infinity = 6
일 때, 해의 일관된 동작이 분명하게 드러나고 의 값은 예측 값에 매우 가깝습니다.
for Bnew = infinity+1:maxinfinity solinit = bvpinit(sol,[0 Bnew]); % Extend solution to new interval sol = bvp4c(@fsode,@fsbc,solinit); x = sol.x; f = sol.y; fprintf('Value computed using infinity = %g is %7.5f.\n', ... Bnew,f(3,1)) plot(x,f(2,:),x(end),f(2,end),'o'); drawnow end
Value computed using infinity = 4 is 0.92774. Value computed using infinity = 5 is 0.92770. Value computed using infinity = 6 is 0.92770.
hold off
로컬 함수(Local Function)
BVP 솔버 bvp4c
가 해를 계산하기 위해 호출하는 로컬 헬퍼 함수는 다음과 같습니다. 또는 이러한 함수를 MATLAB 경로에 있는 디렉터리에 고유의 파일로 저장할 수도 있습니다.
function dfdeta = fsode(x,f) % equation being solved dfdeta = [ f(2) f(3) -f(1)*f(3) - 0.5*(1 - f(2)^2) ]; end %------------------------------------------- function res = fsbc(f0,finf) % boundary conditions res = [f0(1) f0(2) finf(2) - 1]; end %-------------------------------------------
참고 문헌
[1] Cebeci, T. and H. B. Keller. "Shooting and Parallel Shooting Methods for Solving the Falkner-Skan Boundary-layer Equation." J. Comp. Phys., Vol. 7, 1971, pp. 289-300.