ドキュメンテーション センター

  • 評価版
  • 製品アップデート

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

常微分方程式

関数の概要

ODE ソルバー

次の表は、初期値問題のソルバー、各ソルバーを使って解くことのできる問題の種類、各ソルバーが使用する手法を示しています。

ソルバー

対象とする問題

手法

ode45

ノンスティッフな微分方程式

ルンゲ・クッタ

ode23

ノンスティッフな微分方程式

ルンゲ・クッタ

ode113

ノンスティッフな微分方程式

Adams

ode15s

スティッフな微分方程式と微分代数方程式

NDF (BDF)

ode23s

スティッフな微分方程式

Rosenbrock

ode23t

中程度のスティッフな微分方程式と微分代数方程式

台形法

ode23tb

スティッフな微分方程式

TR-BDF2

ode15i

完全陰的微分方程式

BDF

評価と拡張

次の関数を使用して、ODE の解の評価と拡張ができます。

関数

説明

deval

ODE ソルバーの出力を使用して数値解を評価

odextend

ODE 初期値問題の解の拡張

ソルバー オプション

options 構造体は名前の付いたプロパティを含み、その値が ODE ソルバーに渡され、問題の解法に影響を与えます。次の関数を使用して options 構造体を作成、変更、アクセスできます。

関数

説明

odeset

ODE ソルバーへ入力する options 構造体の作成と変更

odeget

関数 odeset を使って作成した options 構造体からプロパティを抽出

出力関数

出力関数を設定すると、ソルバーは積分ステップを計算するたびにその関数を呼び出します。関数 odeset を使って、次のサンプル関数を 「OutputFcn」 プロパティに設定することができ、またサンプル関数を変更することでユーザー独自の関数を作成できます。

関数

説明

odeplotodeplot

時系列プロット

odephas2odephas2

2 次元の位相面図

odephas3odephas3

3 次元の位相面図

odeprintodeprint

コマンド ウィンドウに表示

初期値問題

1 次の ODE

常微分方程式 (ODE) には、1 つの独立変数 t に対する従属変数 y の導関数が 1 つ以上含まれます。time は通常、時間を表します。y の t に対する 1 次導関数は y、2 次導関数は y のように表されます。通常、y(t) は要素 y1, y2, ..., yn をもつベクトルです。

MATLAB® ソルバーは、次の形式の 1 次 ODE を取り扱います。

  • 陽的 ODE の形式 y ' = f (t, y)

  • 線形陰的 ODE の形式 M(t, y) y ' = f (t, y), (M(t, y) は行列)

  • 完全陰的 ODE の形式 f (t, y, y ') = 0, (ode15i のみ)

高次の ODE

MATLAB の ODE ソルバーは、1 次の微分方程式のみに対応しています。高次の ODE ソルバーを使用するには、各方程式を、次に示すように等価な 1 次の微分方程式に書き直す必要があります。

y′ = f(t,y)

任意の常微分方程式

y(n) = f(t,y,y′,...,y(n − 1))

は下記の置き換えを行うことにより、1 次の方程式に書き換えることができます。

y1 = y, y2 = y′,..., yn = y(n − 1)

y1= y, y2 = y', ... , yn = y(n − 1)

その結果、等価な方程式として次に示す n 個の 1 次の ODE が得られます。

たとえば、次に示す 2 次のファン デル ポールの方程式を書き換えることができます。

書き換えには、次の代入を行います。

その結果得られる 1 次の ODEは、次のようになります。

初期値

一般的に、与えられた ODE を満たす関数 y(t) は多数あり、求める解を指定するには付加的な情報が必要となります。「初期値問題」では、ある特定の「初期条件」を満たす解を求めます。この初期条件は、y が初期時刻 t0 において y0 に等しいという形式で指定します。ODE の初期値問題は、次のようになります。

関数 f (t, y) が十分に滑らかな場合、この問題は唯一の解をもちます。一般的にこの解の解析式はないため、数値的な方法によって y(t) を近似する必要があります。

ソルバーのタイプ

ノンスティッフな問題

ノンスティッフな問題に対しては、3 つのソルバーが用意されています。

ode45

陽的ルンゲ・クッタ (4,5) 式をベースにした Dormand-Prince 法です。これは、単段階ソルバーです。つまり、y(tn) の計算では、直前の y(tn–1) における解のみ必要とします。一般に、ode45 は大部分の問題に対して最初の試みとして適用するのに最適な関数です。

ode23

陽的ルンゲ・クッタ (2,3) 式をベースにした Bogacki and Shampine による手法です。粗い許容範囲で計算する場合や少しスティッフなシステムの場合、ode45 より有効になることがあります。ode45 と同様に ode23 は、1 ステップ ソルバーです。

ode113

可変次数の Adams-Bashforth-Moulton PECE 法ソルバーです。厳しい許容範囲が設定された問題や ODE の計算に特に時間のかかるような問題では、ode45 より効率的な場合があります。ode113 は複数ステップ ソルバーです。すなわち、現在の点の解を求めるために、それより前の数点を使います。

スティッフな問題

すべての難解な問題が必ずしもスティッフであるとは言えませんが、スティッフな問題はすべて、専用に設計されたソルバーでなければ解くことが難しくなります。スティッフな問題に対するソルバーは、他のソルバーと同じように使うことができます。問題に付加的な情報を与えることで、多くの場合、スティッフなソルバーの効率を高めることができます 「積分器のオプション」を参照)。

スティッフな問題に対するソルバーには、次の 4 種類があります。

ode15s

数値微分法 (NDF) に基づく可変次数のソルバーです。オプションとして、後退差分法 (BDF) (Gear 法とも呼ばれます) も使います。ode113 と同様に、ode15s も多段ソルバーです。問題がスティッフな可能性がある場合や、ode45 がうまく機能しなかったり効率が悪い場合は、ode15s を使ってください。

ode23s

2 次の Rosenbrock の公式を改善したものをベースにしています。1 ステップ ソルバーなので、粗い許容範囲の場合には ode15s よりも効率的な場合があります。ode15s では効果的でないスティッフな問題を解くことができます。

ode23t

「フリー」の内挿を使って台形法を実行します。問題が適度にスティッフで数値的減衰のない解が必要な場合に、このソルバーを使用してください。

ode23tb

最初の段階で台形法ステップを使い、2 番目の段階で 2 次の後退微分方程式を使う陰的ルンゲ・クッタ法、TR-BDF2 を実行します。ode23s と同様に、このソルバーは粗い許容範囲の場合には、ode15s よりも効率的に機能する場合があります。

