미분 방정식
이 예제에서는 MATLAB®을 사용하여 다른 유형의 여러 미분 방정식을 정식화하고 푸는 방법을 보여줍니다. MATLAB에서는 다양한 미분 방정식을 풀 수 있는 여러 가지 수치 알고리즘을 제공합니다.
초기값 문제
경계값 문제
지연 미분 방정식
편미분 방정식
초기값 문제
vanderpoldemo
는 반데르폴 방정식을 정의하는 함수입니다.
type vanderpoldemo
function dydt = vanderpoldemo(t,y,Mu) %VANDERPOLDEMO Defines the van der Pol equation for ODEDEMO. % Copyright 1984-2014 The MathWorks, Inc. dydt = [y(2); Mu*(1-y(1)^2)*y(2)-y(1)];
이 방정식은 2개의 1계 상미분방정식으로 구성된 연립방정식으로 작성되었습니다. 다양한 값의 파라미터 에 대해 이러한 방정식을 계산합니다. 적분 속도를 높이기 위해, 여기서는 의 값을 기반으로 적절한 솔버를 선택해야 합니다.
인 경우 어느 MATLAB ODE 솔버든지 반데르폴 방정식을 효율적으로 풀 수 있습니다. 한 예로 ode45
솔버를 들 수 있습니다. 이 방정식은 정의역 에서 초기 조건 및 을 사용하여 계산됩니다.
tspan = [0 20]; y0 = [2; 0]; Mu = 1; ode = @(t,y) vanderpoldemo(t,y,Mu); [t,y] = ode45(ode, tspan, y0); % Plot solution plot(t,y(:,1)) xlabel('t') ylabel('solution y') title('van der Pol Equation, \mu = 1')
의 크기가 더 큰 경우 문제의 경직성(Stiff)은 더 강해집니다. 즉, 이러한 문제는 일반적인 기법으로는 계산하기 힘들다는 의미입니다. 대신 적분 속도를 높이기 위해서는 특별한 수치 계산법이 필요합니다. 함수 ode15s
, ode23s
, ode23t
, ode23tb
는 경직성 문제를 효율적으로 풀 수 있습니다.
인 경우를 계산하는 반데르폴 방정식의 해는 동일한 초기 조건을 갖는 ode15s
를 사용합니다. 해의 주기적 움직임을 확인할 수 있으려면 시간 범위를 까지 크게 늘려야 합니다.
tspan = [0, 3000]; y0 = [2; 0]; Mu = 1000; ode = @(t,y) vanderpoldemo(t,y,Mu); [t,y] = ode15s(ode, tspan, y0); plot(t,y(:,1)) title('van der Pol Equation, \mu = 1000') axis([0 3000 -3 3]) xlabel('t') ylabel('solution y')
경계값 문제
bvp4c
와 bvp5c
는 상미분 방정식에 대한 경계값 문제를 풉니다.
예제 함수 twoode
는 2개의 1계 ODE로 구성된 연립미분방정식입니다. 미분 방정식은 다음과 같습니다.
type twoode
function dydx = twoode(x,y) %TWOODE Evaluate the differential equations for TWOBVP. % % See also TWOBC, TWOBVP. % Lawrence F. Shampine and Jacek Kierzenka % Copyright 1984-2014 The MathWorks, Inc. dydx = [ y(2); -abs(y(1)) ];
함수 twobc
는 문제에 대해 경계 조건 및 를 갖습니다.
type twobc
function res = twobc(ya,yb) %TWOBC Evaluate the residual in the boundary conditions for TWOBVP. % % See also TWOODE, TWOBVP. % Lawrence F. Shampine and Jacek Kierzenka % Copyright 1984-2014 The MathWorks, Inc. res = [ ya(1); yb(1) + 2 ];
bvp4c
를 호출하기 전에 메시에 표현할 해에 대한 추측값을 제공해야 합니다. 그러면 솔버가 더 정확한 해를 계산해 나가면서 메시를 적용합니다.
bvpinit
함수는 솔버 bvp4c
에 전달할 수 있는 형식으로 초기 추측값을 조합합니다. [0 1 2 3 4]
로 구성된 메시와 및 에 대한 상수 추측값의 경우, bvpinit
에 대한 호출 구문은 다음과 같습니다.
solinit = bvpinit([0 1 2 3 4],[1; 0]);
초기 추측값이 위와 같은 경우 bvp4c
를 사용하여 문제를 풀 수 있습니다. deval
을 사용하여 일부 점에서 bvp4c
가 반환한 해를 계산한 다음, 결과로 나타나는 값을 플로팅합니다.
sol = bvp4c(@twoode, @twobc, solinit); xint = linspace(0, 4, 50); yint = deval(sol, xint); plot(xint, yint(1,:)); xlabel('x') ylabel('solution y') hold on
이 특정 경계값 문제에는 정확히 두 개의 해가 있습니다. 초기 추측값을 및 으로 변경하여 또 다른 해를 구할 수 있습니다.
solinit = bvpinit([0 1 2 3 4],[-1; 0]); sol = bvp4c(@twoode,@twobc,solinit); xint = linspace(0,4,50); yint = deval(sol,xint); plot(xint,yint(1,:)); legend('Solution 1','Solution 2') hold off
지연 미분 방정식
dde23
, ddesd
, ddensd
는 다양한 지연을 갖는 지연 미분 방정식을 풉니다. 예제 ddex1
, ddex2
, ddex3
, ddex4
, ddex5
는 이 솔버들의 사용 방법에 대한 소형 자습서 역할을 합니다.
ddex1
예제에서는 연립미분방정식을 푸는 방법을 보여줍니다.
익명 함수를 사용하여 이러한 방정식을 표현할 수 있습니다.
ddex1fun = @(t,y,Z) [Z(1,1); Z(1,1)+Z(2,2); y(2)];
이 문제의 내역(인 경우)은 다음과 같이 상수입니다.
내역을 1로 구성된 벡터로 표현할 수 있습니다.
ddex1hist = ones(3,1);
요소를 2개 가진 벡터는 연립방정식에서 지연을 나타냅니다.
lags = [1 0.2];
함수, 지연, 해 내역, 적분 구간 를 솔버에 입력값으로 전달합니다. 솔버는 플로팅에 적합한 전체 적분 구간에 대해 연속해를 생성합니다.
sol = dde23(ddex1fun, lags, ddex1hist, [0 5]); plot(sol.x,sol.y); title({'An example of Wille and Baker', 'DDE with Constant Delays'}); xlabel('time t'); ylabel('solution y'); legend('y_1','y_2','y_3','Location','NorthWest');
편미분 방정식
pdepe
는 하나의 공간 변수와 시간에 대해 편미분 방정식을 풉니다. 예제 pdex1
, pdex2
, pdex3
, pdex4
, pdex5
는 pdepe
의 사용 방법에 대한 소형 자습서 역할을 합니다.
이 예제 문제에서는 함수 pdex1pde
, pdex1ic
, pdex1bc
를 사용합니다.
pdex1pde
는 미분 방정식을 정의합니다.
type pdex1pde
function [c,f,s] = pdex1pde(x,t,u,DuDx) %PDEX1PDE Evaluate the differential equations components for the PDEX1 problem. % % See also PDEPE, PDEX1. % Lawrence F. Shampine and Jacek Kierzenka % Copyright 1984-2014 The MathWorks, Inc. c = pi^2; f = DuDx; s = 0;
pdex1ic
는 초기 조건을 설정합니다.
type pdex1ic
function u0 = pdex1ic(x) %PDEX1IC Evaluate the initial conditions for the problem coded in PDEX1. % % See also PDEPE, PDEX1. % Lawrence F. Shampine and Jacek Kierzenka % Copyright 1984-2014 The MathWorks, Inc. u0 = sin(pi*x);
pdex1bc
는 경계 조건을 설정합니다.
type pdex1bc
function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t) %PDEX1BC Evaluate the boundary conditions for the problem coded in PDEX1. % % See also PDEPE, PDEX1. % Lawrence F. Shampine and Jacek Kierzenka % Copyright 1984-2014 The MathWorks, Inc. pl = ul; ql = 0; pr = pi * exp(-t); qr = 1;
pdepe
에는 해의 스냅샷에 필요한 공간 이산화 x
와 시간 벡터 t
가 필요합니다. 이 예제에서는 20개의 노드로 구성된 메시를 사용하여 이 문제를 푼 후 다섯 개의 t
값에서 해를 구합니다. 그리고 해의 첫 번째 성분을 추출하고 플로팅합니다.
x = linspace(0,1,20); t = [0 0.5 1 1.5 2]; sol = pdepe(0,@pdex1pde,@pdex1ic,@pdex1bc,x,t); u1 = sol(:,:,1); surf(x,t,u1); xlabel('x'); ylabel('t'); zlabel('u');