MATLAB

最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

微分方程式

この例では、MATLAB を使用して 3 つの異なる微分方程式に対する公式と解法を示します。MATLAB は、広範な微分方程式を解くためのさまざまな数値的アルゴリズムを提供します。

開始値問題

VANDERPOLDEMO は、van der Pol 方程式を定義する関数です。

type vanderpoldemo
function dydt = vanderpoldemo(t,y,Mu)
%VANDERPOLDEMO Defines the van der Pol equation for ODEDEMO.

% Copyright 1984-2002 The MathWorks, Inc.  
% $Revision: 1.1.8.13.2.2 $  $Date: 2014/02/05 21:22:25 $ 

dydt = [y(2); Mu*(1-y(1)^2)*y(2)-y(1)];

方程式は 2 つの 1 次 ODE のシステムとして記述されます。これらは、パラメーター Mu のさまざまな値に対して評価されます。より速く積分を行うには、パラメーター Mu の値を元に適切なソルバーを選択します。

Mu = 1 の場合、MATLAB ODE ソルバーのいずれでも van der Pol 方程式を効率的に解くことができます。以下で使用する ODE45 のソルバーは、そのような例の 1 つです。方程式は、領域 [0, 20] 内で求められます。

tspan = [0, 20];
y0 = [2; 0];
Mu = 1;
ode = @(t,y) vanderpoldemo(t,y,Mu);
[t,y] = ode45(ode, tspan, y0);

% Plot of the solution
plot(t,y(:,1))
xlabel('t')
ylabel('solution y')
title('van der Pol Equation, \mu = 1')

Mu の大きなゲインに対して、問題はスティッフになります。より速く積分を行うには、特別な数値メソッドが必要になります。ODE15S、ODE23S、ODE23T、ODE23TB は、スティッフな問題を効率よく解くことができます。

Mu = 1000 の場合の van der Pol 方程式を解くために、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 は、常微分方程式に対する境界値問題を解きます。

例の関数 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-2002 The MathWorks, Inc.  
%   $Revision: 1.1.8.13.2.2 $  $Date: 2014/02/05 21:22:25 $ 

dydx = [ y(2); -abs(y(1)) ];

TWOBC は、TWOODE に対する境界条件をもちます。

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-2002 The MathWorks, Inc.  
%   $Revision: 1.1.8.13.2.2 $  $Date: 2014/02/05 21:22:25 $ 

res = [ ya(1); yb(1) + 2 ];

BVP4C を使用する前に、メッシュで表すために解の推定値を与えなければなりません。ソルバーは、その後、解の推定値を絞り込むときにメッシュを適応させます。

BVPINIT は、ソルバー BVP4C が必要とするフォームで初期推定を作成します。[0 1 2 3 4] の初期メッシュと定数の推定 y(x) = 1, y'(x) = 0 に対して、以下のように BVPINIT を呼び出します。

solinit = bvpinit([0 1 2 3 4],[1; 0]);

この初期推定を使うことで、BVP4C で問題を解くことができます。

(下記の) 解 sol は、その後、DEVAL を使って点 xint で評価され、プロットされます。

sol = bvp4c(@twoode, @twobc, solinit);

xint = linspace(0, 4, 50);
yint = deval(sol, xint);
plot(xint, yint(1,:),'b');
xlabel('x')
ylabel('solution y')
hold on

この特定の境界値問題は、正確には 2 つの解をもちます。もう 1 つの解は、以下の初期推定で得られます。

   y(x) = -1, y'(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,:),'r');
hold off

偏微分方程式

PDEPE は、1 つの空間変数と時間に対する偏微分方程式を解きます。

例 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-2002 The MathWorks, Inc. 
%   $Revision: 1.1.8.13.2.2 $  $Date: 2014/02/05 21:22:25 $

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-2002 The MathWorks, Inc. 
%   $Revision: 1.1.8.13.2.2 $  $Date: 2014/02/05 21:22:25 $

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-2002 The MathWorks, Inc. 
%   $Revision: 1.1.8.13.2.2 $  $Date: 2014/02/05 21:22:25 $

pl = ul;
ql = 0;
pr = pi * exp(-t);
qr = 1;

PDEPE は、空間の離散化 x と、解のスナップショットで必要な時間ベクトル t を必要とします。20 の接点のメッシュを使用してこの問題を解き、t の 5 つの値における解を要求します。最終的に、解の最初の要素を抽出してプロットします。

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

製品評価版の入手または製品の購入