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

μ의 크기가 더 큰 경우 문제의 경직성(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')

경계값 문제

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

이 특정 경계값 문제에는 정확히 두 개의 해가 있습니다. 초기 추측값을 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

지연 미분 방정식

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

편미분 방정식

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

참고 항목

| |

관련 항목