完全陰的 ODE

ソルバー ode15i は、次の形式の完全陰的微分方程式

f(t,y,y') = 0

を可変次数 BDF 法により解きます。ode15i の基本構文は、次のようになります。

[t,y] = ode15i(odefun,tspan,y0,yp0,options)

入力引数は以下になります。

odefun

この関数は、f(t,y,y') = 0 の形式の微分方程式の左辺を計算します。

tspan

積分区間 [t0,tf] を指定するベクトル。特定の時間 (すべて増加またはすべて減少) での解を得るには、tspan = [t0,t1,...,tf] を使用します。

y0, yp0

それぞれ、y(t0) と y'(t0) に対する初期条件のベクトルです。指定した値は、矛盾がないようにします。すなわち f(t0,y0,yp0) = 0 を満たします。

options

関数 odeset を使用して作成するオプション積分引数です。詳細は、関数 odeset のリファレンス ページを参照してください。

出力引数は、離散点で近似された解を含みます。

t

時刻の列ベクトル。

y

解の配列。y の各行は、t の対応する行に返される時刻での解です。

これらの引数の詳細は、関数 ode15i のリファレンス ページを参照してください。

ソルバー構文

ode15i を除くすべての ODE ソルバー関数は、さまざまな数値法を試しやすくするために、同形式の構文を使用しています。さまざまな解法を同じ問題に適用するには、単に ODE ソルバー関数名を変えます。すべてのソルバー関数で共通する最も簡単な構文は、次のとおりです。

[t,y] = solver(odefun,tspan,y0,options)

solver は、これまでに説明した ODE ソルバー関数の 1 つを示します。

基本的な入力引数は、次のとおりです。

odefun

ODE を計算する関数ハンドル。関数は次に示す形式になります。

dydt = odefun(t,y)

t はスカラー、dydty は列ベクトルです。

tspan

積分区間を指定するベクトル。ソルバーは tspan(1) で初期条件を設定し、tspan(1) から tspan(end) まで積分します。

y0

問題の初期条件ベクトル。

初期値問題」も参照してください。

options

既定の積分プロパティを変更するオプション パラメーターの構造体。

積分器のオプション」では、構造体の作成方法と指定できるプロパティについて説明します。

出力引数は、離散点で近似された解を含みます。

t

時刻の列ベクトル。

y

解の配列。y の各行は、t の対応する行に返される時刻での解です。

これらの引数の詳細は、ODE ソルバーのリファレンス ページを参照してください。

積分器のオプション

ODE ソルバーの既定の積分プロパティは、一般的な問題を扱うために選択されています。これらの既定値を変更することにより、ODE ソルバーのパフォーマンスを向上できる場合があります。これは、1 つまたは複数のプロパティ値を指定する options 構造体をソルバーに渡して行うことができます。

たとえばソルバーの相対許容誤差の値を既定値 1e-3 から 1e-4 に変更するには、以下のようにします。

  1. 関数 odeset を使用して options 構造体を作成します。

    options = odeset('RelTol', 1e-4);
  2. options 構造体をソルバーに渡します。

    • ode15i 以外のすべてのソルバーに次の構文を使用します。

      [t,y] = solver(odefun,tspan,y0,options)
      
    • ode15i には次の構文を使用します。

      [t,y] = ode15i(odefun,tspan,y0,yp0,options)

使用可能なオプションの詳細は、関数 odeset のリファレンス ページを参照してください。

van der Pol 方程式 (ノンスティッフ)

この例は、ODE 問題の初期値問題を解くステップを示しています。

  1. 問題を 1 次の ODE 問題に書き換えます。 以下の van der Pol 方程式 (2次) を書き換えます。

    ここで μ > 0 はスカラー パラメーターで、置換により y'1 = y2 です。その結果得られる 1 次の ODEは、次のようになります。

  2. 1 次の ODE をコード化する。1 次の ODE として方程式を表現してから、ODE ソルバーが使用できる関数にコード化します。関数は次の形式になります。

    dydt = odefun(t,y)

    ty が、関数の 2 つの引数になる必要がありますが、関数がそれらを使用しなくてもかまいません。出力 dydt は列ベクトルで、y の導関数です。

    次のコードは関数 vdp1 を使って表した van der Pol 方程式です。関数 vdp1 は μ = 1 と仮定します。変数 y1 と y2 は、2 要素のベクトルのエントリ y(1)y(2) です。

    function dydt = vdp1(t,y)
    dydt = [y(2); (1-y(1)^2)*y(2)-y(1)];

    関数 vdp1 は引数 ty を受け入れる必要がありますが、計算では t を使用しません。

  3. 問題にソルバーを適用します。

    問題にどのソルバーを使用するのかを決めます。次にソルバーを呼び出して、1 次の ODE を記述した関数と、問題を解く時間区間、および初期条件ベクトルを渡します。

    van der Pol 方程式に対して、初期値 y(1) = 2y(2) = 0 を使って、時間区間 [0 20]ode45 を使用します。

    [t,y] = ode45(@vdp1,[0 20],[2; 0]);
    

    この例では関数 vdp1ode45 への関数ハンドルとして渡すために、@ を使います。その結果、出力として時間の列ベクトル t と解の配列 y が得られます。y の各行は、t の対応する行に返される時刻と対応します。y の 1 列目は、y1 に対応し、2 列目は y2 に対応します。

  4. ソルバーの出力を確認します。plot コマンドを使うと、ソルバーの出力をプロットできます。

    plot(t,y(:,1),'-',t,y(:,2),'--')
    title('Solution of van der Pol Equation, \mu = 1');
    xlabel('time t');
    ylabel('solution y');
    legend('y_1','y_2')

ソルバーの出力関数を使って、出力を処理することもできます。ソルバーはタイム ステップの計算するたびに、積分プロパティ OutputFcn に設定した関数を呼び出します。求める関数に OutputFcn を設定するには、関数 odeset を使用します。OutputFcn の詳細は、関数 odeset のリファレンス ページの「ソルバー出力プロパティ」を参照してください。

van der Pol 方程式 (スティッフ)

この例は、スティッフな問題を取り扱います。スティッフな問題では、解が積分区間に比べて非常に短い時間スケールで変化しますが、求める解 (適切な精度を保った数値解) は、はるかに大きい時間スケールでしか変化しません。スティッフな問題用に設計されていない手法は、可能な限り高速な変化を解くため十分に小さなタイム ステップで計算を行うため、解がゆっくりと変化する区間に対しては効果的ではありません。

μ を大きくして 1000 にすると、van der Pol 方程式の解は急激に変化し、より長い時間スケールの振動を示します。初期値問題の解の近似は、さらに難しい問題になります。この特殊な問題はスティッフなので、ode45 のようなスティッフな問題に対応していないソルバーは非常に効率が悪く、実用に適しません。ソルバー ode15s は、このようなスティッフな問題を解くために設計されています。

関数 vdp1000 は、前の例で示した van del Pol 方程式を計算しますが、μ = 1000 を使用します。

function dydt = vdp1000(t,y)
dydt = [y(2); 1000*(1-y(1)^2)*y(2)-y(1)];

    メモ:   この例では、μ を ODE 関数でハードコード化しています。vdpodevdpode の例でも同じ問題を解いていますが、ユーザー定義の μ を ODE 関数へのパラメーターとして渡しています。

ode15s を使用して問題を解きます。初期条件ベクトルには [2; 0] をそのまま使い、計算時間は [0 3000] にします。スケーリング (積分ステップ間隔とシステムの振幅周期を調べる) のために、y(t) の最初の成分をプロットします。

[t,y] = ode15s(@vdp1000,[0 3000],[2; 0]);
plot(t,y(:,1),'-');
title('Solution of van der Pol Equation, \mu = 1000');
xlabel('time t');
ylabel('solution y_1');

van der Pol 方程式 (ODE のパラメーター化)

前説ではパラメーター µ を 2 つの値に変えて van der Pol 方程式を解く方法を示しました。この2つの例では、µ の値が 1 もしくは 1000 とハードコードされています。さまざまなパラメーター値について ODE を解く場合、ODE 関数にそのパラメーターを指定し、ODE ソルバーを実行するたびにそのパラメーターに値を代入すると、より便利になります。この節では、van der Pol 方程式に対してこの処理を行う方法を説明します。

ODE 関数にパラメーター値を与える方法の 1 つは、次のような MATLAB ファイルを記述することです。

  • 入力としてパラメーターを受け入れる

  • ODE を入れ子関数として含み、内部に入力パラメーターを含める

  • ODE ソルバーを呼び出す

以下のコードは、この処理を説明しています。

function [t,y] = solve_vdp(mu)
tspan = [0 max(20, 3*mu)];
y0 = [2; 0];

% Call the ODE solver ode15s.
[t,y] = ode15s(@vdp,tspan,y0);

    % Define the ODE function as nested function, 
    % using the parameter mu.
    function dydt = vdp(t,y)
    dydt = [y(2); mu*(1-y(1)^2)*y(2)-y(1)];
    end
end

ODE を表す関数 vdp は入れ子関数であるため、パラメーター mu の値を使用することができます。

mu = 1 として MATLAB ファイルを実行するには、次のように入力します。

[t,y] = solve_vdp(1);

µ = 1000 としてコードを実行するには、次のように入力します。

[t,y] = solve_vdp(1000);

これらの関数に基づく例は、関数 vdpodevdpode のコードを参照してください。

van der Pol 方程式 (解の評価)

ODE ソルバーで採用されている数値解法は、積分区間 [a,b] に対し連続解を作成します。関数 deval とソルバーで返される構造体 sol を使って、近似解 S(x)を [a,b] 内の任意の点で評価できます。たとえば van der Pol 方程式 (ノンスティッフ) の問題を、次のように 1 つの出力引数 sol

sol = ode45(@vdp1,[0 20],[2; 0]);

ode45 を呼び出して解く場合、ode45 は解を構造体として返します。ベクトル xint = 1:5 の各点で、次のように近似解を評価できます。

xint = 1:5;
Sxint = deval(sol,xint)

Sxint =

    1.5081    0.3235   -1.8686   -1.7407   -0.8344
   -0.7803   -1.8320   -1.0220    0.6260    1.3095

関数 deval はベクトル化されています。ベクトル xint に対して、Sxinti 番目の列は解 y(xint(i)) を近似します。

オイラー方程式 (ノンスティッフ)

rigidode は、ノンスティッフな問題 [8] を解くために Krogh が示した標準テスト問題の解法を示します。

この ODE 方程式は、外力のない剛体のオイラー方程式から構成されます。

便宜上、1 つの MATLAB ファイルで問題全体を定義して解きます。微分方程式は、ローカル関数 f としてコード化されています。この例では ode45 ソルバーを出力引数なしで呼び出すので、既定の出力関数 odeplotodeplot が使われて解成分をプロットします。

この例を実行するには、コマンド ラインで「rigidoderigidode」と入力してください。

function rigidode 
%RIGIDODE  Euler equations: rigid body without external forces
tspan = [0 12];
y0 = [0; 1; 1];

% Solve the problem using ode45
ode45(@f,tspan,y0);
% ------------------------------------------------------------
function dydt = f(t,y)
dydt = [ y(2)*y(3) 
        -y(1)*y(3) 
        -0.51*y(1)*y(2) ];

完全陰的 ODE

次の例では、Weissinger の式により定義される陰的 ODE 問題を解くための関数 ode15i の使用方法を示します。

ty2(y')3 – y3(y')2 + t(t2 + 1)y' – t2y = 0

初期値は y(1) = (3/2)1/2 とします。ODE の厳密解は以下になります。

y(t) = (t2 + 0.5)1/2

例では、関数 weissinger を使用します。この関数は MATLAB で提供されており、方程式の左辺を計算します。

例では ode15i を呼び出す前に、y'(t0) に対して矛盾のない初期値を計算するために補助関数 decic を使用します。次の呼び出しでは、与えられた初期値 y(1) = (3/2)1/2 は固定され、y'(1) に対して推定値 0 が指定されます。詳細は、decic のリファレンス ページを参照してください。

t0 = 1;
y0 = sqrt(3/2);
yp0 = 0;
[y0,yp0] = decic(@weissinger,t0,y0,1,yp0,0);

ode15i を呼び出して ODE を解くことができ、次のコマンドで解析解に対する数値解をプロットすることができます。

[t,y] = ode15i(@weissinger,[1 10],y0,yp0);
ytrue = sqrt(t.^2 + 0.5);
plot(t,y,t,ytrue,'o');

有限要素の離散化

fem1ode は、偏微分方程式の有限要素の離散化の結果である ODE 方程式の解を示します。fem1ode(N)N の値は、離散化を制御するもので求まる結果の方程式系は、N 個の方程式から構成されます。既定は、N = 19 です。

この例は、質量行列を含んでいます。この ODE 方程式は、次の偏微分方程式

と初期条件

u(0,x) = sin(x)

および以下の境界条件を満たします。

u(t,0) = u(t,π) = 0.

整数 N が選択され、h が π/(N + 1) として定義されます。そして偏微分方程式の解が、k = 0, 1, ..., N+1 のときの xk = kh で次の式によって近似されます。

ここで、ϕk(x) は、xk で 1、その他すべての xj で 0 となる区分線形関数です。Galerkin の離散化は、次の ODE 方程式を導きます。

また、三重対角行列 M(t) と J は次の式で与えられます。

初期値 c(0) は、偏微分方程式の初期条件を使用します。問題は、[0,π] の時間間隔で解かれます。

fem1ode の例の中で、プロパティ

options = odeset('Mass',@mass,'MStateDep','none','Jacobian',J)

M(t)y' = Jy のように表現されることが問題であると示されています。入れ子関数 mass(t) は、時間依存の質量行列 M(t) を計算し、J は定数 Jacobian です。

この例を実行するには、コマンド ラインで「fem1odefem1ode」と入力してください。コマンド ラインで N の値を、fem1ode の引数として設定することができます。既定値は N = 19 です。

function fem1ode(N)
%FEM1ODE  Stiff problem with a time-dependent mass matrix 

if nargin < 1
  N = 19;
end
h = pi/(N+1);
y0 = sin(h*(1:N)');
tspan = [0; pi];

% The Jacobian is constant.
e = repmat(1/h,N,1);    %  e=[(1/h) ... (1/h)];
d = repmat(-2/h,N,1);   %  d=[(-2/h) ... (-2/h)]; 
% J is shared with the derivative function.
J = spdiags([e d e], -1:1, N, N);

d = repmat(h/6,N,1);  
% M is shared with the mass matrix function.
M = spdiags([d 4*d d], -1:1, N, N);

options = odeset('Mass',@mass,'MStateDep','none', ...
                 'Jacobian',J);

[t,y] = ode15s(@f,tspan,y0,options);

figure;
surf((1:N)/(N+1),t,y);
set(gca,'ZLim',[0 1]);
view(142.5,30);
title(['Finite element problem with time-dependent mass ' ...
       'matrix, solved by ODE15S']);
xlabel('space ( x/\pi )');
ylabel('time');
zlabel('solution');
%--------------------------------------------------------------
function yp = f(t,y)
% Derivative function.
   yp = J*y;    % Constant Jacobian provided by outer function
end             % End nested function f
%--------------------------------------------------------------
function Mt = mass(t)
% Mass matrix function.
   Mt = exp(-t)*M;    % M is provided by outer function
end                   % End nested function mass
%--------------------------------------------------------------
end

大規模でスティッフなスパース問題

brussode は、(潜在的に) 大規模でスティッフなスパース問題の解を示します。この問題は、古典的な "Brusselator" 方程式 [3] で、化学反応の拡散モデルです。

解の計算は、時間間隔 [0,10] で α = 1/50

ui(0) = 1 + sin(2πxi)
vi(0) = 3

ここで、i = 1,...,N は xi = i/(N + 1) となります。また、2N 個の方程式がありますが、方程式が u1、v1、u2、v2、... のように順序付けられている場合、ヤコビは、定数幅 5 をもつ帯域になります。

呼び出し brussode(N) では、N は N に相当し、パラメーター N2 はグリッド ポイントの数を指定します。結果のシステムは、2N 個の方程式から構成されます。既定では N20 です。この問題は、N が大きくなるにつれてスティッフ性が強くなると共に、ヤコビのスパース性が増大します。

入れ子関数 f(t,y) は、Brusselator 問題用の導関数ベクトルを返します。ローカル関数 jpattern(N) は、次に示すヤコビ行列 ∂f/∂y 内の非ゼロの位置を示す 10 のスパース行列を返します。この例は、この行列を JPattern プロパティに代入し、ソルバーはこのスパース パターンを使用してヤコビをスパース行列として数値的に生成します。スパースパターンを与えることは、ヤコビを作成するために必要な関数値の算出回数を非常に減らし、積分を速くすることになります。

Brusselator 問題にスパース パターンを与えなければ、2N 行 2N 列のヤコビ行列を計算するために、2N 回の関数計算が必要となります。スパース パターンが与えられた場合、N の値にかかわらず、4 回の計算のみが必要となります。

この例を実行するには、コマンド ラインで「brussodebrussode」と入力してください。コマンド ラインで N の値を、brussode の引数として設定することができます。既定値は N = 20 です。

function brussode(N)
%BRUSSODE  Stiff problem modeling a chemical reaction 

if nargin < 1
  N = 20;
end

tspan = [0; 10];
y0 = [1+sin((2*pi/(N+1))*(1:N));
repmat(3,1,N)];

options = odeset('Vectorized','on','JPattern',jpattern(N));

[t,y] = ode15s(@f,tspan,y0,options);

u = y(:,1:2:end);
x = (1:N)/(N+1);
surf(x,t,u);
view(-40,30);
xlabel('space');
ylabel('time');
zlabel('solution u');
title(['The Brusselator for N = ' num2str(N)]);
% --------------------------------------------------------------
function dydt = f(t,y)
c = 0.02 * (N+1)^2;
dydt = zeros(2*N,size(y,2));      % preallocate dy/dt
% Evaluate the two components of the function at one edge of 
% the grid (with edge conditions).
i = 1;
dydt(i,:) = 1 + y(i+1,:).*y(i,:).^2 - 4*y(i,:) + ...
            c*(1-2*y(i,:)+y(i+2,:));
dydt(i+1,:) = 3*y(i,:) - y(i+1,:).*y(i,:).^2 + ...
              c*(3-2*y(i+1,:)+y(i+3,:));
% Evaluate the two components of the function at all interior 
% grid points.
i = 3:2:2*N-3;
dydt(i,:) = 1 + y(i+1,:).*y(i,:).^2 - 4*y(i,:) + ...
            c*(y(i-2,:)-2*y(i,:)+y(i+2,:));
dydt(i+1,:) = 3*y(i,:) - y(i+1,:).*y(i,:).^2 + ...
              c*(y(i-1,:)-2*y(i+1,:)+y(i+3,:));
% Evaluate the two components of the function at the other edge 
% of the grid (with edge conditions).
i = 2*N-1;
dydt(i,:) = 1 + y(i+1,:).*y(i,:).^2 - 4*y(i,:) + ...
            c*(y(i-2,:)-2*y(i,:)+1);
dydt(i+1,:) = 3*y(i,:) - y(i+1,:).*y(i,:).^2 + ...
              c*(y(i-1,:)-2*y(i+1,:)+3);
end   % End nested function f
end   % End function brussode
% --------------------------------------------------------------
function S = jpattern(N)
B = ones(2*N,5);
B(2:2:2*N,2) = zeros(N,1);
B(1:2:2*N-1,4) = zeros(N,1);
S = spdiags(B,-2:2,2*N,2*N);
end;

イベント位置

ballode は、バウンドするボールの運動をモデル化します。この例は、ODE ソルバーがイベントの位置を検出する機能を示します。

バウンドするボールの方程式は、次のように表せます。

y'1 = y2
y'2 = –9.8

この例で、イベント関数はローカル関数 events でコード化されています。

[value,isterminal,direction] = events(t,y)

これは、次の値を返します。

  • イベント関数の値

  • value = 0 になった場合にシミュレーションを停止すべきかどうかの情報 (isterminal = 1 または 0)

  • 求めるゼロ交差の方向

-1

負の方向へのゼロ交差のみを検出

0

すべてのゼロ交差の検出

1

正の方向へのゼロ交差のみを検出

valueisterminaldirection の長さは、event 関数の数と同じになります。各ベクトルの i 番目の要素は、i 番目のイベント関数に対応します。高度なイベントの発生位置を検出する例は、orbitode (「高度なイベント位置検出問題」) を参照してください。

ballodeEvents プロパティを @events に設定すると、ボールが落ちる間 (direction = -1)、ボールが地面にヒットした (高さ y(1)0 になった) とき、シミュレーションを停止 (isterminal = 1) することになります。この例では跳ね返ったボールに、対応する初期条件を与えて積分を再開します。

この例を実行するには、コマンド ラインで「ballodeballode」と入力してください。

function ballode

tstart = 0;
tfinal = 30;
y0 = [0; 20];
refine = 4;
options = odeset('Events',@events,'OutputFcn', @odeplot,...
                 'OutputSel',1,'Refine',refine);

set(gca,'xlim',[0 30],'ylim',[0 25]);
box on
hold on;

tout = tstart;
yout = y0.';
teout = [];
yeout = [];
ieout = [];
for i = 1:10
  % Solve until the first terminal event.
  [t,y,te,ye,ie] = ode23(@f,[tstart tfinal],y0,options);
  if ~ishold
    hold on
  end
  % Accumulate output.  
  nt = length(t);
  tout = [tout; t(2:nt)];
  yout = [yout; y(2:nt,:)];
  teout = [teout; te];    % Events at tstart are never reported.
  yeout = [yeout; ye];
  ieout = [ieout; ie];

  ud = get(gcf,'UserData');
  if ud.stop
    break;
  end
  
  % Set the new initial conditions, with .9 attenuation.
  y0(1) = 0;
  y0(2) = -.9*y(nt,2);

  % A good guess of a valid first time step is the length of 
  % the last valid time step, so use it for faster computation.
  options = odeset(options,'InitialStep',t(nt)-t(nt-refine),...
                           'MaxStep',t(nt)-t(1));
  tstart = t(nt);
end

plot(teout,yeout(:,1),'ro')
xlabel('time');
ylabel('height');
title('Ball trajectory and the events');
hold off
odeplot([],[],'done');
% --------------------------------------------------------------
function dydt = f(t,y)
dydt = [y(2); -9.8];
% --------------------------------------------------------------
function [value,isterminal,direction] = events(t,y)
% Locate the time when height passes through zero in a 
% decreasing direction and stop integration.
value = y(1);     % Detect height = 0
isterminal = 1;   % Stop the integration
direction = -1;   % Negative direction only

高度なイベント位置検出問題

orbitode は、ノンスティッフな問題に使用するソルバーの標準テスト問題の解法を示します。これは月の周りを回り、地球に戻ってくる宇宙船の軌跡を示すものです (Shampine and Gordon [8], p.246)。

orbitode 問題は、次に示す 4 つの方程式から構成されるシステムです。

ここで、

最初の 2 つの解要素は無限小の質量をもつ質点の座標系で、一方の解要素を、他方の解要素に対してプロットしたものが、その質点の軌跡になります。初期条件は、軌跡が周期性をもつように選択されています。μ の値は、宇宙船が月の周りを回り、地球に戻るように設定します。軌道の定量的な動作の再生成のためには、適度に厳密な許容範囲が必要です。適切な値は、RelTol1e-5 で、AbsTol1e-4 です。

入れ子関数 events は、開始点から最も遠い点の位置を示すイベント関数と、宇宙船が開始点に戻る時刻を示すイベント関数を含みます。積分に使われるステップ サイズはイベント位置には依存しませんが、イベント位置自体は正確に検出されることに注意してください。この例では、ゼロクロッシングの方向を決める機能が重要な役割を果たします。初期点の折り返し点と最大距離の点は同じイベント関数値をもちますが、交差の方向を使って、それらを区別します。

この例を実行するには、コマンド ラインで「orbitodeorbitode」と入力してください。この例は、出力関数「odephas2odephas2」を使って 2 次元の位相平面プロットを作成するので、積分の進行状況を見ることができます。

function orbitode
%ORBITODE  Restricted three-body problem

mu = 1 / 82.45;
mustar = 1 - mu;
y0 = [1.2; 0; 0; -1.04935750983031990726];
tspan = [0 7];

options = odeset('RelTol',1e-5,'AbsTol',1e-4,...
                 'OutputFcn',@odephas2,'Events',@events);

[t,y,te,ye,ie] = ode45(@f,tspan,y0,options);

plot(y(:,1),y(:,2),ye(:,1),ye(:,2),'o');
title ('Restricted three body problem')
ylabel ('y(t)')
xlabel ('x(t)')
% --------------------------------------------------------------
function dydt = f(t,y)
r13 = ((y(1) + mu)^2 + y(2)^2) ^ 1.5;
r23 = ((y(1) - mustar)^2 + y(2)^2) ^ 1.5;
dydt = [ y(3)
         y(4)
         2*y(4) + y(1) - mustar*((y(1)+mu)/r13) - ...
                         mu*((y(1)-mustar)/r23)
        -2*y(3) + y(2) - mustar*(y(2)/r13) - mu*(y(2)/r23) ];
end   % End nested function f
% --------------------------------------------------------------
function [value,isterminal,direction] = events(t,y)
% Locate the time when the object returns closest to the
% initial point y0 and starts to move away; stop integration.
% Also locate the time when the object is farthest from the 
% initial point y0 and starts to move closer.
% 
% The current distance of the body is
% 
%   DSQ = (y(1)-y0(1))^2 + (y(2)-y0(2))^2 
%       = <y(1:2)-y0(1:2),y(1:2)-y0(1:2)>
% 
% A local minimum of DSQ occurs when d/dt DSQ crosses zero 
% heading in the positive direction. Compute d(DSQ)/dt as
% 
%  d(DSQ)/dt = 2*(y(1:2)-y0(1:2))'*dy(1:2)/dt = ... 
%                 2*(y(1:2)-y0(1:2))'*y(3:4)
% 
dDSQdt = 2 * ((y(1:2)-y0(1:2))' * y(3:4));
value = [dDSQdt; dDSQdt];
isterminal = [1; 0];            % Stop at local minimum
direction = [1; -1];            % [local minimum, local maximum]
end   % End nested function events
end 

微分代数方程式

hb1dae は、hb1odehb1ode の例を微分代数方程式 (DAE) 問題として作り替えたものです。hb1ode の例でコード化されている Robertson 問題は、スティッフな ODE 問題を解くコードに対する古典的なテスト問題です。

    メモ:   Robertson 問題は、LSODI   [4] の PROLOG (論理型プログラミング言語の一種) の例です。

hb1ode の中で、この問題は初期条件 y1(0) = 1、y2(0) = 0、y3(0) = 0 を与えることにより、定常状態の解が計算されます。これらの微分方程式は、線形保存則を満たし、DAE として問題を再定式化するのに使います。

これらの式は成分の総和が 1 とならない y(0) の解をもちません。この問題は、My' = f(t,y) の型をしており、

M は特異ですが、hb1dae はこの情報をソルバーに伝えません。ソルバーは問題が ODE ではなく、DAE であることを識別しなければなりません。同様に、初期条件に矛盾がないことは明らかですが、例では矛盾する値 y3(0) = 10–3 を使用して、矛盾のない初期条件の計算を示しています。

この例を実行するには、コマンド ラインで「hb1daehb1dae」と入力してください。hb1dae は、次の特徴をもっています。

  • 他の成分と比べて小さな絶対許容誤差を y2 に課します。これは、y2 が他の要素より非常に小さく、主要な変化が相対的に短い時間で生じるためです。

  • y2 の動作をより明確に示すために、解を計算する点を追加指定します。

  • y2 に 104 を掛けることによって、y2 が他の解と一緒にプロットしたときでも見えるようにします。

  • 対数スケールを使用して、長時間にわたる解のプロットを作成します。

    function hb1dae
    %HB1DAE  Stiff differential-algebraic equation (DAE)
    
    % A constant, singular mass matrix
    M = [1 0 0
         0 1 0 
         0 0 0];
    
    % Use inconsistent initial condition to test initialization.
    y0 = [1; 0; 1e-3];
    tspan = [0 4*logspace(-6,6)];
    
    % Use the LSODI example tolerances.'MassSingular' is left
    % at its default 'maybe' to test the automatic detection
    % of a DAE.
    options = odeset('Mass',M,'RelTol',1e-4,...
                     'AbsTol',[1e-6 1e-10 1e-6],...
                     'Vectorized','on');
    
    [t,y] = ode15s(@f,tspan,y0,options);
    
    y(:,2) = 1e4*y(:,2);
    
    semilogx(t,y);
    ylabel('1e4 * y(:,2)');
    title(['Robertson DAE problem with a Conservation Law, '...
           'solved by ODE15S']);
    xlabel('This is equivalent to the stiff ODEs in HB1ODE.');
    % --------------------------------------------------------
    function out = f(t,y)
    out = [ -0.04*y(1,:) + 1e4*y(2,:).*y(3,:)
             0.04*y(1,:) - 1e4*y(2,:).*y(3,:) - 3e7*y(2,:).^2
             y(1,:) + y(2,:) + y(3,:) - 1 ];

非負の解

解のある要素が非負でなければならない場合、関数 odeset を使用して、これらの要素のインデックスに NonNegative プロパティを設定します。

    メモ:   このオプションは、ode23sode15i、および、質量行列がある問題に適用される陰的ソルバー (ode15sode23tode23tb) には適用されません。

非負性が与える影響は小さくはありません。このオプションは必要な場合のみ、たとえば解の適用や積分ができなくなる場合に使用することをお勧めします。

区間 [0, 40] 上で解かれる次の初期値問題を考えます。

y' = - |y|, y(0) = 1

この問題の解はゼロまで減少します。ソルバーが負の近似解を生成する場合、この値を通り ODE の解の追跡を始めます。解は負の無限大になり、計算できなくなります。NonNegative プロパティを使用すると、これが起こらなくなります。

この例では、ode45 の最初の呼び出しは、ソルバー パラメーターに既定を使用します。

ode = @(t,y) -abs(y);
[t0,y0] = ode45(ode,[0, 40], 1);

2 番目の呼び出しでは、options を使用して非負の制約を課します。

options = odeset('NonNegative',1);
[t1,y1] = ode45(ode,[0, 40], 1, options);

次のプロットは、数値解を厳密解と比較します。

以下は、このプロットを得るために使用されたコードの詳細を表示しています。

ode = @(t,y) -abs(y);
options = odeset('Refine',1);
[t0,y0] = ode45(ode,[0, 40], 1,options);
options = odeset(options,'NonNegative',1);
[t1,y1] = ode45(ode,[0, 40], 1, options);
t = linspace(0,40,1000);
y = exp(-t);
plot(t,y,'b-',t0,y0,'ro',t1,y1,'b*');
legend('Exact solution','No constraints','Nonnegativity', ...
       'Location','SouthWest')

kneeode の例-  kneeode の例では、数値解に非負の制約を課すことによって "knee problem" を解きます。次の初期値問題を考えます。

ε*y' = (1-x)*y - y^2,    y(0) = 1 

0 < ε < 1 に対し、この問題の解は、x < 1 および x > 1 それぞれに対して、ゼロ アイソクライン y = 1 - xy = 0 に近づきます。既定の許容誤差を使用して計算される場合の数値解は、積分の全区間に対して、等斜褶曲 y = 1 - x になります。非負という制約を課すと正確な解になります。

以下は、kneeode の例のコードです。

function kneeode
%KNEEODE  The "knee problem" with Nonnegativity constraints.

% Problem parameter
epsilon = 1e-6;

y0 = 1;
xspan = [0, 2];
 
% Solve without imposing constraints 
options = [];
[x1,y1] = ode15s(@odefcn,xspan,y0,options);
 
% Impose nonnegativity constraint
options = odeset('NonNegative',1);
[x2,y2] = ode15s(@odefcn,xspan,y0,options);
 
figure
plot(x1,y1,'b.-',x2,y2,'g-')
axis([0,2,-1,1]);
title('The "knee problem"');
legend('No constraints','nonnegativity')
xlabel('x');
ylabel('solution y')

   function yp = odefcn(x,y)
      yp = ((1 - x)*y - y^2)/epsilon;
   end  
end  % kneeode

導関数は、入れ子関数 odefcn 内で定義されています。odefcn 内で使用される epsilon の値は、外部関数から得られます。

function yp = odefcn(x,y)
yp = ((1 - x)*y - y^2)/epsilon;
end

この例では、関数 ode15s で最初に既定のオプションを使用し、次に非負であるという制約を課すことによって問題を解きます。例を実行するには、MATLAB コマンド プロンプトで「kneeode」と入力してください。

出力のプロットは次のようになります。このプロットから制約を課した後の解の動作が正しいことが確かめられます。

その他の例

以下の例を使用できます。次のように入力します。

edit examplename

を入力してコードを表示し、

examplename

で例を実行します。

例の名前

説明

amp1dae

スティッフな DAE -電気回路

ballode

簡単なイベントの位置検出問題-バウンドするボール

batonode

時間と状態に依存した質量行列をもつ ODE -バトンの運動

brussode

スティッフな大規模問題-化学反応の拡散作用 (Brusselator)

burgersode

状態量に強く依存した質量行列をもつ ODE -移動メッシュ法を使って Burger の方程式を解く

fem1ode

時変の質量行列をもつスティッフな問題-有限要素法

fem2ode

定質量行列をもつスティッフな問題-有限要素法

hb1ode

非常に長い区間で解くスティッフな ODE 問題- Robertson 化学反応

hb1dae

Robertson の問題-保存則下でのスティッフな線形陰的 DAE 問題

ihb1dae

Robertson の問題-スティッフ、完全陰的 DAE 問題

iburgersode

Burger 方程式を陰的 ODE として解く問題

kneeode

非負の制約をもつ "knee problem"

orbitode

高度なイベント位置検出問題-制約付きの 3 体問題

rigidode

ノンスティッフな問題-外力のない剛体のオイラー方程式

vdpode

パラメーター化した van der Pol 方程式 (μ が大きいときスティッフ)

トラブルシューティング

一般的な質問

質問

回答

ODE ソルバーは integral とどのように異なりますか。

integral は y' = f(t) の形式の問題を解きます。各種の ODE ソルバーは、さらに一般的な問題 y' = f(t,y)、質量行列を含む線形陰的な問題 M(t,y)y' = f(t,y)、および完全に陰的な問題 f(t,y,y') = 0 を解きます。

未知数の数より方程式が多い場合やその逆の場合 (未知数の数より方程式が少ない時)、ODE 方程式を解くことができますか。

いいえ。

メモリと計算効果

質問

回答

ODE 関連関数で、どの程度のサイズの問題を解くことができますか。

基本的な制限要因は、メモリ サイズとシミュレーション時間です。ステップごとにシステムの方程式の数を n で表すとすると、ノンスティッフな問題に対するソルバーは、長さ n のベクトルを代入します。スティッフな問題に対するソルバーは、n のベクトルを割り当てますが、nn 列のヤコビ行列も割り当てます。これらのソルバーでは、スパース オプションを使用すると良い場合があります。

問題がノンスティッフな場合やスパースオプションを使った場合は、数千の未知数を含む問題が解けることもあります。しかし、この場合も結果を維持する保存スペースが問題になります。ソルバーに特定の点のみで解を計算させたり、または出力引数を設定しないでソルバーを呼び出したり、解の観測に出力関数を使用することを試みてください。

非常に大きなシステムを解く予定ですが、y のほんの 2、3 の要素しか求める必要がありません。不要な要素を保存しない方法はありますか。

はい。「ユーザーが設定できる出力関数」の機能が搭載されています。ソルバーの読み込みに出力引数を設定しなかった場合、ソルバーは、すべての要素の履歴保存用に領域を割り当てません。その代わりに、ソルバーは OutputFcn(t,y,flag) をタイム ステップごとに呼び出します。特定の要素の履歴を保持するには、対象となる要素だけを保存するか、プロットする出力関数を記述します。

積分計算の立ち上がり部分でどのくらいの CPU 時間がかかりますか。またその時間を減らす方法とは?

立ち上がりで最も時間がかかるのは、ソルバーが問題のスケールに合ったステップ サイズを見つけようとする部分です。適切なステップ サイズがわかっている場合は、InitialStepプロパティを設定してください。たとえば、イベントの位置検出ループ内で積分器を繰り返し読み込むとき、イベントの前の最後のステップ サイズは、おそらく次の積分計算にとっても適切なスケールと考えられます。例「ballodeballode」を参照してください。

積分計算のタイム ステップ

質問

回答

積分計算で使われる最初のステップ サイズが大きすぎて、重要な動作を検出できません。

InitialStepプロパティを使えば、最初のステップ サイズを設定できます。積分器はまずこの値を使用し、必要に応じてさらに小さくします。

一定のステップ幅で積分計算できますか。

いいえ。

許容誤差とその他のオプション

質問

回答

RelTolAbsTolはどのようにして決めればよいでしょうか。

RelTol は相対許容誤差であり、解の正しい桁数を制御します。AbsTol は絶対許容誤差であり、真の解と計算された値との差を制御します。各タイム ステップで、解の i 番目の要素の誤差 e は、次の条件を満たす必要があります。

|e(i)| ≤max(RelTol*abs(y(i)),AbsTol(i))

したがって、おおまかに言うと、しきい値 AbsTol(i) より小さいものを除いて、すべての解の要素は RelTol 程度の正確な桁数が必要です。要素 y(i) が小さく、重要ではない場合でも、AbsTol(i) を小さくし、y(i) が正しい桁数を得るようにします。これにより、重要な要素を十分な精度で計算することができます。

コンピューターの計算精度と同等の解を求めたいのですが、RelToleps にしてはいけないのでしょうか。

マシンの精度に近い値を得ることはできますが、ぎりぎりまで近づけることはできません。ソルバーは連続関数を近似しようとするため、eps に近い値を RelTol として使用できません。許容範囲を eps に近い値にすると、マシンの計算は、すべての関数が不連続であるとみなすことになります。

ある 1 要素に対しては正確な解が必要ないことを設定するには、どのようにすればよいでしょうか。

その要素に対応する「絶対許容誤差」を大きくすることができます。許容誤差がその要素よりも大きい場合、その要素は正しい桁数には設定されません。他の要素を正確に計算するために、場合によってはこの要素にも正しい桁数を設定する必要がありますが、一般的に、このような操作は自動的に行われます。

方程式の他の形式

質問

回答

ライン手法で離散化した PDE (偏微分方程式) を扱うことができますか。

PDE を離散化すると ODE 方程式が生成されるので可能です。離散化の結果によっては、質量行列を含んだ形になるかもしれませんが、それには ODE ソルバーが使用できます。多くの場合、このシステムはスティッフになります。これは、PDE が放物型である場合や、流体内での化学反応のような非常に異なる時間スケールで生じる現象の場合に予想されます。このような場合、次の 4 つのソルバーのうちいずれかを使用してください。ode15s, ode23s, ode23t, ode23tb.

多くの方程式がある場合は JPattern プロパティを設定してください。これを設定することにより、非常にコストのかかる計算過程が成功するか失敗するかに影響を及ぼすこともあります。JPattern を使用する例として、「大規模でスティッフなスパース問題」を参照してください。システムがスティッフでない場合、あるいはそれほどスティッフではない場合には、ode23ode45 の方が、ode15sode23sode23tode23tb より有効です。

1次の放物型-楕円型偏微分方程式は、MATLAB PDEソルバー pdepe で直接解くことができます。詳細は、「偏微分方程式」を参照してください。

微分代数方程式 (DAE) 系を解けますか。

はい。ソルバー ode15sode23t は、M(t,y) が特異である M(t,y)y' = f(t,y) の形式の DAE を解くことができます。DAE 問題はインデックスが 1 でなければなりません。ode15i はインデックスが 1 の完全陰的 DAE f(t,y,y') = 0 を解くことができます。例は、「amp1dae」、「hb1daehb1dae」または「ihb1daeihb1dae」を参照してください。

サンプリングされたデータによるシステムを積分できますか。

直接的には無理です。内挿やその他のデータ近似の手法を使って、関数としてデータを表現しなければなりません。この関数については、滑らかさが重要です。スプラインのような区分的多項式で置き換えた場合、一見滑らかにみえますが、ソルバーにとっては不連続になっています。微分値が不連続になる点で、ソルバーは非常に細かいステップを使用します。データを表現するのに滑らかな関数を使用するか、感度的に落ちる低次のソルバー (ode23ode23sode23tode23tb) を使用してください。

最終時間での値はわかっているが、初期値がわからない場合どうしますか。

ODE 関連関数のすべてのソルバーは、時間に対して後ろ向きにも、前向きにも解くことができます。ソルバーの構文は [t,y] = ode45(odefun,[t0 tf],y0); となります。ここでは、t0 > tf も使用できます。

他の一般的な問題

質問

回答

計算結果が予想していたものと異なります。

予想の方が正しいのであれば、許容誤差を既定より小さくする必要があります。積分する時間間隔が長すぎる問題や、少し不安定な問題の解を正確に計算するには、より小さな相対許容誤差を使って計算する必要があります。

いずれかの時点において、求めた解の要素が、それに対応する絶対許容誤差よりも小さくなっていないかをチェックする必要があります。小さくなっている場合は、これらの要素について正しい数値をもつ桁数が確保されていないことになります。これらの要素については問題がなくても、正確に計算できていないために、それに依存する他の要素の精度が落ちてしまうこともあります。

プロットが十分滑らかではありません。

Refine の値を、既定値 (ode454、その他で 1) より大きい値にしてください。Refine の値を大きくするほど、出力点数が多くなります。Refine の値は、計算速度にはほとんど影響しません。

計算しながらプロットすると途中までは問題ないのですが、あるところで計算が止まってしまいます。

まず、計算が停止する点の近くで ODE の解が滑らかであるかを確認してください。滑らかでない場合は、その部分で小さな時間刻みが必要です。tspan の値を ODE 関数の滑らかな部分でいくつかに分割するとよいかもしれません。

関数が滑らかで計算が非常に小さなステップ幅で進んでいるのであれば、スティッフな問題に未対応のソルバーでスティッフな問題を解決しようとしているため、このようなことが起きている可能性があります。ode15sode23sode23tode23tb のいずれかのソルバーに変えてシミュレーションしてください。

非常に多くのタイム ステップを使うので、積分計算が非常に遅くなります。

はじめに、tspan に必要以上に多くの点を設定していないかどうかを調べてください。ソルバーは、滑らかな計算結果を出力するのに必要な数のタイム ステップで計算します。ODEの解が tspan に比べて非常に小さい時間スケールで変化している場合、多くのタイム ステップが必要になります。問題は積分計算が長いほど、解くことが難しくなります。tspan を小さく分割して計算してみてください。

tspan の上で、ODE 関数の解が目立って変化していないようであれば、その問題がスティッフであることが考えられます。ode15sode23sode23tode23tb のいずれかを使ってみてください。

最後に、ODE が効率よく記述されているかを確かめてください。ソルバーは ODE の導関数を何回も計算します。数値積分に要する時間は、ODE の計算に大きく依存します。複雑な計算になっている定数値を ODE の計算のたびに計算し直さずに、これらをグローバル値として保存しておくか、一度計算して入れ子関数に渡してみてください。

次の式を満たす時間 t で解が急変することがわかっています。

t0 ≤ t ≤ tf

しかし、積分のタイム ステップがそれを見逃しています。

時間 t で急変することがわかっている場合は、tspan[t0 t][t tf] の 2 つに分けて、2 回計算するとよいかもしれません。

微分方程式の係数または解に周期性がある場合、最大ステップサイズを周期長に設定し、ステップが周期長より大きくならないようにしてください。

この情報は役に立ちましたか?