Contex Contex_r Contr5 Contr5_r Contra Contra_r Contro Contro_r Contrp Contrp_r Dassl Dassl_r Radau Radau5 Radau5_r Radau_r Radaup Radaup_r Rodas Rodas_r Seulex Seulex_r

Dassl() Sub Dassl ( N As  Long, Res As  LongPtr, T As  Double, Y() As  Double, Yp() As  Double, Tout As  Double, RTol() As  Double, ATol() As  Double, Info As  Long, Optional Jac As  LongPtr = NullPtr, Optional Ml As  Long = -1, Optional Mu As  Long, Optional Neval As  Long, Optional Njac As  Long, Optional Nstep As  Long, Optional Mode As  Long, Optional ITstop As  Long, Optional Tstop As  Double, Optional Hmax As  Double, Optional H0 As  Double, Optional Maxord As  Long, Optional NonNeg As  Long, Optional NoInit As  Long  ) 微分代数方程式(DAE) (1〜5次 後退微分公式 (BDF)) 注 - 本プログラムは次バージョンで廃止予定です. 目的 本ルーチンは陰的な微分代数方程式(DAE) f(t, y, y') = 0, ただし t = t0 において y = y0 および y' = y0' の解を求める. ただし, t0, y0およびy0'は既知でそれぞれt, yおよびy'の初期値である. 上の方程式が連立微分代数方程式であれば, yおよびy'はベクトルで表される. 本ルーチンでは, 微分y'を1〜5次の後退微分公式(BDF)によって近似し, 各時間ステップにおいて得られた非線形連立方程式をニュートン法により解く. 引数 [in] N 微分方程式の数. (N >= 1) [in] Res 微分代数方程式の残差 delta = f(t, y, y') を求めるユーザー定義サブルーチンで, 次のように定義すること. Sub Res(N As Long, T As Double, Y() As Double, Yp() As Double, Delta() As Double, IRes As Long) Delta(i) = T, Y()およびYp()における残差 (i = 0 〜 N-1) End Sub ただし, Nは方程式の数, Y()はTにおける関数値, Yp()はTにおける微分係数である. IResは入力時常に0の整数フラグである. サブルーチンResでは, Y()の値が正しくなかった場合あるいは停止条件に会った場合にのみIResを変更しなけらばならない. 入力値が正しくなかった場合 IRes = -1 に設定すると, Dasslは IRes = -1 を避けて問題を解こうとする. IRes = -2 とすると, Dasslは呼び出し元プログラムに Info = 11 で戻る. 変数 N, T, Y() および Yp()をResの中で変更してはならない. [in,out] T 本ルーチンはTからToutまでの積分を行う. 積分を開始する点を与え, 最終ステップの最後の点が返される. 新たなToutにおける解を求めるために積分を継続することができる(インターバル・モードとよぶ). その場合, 2回目以降の呼び出しでは, Tは前回のToutに等しくなければならない. また, Toutまでの各中間ステップにおいて戻ることもできる(中間結果出力モードとよぶ). 本モードは解の挙動をみたい場合に使うとよい. モードはパラメータModeで指定できる. [in] 独立変数Tの初期値. [out] 独立変数Tの最終ステップの最後の点の値 (インターバル・モードでは通常Toutに等しい). 解がこの点まで正常に求められたことを示す. [in,out] Y() 配列 Y(LY - 1) (LY >= N) [in] Tの初期値における従属変数の初期値. [out] 最終のT(インターバル・モードにおいてはToutに等しい)において求められた解(の近似値). [in,out] Yp() 配列 Yp(LYyp - 1) (LYp >= N) [in] Tの初期値における従属変数の微分値の初期値. [out] 最終のT(インターバル・モードにおいてはToutに等しい)において求められた微分値(の近似値). [in] Tout 解を求めたい点を設定する. Tout = T とはできない. 積分を行うのはtについて前進方向(Tout > T)でも後退方向(Tout < T)でもよい. ただし, 初期化時でなければ積分方向を変更することはできない. 要求精度を満たすように自動的に選ばれたステップ幅を用いてTからToutに向かって解が求められる. 必要であれば, 各中間ステップにおける解とその微分係数を確認するために戻ることができる(中間結果出力モード). ただし, その場合であっても本来のとおりToutを指定しなければならない. [in] RTol RTol(LRTol - 1) (LRTol >= 1) (RTol()の全要素 >= 0) 求める解の精度を指定する相対誤差許容値. 本パラメータはスカラー(LRTol = 1)または配列(LRTol = N)のどちらでもよい. 配列は誤差テストをより細かく制御したい場合に使うことができる. LRTol = 2, ... または N-1の場合, LRTol = 1 とみなす. LRTol > Nの場合, LRTol = N とみなす. 許容値は各ステップにおける局所的な誤差テストに使われ, Y()の各要素がおおむね次式を満たすようにする (i = 0 〜 N-1).   abs(局所誤差) <= RTol(0)*abs(Y(i)) + ATol(0) (LRTolまたはLATol = 1 の場合)     または   abs(局所誤差) <= RTol(i)*abs(Y(i)) + ATol(i) (LRTolおよびLATol = N の場合) RTol(i) = 0 と設定するとその要素については純粋に絶対誤差テストとなる. [in] ATol ATol(LATol - 1) (LATol >= 1) (ATol()の全要素 >= 0) 求める解の精度を指定する絶対誤差許容値. 本パラメータはスカラー(LATol = 1)または配列(LATol = N)のどちらでもよい. 配列は誤差テストをより細かく制御したい場合に使うことができる. LATol = 2, ... または N-1の場合, LATol = 1 とみなす. LATol > Nの場合, LATol = N とみなす. 許容値は各ステップにおける局所的な誤差テストに使われ, Y()の各要素がおおむね次式を満たすようにする (i = 0 〜 N-1).   abs(局所誤差) <= RTol(0)*abs(Y(i)) + ATol(0) (LRTolまたはLATol = 1 の場合)     または   abs(局所誤差) <= RTol(i)*abs(Y(i)) + ATol(i) (LRTolおよびLATol = N の場合) ATol(i) = 0 と設定するとその要素については純粋に相対誤差テストとなる. [in,out] Info [in] 制御コード. = 0: 問題の最初の呼び出し時(新たに問題を開始する場合)には Info = 0 と設定せよ. 初期化が行われ, 新たな問題についての計算が開始される. = 1, 2, 11〜13: Info = 1, 2, 11〜13 で戻った場合, 計算を継続するために, Infoの値を変えずに再呼び出しすることができる(下記説明を参照せよ). [out] リターンコード. ユーザーはこの値をチェックし, Info = 1, 2, 11〜13 であれば必要により次のアクションとして本ルーチンを再度呼び出すことができる. = -1: パラメータ N の誤り. (N < 1) = -3: パラメータ T の誤り. (継続呼び出しで T = Tout または T <> 前回のTout) = -4: パラメータ Y() の誤り. = -5: パラメータ Yp() の誤り. = -6: パラメータ Tout の誤り. (積分の方向を変更した) = -7: パラメータ Rtol() の誤り. (RTol(i) < 0) = -8: パラメータ Atol() の誤り. (ATol(i) < 0) = -9: パラメータ Info の誤り. (Info <> 0, 1, 2, 11〜13) = -18: パラメータ Tstop の誤り. (TstopがTまたはToutより後ろ) = 1: 正常終了 (T = Tout). Toutを再設定して再呼び出しすることができる. ToutはTと異なる値であること. = 2: 中間結果出力モードのため戻った (まだ Toutに達していない). Tout方向の次のステップを続行するために再呼び出しすることができる. = 11: 最大ステップ数(500)を超えた. 再呼び出しして続行することができる. さらに500ステップ続けることができる. = 12: 誤差許容値が厳しすぎる. 再呼び出しして, 自動的に調節された許容値を使って続行することができる. 許容値を自分で変更してから再呼び出しすることもできる. = 13: 求められた解が0であるため, 純粋な相対誤差テスト(ATol = 0)ができない. Atolを正の値にして再呼び出しして続行することができる. = 14: 誤差テストの失敗を繰り返した. = 15: 修正子の収束テストの失敗を繰り返した. = 16: 偏微分行列が特異である. = 17: 最後の試行ステップにおいて収束テストの失敗を繰り返した. = 18: IRes = -1 のため修正子が収束しない. = 19: IRes = -2 のため終了した. = 20: 微分の初期値の計算に失敗した. = 21: 無限ループが検出された. [in] Jac (省略可) 偏微分行列を求めるユーザー定義サブルーチンで, 次のように定義すること. (省略時 = NullPtr) Sub Jac(N As Long, T As Double, Y() As Double, Yp() As Double, Ypd() As Double, Cj As Double) dfi/dyj + Cj*dfi/dyj' (i = 0〜N-1, j = 0〜N-1) をYpd()に設定. End Sub TおよびY()における偏微分行列の計算値をYpd()に設定すること. Ml = N の場合 Ypd()は通常のN×N行列, 0 <= Ml < N の場合(2Ml+Mu+1×N)帯行列形式である. 帯行列形式の詳細は下記を参照のこと. Ypd()の全ての要素はJacの呼び出し前に0に初期化されるので, 0でない要素だけ設定すればよい. 変数N, T, Y(), Yp() および CjをJacの中で変更してはならない. これを省略した場合(Djac = NullPtrの場合), 偏微分行列は差分近似により計算される. [in] Ml (省略可) Ml = N の場合, 偏微分行列は通常のN×N行列として扱われる. 0 <= Ml < N であれば帯行列形式で扱う(2*Ml + Mu < N のときに効果がある). その場合, Mlは下帯幅を表す. (Ml >= 0, Ml <= N) (省略時 = N) [in] Mu (省略可) 偏微分行列を帯行列として扱う場合に上帯幅を表す. Ml = N の場合には参照されない. (Mu >= 0, Mu < N) (省略時 = 0) [out] Neval (省略可) 関数評価回数. (偏微分行列の計算のためのものは含まない) [out] Njac (省略可) 偏微分行列評価回数. (有限差分の場合を含む) [out] Nstep (省略可) ステップ数. [in] Mode (省略可) 動作モード. (省略時 = 0) = 0: Toutで戻る (インターバル・モード). = 1: ステップごとに戻る (中間結果出力モード) (他の値であれば省略時の既定値とみなす) [in] ITstop (省略可) 本ルーチンはToutを超えて積分を行い内挿によりToutにおける値を求めることがあるが, ある点において不連続になったり, その点を超えて関数あるいは微分が定義されていなかったりして, その点を超えて積分を行うのが許されないことがある. そのような場合, ストップポイント Tstop を設定することができ, その点を超えての積分を行わないようにすることができる. (省略時 = 0) = 0: ストップポイントを設定しない = 1: Tstopをストップポイントとして設定する (上記以外の値であれば省略時の既定値とみなす) [in] Tstop (ITstop = 0 のとき省略可) ストップポイントの値. [in] Hmax (省略可) 最大ステップ幅. (Hmax > 0) (省略時 = 自動選択する) (Hmax <= 0 であれば省略時の既定値とみなす) [in] H0 (省略可) ステップ幅の初期値. (H0 > 0) (省略時 = 自動選択する) (H0 <= 0 であれば省略時の既定値とみなす) [in] Maxord (省略可) 最大次数. (1 <= Maxord <= 5) (省略時 = 5) (Maxord <= 0 または Maxord >= 6 であれば省略時の既定値とみなす) [in] NonNeg (省略可) 方程式の解が常に非負であるかどうかを指定する. (省略時 = 0) = 0: 非負の制約を課することなく問題を解く. = 1: 特別に非負の制約を課して問題を解く. (上記以外の値であれば省略時の既定値とみなす) [in] NoInit (省略可) 初期値について矛盾がない (f(t, y, y') = 0 が成り立っている) かどうかを指定する. (省略時 = 0) = 0: t, y および y'の初期値は矛盾がない. = 1: 初期値が正確にわからないためプログラムでy'の計算を試行する. Yp()に初期近似を設定しておくこと. もし全く不明であれば, Yp()に0を設定しておく. (上記以外の値であれば省略時の既定値とみなす) 詳細 次の例は, N = 6, Ml = 2, Mu = 1 の場合の, 本ルーチンで使われる2*Ml+Mu+1行×N列の帯行列形式を表す. 一般行列形式: a11 a12 0 0 0 0 a21 a22 a23 0 0 0 a31 a32 a33 a34 0 0 0 a42 a43 a44 a45 0 0 0 a53 a54 a55 a56 0 0 0 a64 a65 a66 帯行列形式: * * * + + + * * + + + + * a12 a23 a34 a45 a56 a11 a22 a33 a44 a55 a66 a21 a32 a43 a54 a65 * a31 a42 a53 a64 * * +で示された最初の2行は作業領域として使用される. *で示された配列要素は使用されない. 出典 SLATEC (DASSL) 使用例 次の微分代数方程式を解く(Robertsonの問題). -0.04*y1 + 10000*y2*y3 - y1' = 0 0.04*y1 - 10000*y2*y3 - 30000000*y2^2 - y2' = 0 y1 + y2 + y3 -1 = 0 (t = 0 において y1 = 1, y2 = 0, y3 = 0, y1' = y2' = y3' = 0) Sub R2(N As Long, T As Double, Y() As Double, Yp() As Double, Delta() As Double, Ires As Long) Delta(0) = -0.04 * Y(0) + 10000 * Y(1) * Y(2) - Yp(0) Delta(1) = 0.04 * Y(0) - 10000 * Y(1) * Y(2) - 30000000 * Y(1) ^ 2 - Yp(1) Delta(2) = Y(0) + Y(1) + Y(2) - 1 End Sub Sub Ex_Dassl() Const N = 3 Dim T As Double, Y(N - 1) As Double, Yp(N - 1) As Double, Tout As Double Dim RTol(0) As Double, ATol(0) As Double, ITol As Long Dim Tend As Double, Info As Long '--- Parameters RTol(0) = 0.0000001: ATol(0) = 0.0000001: ITol = 0 '--- Initial values T = 0 Y(0) = 1: Y(1) = 0: Y(2) = 0 Yp(0) = 0: Yp(1) = 0: Yp(2) = 0 '--- Start calculation Tend = 10 Info = 0 Do Tout = T + 1 Call Dassl(N, AddressOf R2, T, Y(), Yp(), Tout, RTol(), ATol(), Info) If T >= Tout Then Debug.Print T, Y(0), Y(1), Y(2) End If If T >= Tend Then Exit Do Loop While Info = 1 If Info <> 1 Then Debug.Print "Error in Dassl: Info =", Info End Sub 実行結果 1 0.966459990474207 3.07463058949132E-05 3.35092632198979E-02 2 0.941609675690738 2.7017862636127E-05 5.83633064466271E-02 3 0.921884807943058 2.43833786001128E-05 7.80908086783409E-02 4 0.905518993577723 2.24047935986599E-05 9.44586016286797E-02 5 0.891518109563595 2.08527025402239E-05 0.108461037733863 6 0.87927028562438 1.95949526270386E-05 0.120710119422982 7 0.868373377809956 1.85498046053848E-05 0.13160807238542 8 0.85854920043329 1.76638124059415E-05 0.141433135754299 9 0.849597711635835 1.69004734510369E-05 0.150385387890728 10 0.841370267657359 1.62339372132848E-05 0.158613498405457