Main Content

미분 방정식

이 예제에서는 MATLAB®을 사용하여 다른 유형의 여러 미분 방정식을 정식화하고 푸는 방법을 보여줍니다. MATLAB에서는 다양한 미분 방정식을 풀 수 있는 여러 가지 수치 알고리즘을 제공합니다.

  • 초기값 문제

  • 경계값 문제

  • 지연 미분 방정식

  • 편미분 방정식

초기값 문제

vanderpoldemo는 반데르폴 방정식을 정의하는 함수입니다.

d2ydt2-μ(1-y2)dydt+y=0.

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계 상미분방정식으로 구성된 연립방정식으로 작성되었습니다. 다양한 값의 파라미터 μ에 대해 이러한 방정식을 계산합니다. 적분 속도를 높이기 위해, 여기서는 μ의 값을 기반으로 적절한 솔버를 선택해야 합니다.

μ=1인 경우 어느 MATLAB ODE 솔버든지 반데르폴 방정식을 효율적으로 풀 수 있습니다. 한 예로 ode45 솔버를 들 수 있습니다. 이 방정식은 정의역 [0,20]에서 초기 조건 y(0)=2dydt|t=0=0을 사용하여 계산됩니다.

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')

Figure contains an axes object. The axes object with title van der Pol Equation, mu blank = blank 1, xlabel t, ylabel solution y contains an object of type line.

μ의 크기가 더 큰 경우 문제의 경직성(Stiff)은 더 강해집니다. 즉, 이러한 문제는 일반적인 기법으로는 계산하기 힘들다는 의미입니다. 대신 적분 속도를 높이기 위해서는 특별한 수치 계산법이 필요합니다. 함수 ode15s, ode23s, ode23t, ode23tb는 경직성 문제를 효율적으로 풀 수 있습니다.

μ=1000인 경우를 계산하는 반데르폴 방정식의 해는 동일한 초기 조건을 갖는 ode15s를 사용합니다. 해의 주기적 움직임을 확인할 수 있으려면 시간 범위를 [0,3000]까지 크게 늘려야 합니다.

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')

Figure contains an axes object. The axes object with title van der Pol Equation, mu blank = 1000, xlabel t, ylabel solution y contains an object of type line.

경계값 문제

bvp4cbvp5c는 상미분 방정식에 대한 경계값 문제를 풉니다.

예제 함수 twoode는 2개의 1계 ODE로 구성된 연립미분방정식입니다. 미분 방정식은 다음과 같습니다.

d2ydt2+|y|=0.

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는 문제에 대해 경계 조건 y(0)=0y(4)=-2를 갖습니다.

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]로 구성된 메시와 y(x)=1y'(x)=0에 대한 상수 추측값의 경우, 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

Figure contains an axes object. The axes object with xlabel x, ylabel solution y contains an object of type line.

이 특정 경계값 문제에는 정확히 두 개의 해가 있습니다. 초기 추측값을 y(x)=-1y'(x)=0으로 변경하여 또 다른 해를 구할 수 있습니다.

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

Figure contains an axes object. The axes object with xlabel x, ylabel solution y contains 2 objects of type line. These objects represent Solution 1, Solution 2.

지연 미분 방정식

dde23, ddesd, ddensd는 다양한 지연을 갖는 지연 미분 방정식을 풉니다. 예제 ddex1, ddex2, ddex3, ddex4, ddex5는 이 솔버들의 사용 방법에 대한 소형 자습서 역할을 합니다.

ddex1 예제에서는 연립미분방정식을 푸는 방법을 보여줍니다.

y1'(t)=y1(t-1)y2'(t)=y1(t-1)+y2(t-0.2)y3'(t)=y2(t).

익명 함수를 사용하여 이러한 방정식을 표현할 수 있습니다.

ddex1fun = @(t,y,Z) [Z(1,1); Z(1,1)+Z(2,2); y(2)];

이 문제의 내역(t0인 경우)은 다음과 같이 상수입니다.

y1(t)=1y2(t)=1y3(t)=1.

내역을 1로 구성된 벡터로 표현할 수 있습니다.

ddex1hist = ones(3,1);

요소를 2개 가진 벡터는 연립방정식에서 지연을 나타냅니다.

lags = [1 0.2];

함수, 지연, 해 내역, 적분 구간 [0,5]를 솔버에 입력값으로 전달합니다. 솔버는 플로팅에 적합한 전체 적분 구간에 대해 연속해를 생성합니다.

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');

Figure contains an axes object. The axes object with title An example of Wille and Baker DDE with Constant Delays, xlabel time t, ylabel solution y contains 3 objects of type line. These objects represent y_1, y_2, y_3.

편미분 방정식

pdepe는 하나의 공간 변수와 시간에 대해 편미분 방정식을 풉니다. 예제 pdex1, pdex2, pdex3, pdex4, pdex5pdepe의 사용 방법에 대한 소형 자습서 역할을 합니다.

이 예제 문제에서는 함수 pdex1pde, pdex1ic, pdex1bc를 사용합니다.

pdex1pde는 미분 방정식을 정의합니다.

π2ut=x(ux).

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는 초기 조건을 설정합니다.

u(x,0)=sinπx.

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는 경계 조건을 설정합니다.

u(0,t)=0,

πe-t+xu(1,t)=0.

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');

Figure contains an axes object. The axes object with xlabel x, ylabel t contains an object of type surface.

참고 항목

| |

관련 항목