[doc][equadif + polynomials] Unfold synopses + fix + improvements
[scilab.git] / scilab / modules / differential_equations / help / ja_JP / dae.xml
1 <?xml version="1.0" encoding="UTF-8"?>
2
3 <!--
4  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
5  * Copyright (C) 2008 - INRIA
6  * ...
7  *
8  * Copyright (C) 2012 - 2016 - Scilab Enterprises
9  *
10  * This file is hereby licensed under the terms of the GNU GPL v2.0,
11  * pursuant to article 5.3.4 of the CeCILL v.2.1.
12  * This file was originally licensed under the terms of the CeCILL v2.1,
13  * and continues to be available under such terms.
14  * For more information, see the COPYING file which you should have received
15  * along with this program.
16  *
17  -->
18
19 <refentry xmlns="http://docbook.org/ns/docbook" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:svg="http://www.w3.org/2000/svg" xmlns:ns3="http://www.w3.org/1999/xhtml" xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:db="http://docbook.org/ns/docbook" xmlns:scilab="http://www.scilab.org" xml:id="dae" xml:lang="ja">
20     <refnamediv>
21         <refname>dae</refname>
22         <refpurpose>微分代数方程式ソルバ</refpurpose>
23     </refnamediv>
24     <refsynopsisdiv>
25         <title>呼び出し手順</title>
26         <synopsis> y = dae(initial, t0, t, res)
27             [y [,hd]] = dae(initial, t0, t [[,rtol], atol], res [,jac] [,hd])
28             [y, rd] = dae("root", initial, t0, t, res, ng, surface)
29             [y, rd [,hd]] = dae("root", initial, t0, t [[,rtol], atol], res [,jac], ng, surface [,hd])
30             [y, rd] = dae("root2", initial, t0, t, res, ng, surface)
31             [y, rd [,hd]] = dae("root2", initial, t0, t [[,rtol], atol], res [,jac], ng, surface [, psol, pjac] [, hd])
32         </synopsis>
33     </refsynopsisdiv>
34     <refsection>
35         <title>引数</title>
36         <variablelist>
37             <varlistentry>
38                 <term>initial</term>
39                 <listitem>
40                     <para>
41                         列ベクトル. <literal>x0 </literal>または
42                         <literal>[x0;xdot0]</literal> とします.
43                         ただし, <literal>x0</literal>は
44                         初期時刻<literal>t0</literal>における状態量の値,
45                         <literal>ydot0</literal>は状態量の微分の初期値
46                         またはその推定値(下記参照)です.
47                     </para>
48                 </listitem>
49             </varlistentry>
50             <varlistentry>
51                 <term>t0</term>
52                 <listitem>
53                     <para>実数, 初期時刻.</para>
54                 </listitem>
55             </varlistentry>
56             <varlistentry>
57                 <term>t</term>
58                 <listitem>
59                     <para>実数のスカラーまたはベクトル. 解を計算する時刻を指定します.
60                         <literal>
61                             <link linkend="daeoptions">%DAEOPTIONS</link>(2)=1
62                         </literal>
63                         と設定することにより
64                         DAEの各ステップで解が得られることに注意してください.
65                     </para>
66                 </listitem>
67             </varlistentry>
68             <varlistentry>
69                 <term>atol</term>
70                 <listitem>
71                     <para>
72                         実数のスカラーまたは<literal>x0</literal>と同じ大きさの
73                         列ベクトル. 解の絶対許容誤差.
74                         <literal>atol</literal>がベクトルの場合,
75                         状態量の各要素毎に許容誤差が指定されます.
76                     </para>
77                 </listitem>
78             </varlistentry>
79             <varlistentry>
80                 <term>rtol</term>
81                 <listitem>
82                     <para>
83                         実数のスカラーまたは<literal>x0</literal>と同じ大きさの
84                         列ベクトル. 解の相対許容誤差.
85                         <literal>rtol</literal>がベクトルの場合,
86                         状態量の各要素毎に許容誤差が指定されます.
87                     </para>
88                 </listitem>
89             </varlistentry>
90             <varlistentry>
91                 <term>res</term>
92                 <listitem>
93                     <para>
94                         <link linkend="external" role="" version="">外部ルーチン</link>.
95                         <literal>g(t,y,ydot)</literal>の値を計算します. 以下のようになります
96                     </para>
97                     <variablelist>
98                         <varlistentry>
99                             <term>Scilab関数</term>
100                             <listitem>
101                                 <para>この場合, その呼出し手順を
102                                     <literal>[r,ires]=res(t,x,xdot)</literal> とする
103                                     必要があり,<literal>res</literal> は
104                                     残差<literal>r=g(t,x,xdot)</literal> とエラーフラグ
105                                     <literal>ires</literal>を返す必要があります.
106                                     <literal>res</literal>が<literal>r</literal>の計算に
107                                     成功した場合には<literal>ires = 0</literal>,
108                                     <literal>(t,x,xdot)</literal>のローカルな残差が定義されない
109                                     場合には <literal>=-1</literal>,
110                                     パラメータが許容範囲外の場合は <literal>=-2</literal> となります.
111                                 </para>
112                             </listitem>
113                         </varlistentry>
114                         <varlistentry>
115                             <term>リスト</term>
116                             <listitem>
117                                 <para>外部ルーチンのこの形式は
118                                     関数にパラメータを指定する際に使用されます.
119                                     以下のような形式とします:
120                                 </para>
121                                 <screen><![CDATA[
122 list(res,p1,p2,...)
123 ]]></screen>
124                                 <para>
125                                     ただし,ここで関数<literal>res</literal>の呼び出し手順は以下のようになります
126                                 </para>
127                                 <screen><![CDATA[
128 r = res(t,y,ydot,p1,p2,...)
129 ]]></screen>
130                                 <para>
131                                     この場合も<literal>res</literal> は
132                                     <literal>(t,x,xdot,x1,x2,...)</literal>の関数として残差の値を返し,
133                                     p1,p2,... は関数のパラメータです.
134                                 </para>
135                             </listitem>
136                         </varlistentry>
137                         <varlistentry>
138                             <term>文字列</term>
139                             <listitem>
140                                 <para>CまたはFortranルーチンの名前を指します.
141                                     &lt;r_name&gt; が指定した名前であると仮定します.
142                                 </para>
143                                 <itemizedlist>
144                                     <listitem>
145                                         <para>Fortranの呼出し手順は以下のようにします</para>
146                                         <para>
147                                             <literal>&lt;r_name&gt;(t,x,xdot,res,ires,rpar,ipar)</literal>
148                                         </para>
149                                         <para>
150                                             <literal>double precision
151                                                 t,x(*),xdot(*),res(*),rpar(*)
152                                             </literal>
153                                         </para>
154                                         <para>
155                                             <literal>integer ires,ipar(*)</literal>
156                                         </para>
157                                     </listitem>
158                                     <listitem>
159                                         <para>Cの呼出し手順は以下のようにします</para>
160                                         <para>C2F(&lt;r_name&gt;)(double *t, double *x, double
161                                             *xdot, double *res, integer *ires, double *rpar, integer
162                                             *ipar)
163                                         </para>
164                                     </listitem>
165                                 </itemizedlist>
166                                 <para>ただし,</para>
167                                 <itemizedlist>
168                                     <listitem>
169                                         <para>
170                                             <literal>t</literal> はカレントの時刻
171                                         </para>
172                                     </listitem>
173                                     <listitem>
174                                         <para>
175                                             <literal>x</literal>は状態量の配列
176                                         </para>
177                                     </listitem>
178                                     <listitem>
179                                         <para>
180                                             <literal>xdot</literal> は状態量の微分の配列
181                                         </para>
182                                     </listitem>
183                                     <listitem>
184                                         <para>res は残差の配列</para>
185                                     </listitem>
186                                     <listitem>
187                                         <para>
188                                             <literal>ires</literal>は実行インジケータ
189                                         </para>
190                                     </listitem>
191                                     <listitem>
192                                         <para>
193                                             <literal>rpar</literal> は浮動小数点数の
194                                             パラメータ値の配列で, 必要ですが <literal>dae</literal> 関数により
195                                             設定できないもの
196                                         </para>
197                                     </listitem>
198                                     <listitem>
199                                         <para>
200                                             <literal>ipar</literal> は整数の
201                                             パラメータ値の配列で, 必要ですが <literal>dae</literal> 関数により
202                                             設定できないもの
203                                         </para>
204                                     </listitem>
205                                 </itemizedlist>
206                             </listitem>
207                         </varlistentry>
208                     </variablelist>
209                 </listitem>
210             </varlistentry>
211             <varlistentry>
212                 <term>jac</term>
213                 <listitem>
214                     <para>
215                         <link linkend="external">外部ルーチン</link>.
216                         指定したパラメータの値 <literal>cj</literal>を用いて
217                         <literal>dg/dx+cj*dg/dxdot</literal>の値を計算します. 以下のようになります
218                     </para>
219                     <variablelist>
220                         <varlistentry>
221                             <term>Scilab関数</term>
222                             <listitem>
223                                 <para>呼出し手順を
224                                     <literal>r=jac(t,x,xdot,cj)</literal> とする
225                                     必要があり,
226                                     <literal>jac</literal>関数は
227                                     <literal>r=dg(t,x,xdot)/dy+cj*dg(t,x,xdot)/dxdot</literal>
228                                     を返す必要があります.
229                                     ただし,<literal>cj</literal>は実数スカラーです.
230                                 </para>
231                             </listitem>
232                         </varlistentry>
233                         <varlistentry>
234                             <term>リスト</term>
235                             <listitem>
236                                 <para>外部ルーチンのこの形式は
237                                     関数にパラメータを指定する際に使用されます.
238                                     以下のような形式とします:
239                                 </para>
240                                 <screen><![CDATA[
241 list(jac,p1,p2,...)
242 ]]></screen>
243                                 <para>
244                                     ただしこの場合の関数<literal>jac</literal>の呼び出し手順は
245                                     以下となります
246                                 </para>
247                                 <screen><![CDATA[
248 r = jac(t,x,xdot,p1,p2,...)
249 ]]></screen>
250                                 <para>
251                                     この場合でも<literal>jac</literal> は,
252                                     <literal>(t,x,xdot,cj,p1,p2,...)</literal>の関数として
253                                     <literal>dg/dx+cj*dg/dxdot</literal>を返します.
254                                 </para>
255                             </listitem>
256                         </varlistentry>
257                         <varlistentry>
258                             <term>文字列</term>
259                             <listitem>
260                                 <para>CまたはFortranルーチンの名前を指します.
261                                     &lt;j_name&gt; が指定した名前であると仮定します.
262                                 </para>
263                                 <itemizedlist>
264                                     <listitem>
265                                         <para>Fortranの呼出し手順は以下のようにします</para>
266                                         <para>
267                                             <literal>&lt;j_name&gt;(t, x, xdot, r, cj, ires,
268                                                 rpar, ipar)
269                                             </literal>
270                                         </para>
271                                         <para>
272                                             <literal>double precision t, x(*), xdot(*), r(*),
273                                                 ci, rpar(*)
274                                             </literal>
275                                         </para>
276                                         <para>
277                                             <literal>integer ires, ipar(*)</literal>
278                                         </para>
279                                     </listitem>
280                                     <listitem>
281                                         <para>Cの呼出し手順は以下のようにします</para>
282                                         <para>C2F(&lt;j_name&gt;)(double *t, double *x, double
283                                             *xdot, double *r, double *cj,
284                                         </para>
285                                         <para>integer *ires, double *rpar, integer *ipar)</para>
286                                     </listitem>
287                                 </itemizedlist>
288                                 <para>
289                                     ただし <literal>t</literal>, x, xdot, ires, rpar, ipar
290                                     は上記と似た定義を有し, r は配列の結果です.
291                                 </para>
292                             </listitem>
293                         </varlistentry>
294                     </variablelist>
295                 </listitem>
296             </varlistentry>
297             <varlistentry>
298                 <term>surface</term>
299                 <listitem>
300                     <para>
301                         <link linkend="external">外部ルーチン</link>.
302                         <literal>ng</literal>個の要素を有する列ベクトル
303                         <literal>surface(t,x)</literal>の値を計算します.
304                         各要素は面(surface)を定義します.
305                     </para>
306                     <variablelist>
307                         <varlistentry>
308                             <term>Scilab関数</term>
309                             <listitem>
310                                 <para>
311                                     その呼び出し手順は<literal>r=surface(t,x)</literal>とする
312                                     必要があり,この関数は <literal>ng</literal>個の要素を有する
313                                     ベクトルを返す必要があります.
314                                 </para>
315                             </listitem>
316                         </varlistentry>
317                         <varlistentry>
318                             <term>リスト</term>
319                             <listitem>
320                                 <para>外部ルーチンのこの形式は
321                                     関数にパラメータを指定する際に使用されます.
322                                     以下のような形式とします:
323                                 </para>
324                                 <screen><![CDATA[
325 list(surface,p1,p2,...)
326 ]]></screen>
327                                 <para>
328                                     ただしこの場合の関数<literal>surface</literal>の呼び出し手順は
329                                     以下となります
330                                 </para>
331                                 <screen><![CDATA[
332 r = surface(t,x,p1,p2,...)
333 ]]></screen>
334                             </listitem>
335                         </varlistentry>
336                         <varlistentry>
337                             <term>文字列</term>
338                             <listitem>
339                                 <para>CまたはFortranルーチンの名前を指します.
340                                     &lt;s_name&gt; が指定した名前であると仮定します,
341                                 </para>
342                                 <itemizedlist>
343                                     <listitem>
344                                         <para>Fortranの呼出し手順は以下のようにします</para>
345                                         <para>
346                                             <literal>&lt;r_name&gt;(nx, t, x, ng, r, rpar,
347                                                 ipar)
348                                             </literal>
349                                         </para>
350                                         <para>
351                                             <literal>double precision t, x(*), r(*),
352                                                 rpar(*)
353                                             </literal>
354                                         </para>
355                                         <para>
356                                             <literal>integer nx, ng,ipar(*)</literal>
357                                         </para>
358                                     </listitem>
359                                     <listitem>
360                                         <para>Cの呼出し手順は以下のようにします</para>
361                                         <para>C2F(&lt;r_name&gt;)(double *t, double *x, double
362                                             *xdot, double *r, double *cj,
363                                         </para>
364                                         <para>integer *ires, double *rpar, integer *ipar)</para>
365                                     </listitem>
366                                 </itemizedlist>
367                                 <para>
368                                     ただし,<literal>t</literal>, x, rpar, ipar は
369                                     上記と同じ定義を有し,<literal>ng </literal> は
370                                     surfacesの数, <literal>nx</literal> は状態量の次元,
371                                     r は結果の配列です.
372                                 </para>
373                             </listitem>
374                         </varlistentry>
375                     </variablelist>
376                 </listitem>
377             </varlistentry>
378             <varlistentry>
379                 <term>rd</term>
380                 <listitem>
381                     <para>
382                         2個のエントリ <literal>[times num]</literal>を有するベクトル.
383                         <literal>times</literal> は面が交差した時刻の値,
384                         <literal>num</literal> は交差した面の数です
385                     </para>
386                 </listitem>
387             </varlistentry>
388             <varlistentry>
389                 <term>pjac</term>
390                 <listitem>
391                     <para>
392                         外部 (関数, リストまたは文字列).
393                         <link linkend="daskr">daskr</link>参照.
394                     </para>
395                 </listitem>
396             </varlistentry>
397             <varlistentry>
398                 <term>psol</term>
399                 <listitem>
400                     <para>
401                         外部 (関数, リストまたは文字列).
402                         <link linkend="daskr">daskr</link>参照.
403                     </para>
404                 </listitem>
405             </varlistentry>
406             <varlistentry>
407                 <term>hd</term>
408                 <listitem>
409                     <para>実数のベクトル,
410                         <literal>dae</literal> コンテキストを保持する出力.
411                         (ホットスタートで)積分を再開するための入力引数として使用可能です.
412                     </para>
413                 </listitem>
414             </varlistentry>
415             <varlistentry>
416                 <term>y</term>
417                 <listitem>
418                     <para>
419                         実数の行列.
420                         <literal>
421                             <link linkend="daeoptions">%DAEOPTIONS</link>(2)=1
422                         </literal>
423                         の場合, 各列
424                         はベクトル <literal>[t;x(t);xdot(t)]</literal> です.
425                         ただし, <literal>t</literal> は解が計算された時刻です.
426                         それ以外の場合, <literal>y</literal> はベクトル
427                         <literal>[x(t);xdot(t)]</literal>です.
428                     </para>
429                 </listitem>
430             </varlistentry>
431         </variablelist>
432     </refsection>
433     <refsection>
434         <title>説明</title>
435         <para>
436             <literal>dae</literal> 関数は,
437             陰的な微分方程式積分用に設計された
438             <link linkend="dassl">dassl</link> および <link linkend="dasrt">dasrt</link>
439             関数の上位に構築されたゲートウエイです.
440         </para>
441         <para>
442             オプション <literal>"root"</literal> は
443             <link linkend="dasrt">dasrt</link> ルーチンをコールし,
444             <literal>"root2"</literal> は <link linkend="dasrt">daskr</link>をコールします.
445         </para>
446         <screen><![CDATA[
447 g(t,x,xdot) = 0
448 x(t0) = x0  および   xdot(t0)=xdot0
449 ]]></screen>
450         <para>
451             <emphasis>initial</emphasis>引数に<literal>xdot0</literal> が指定されない場合,
452             <literal>dae</literal>関数は <literal>g(t,x0,xdot0)=0</literal>
453             を解くことによりこれを計算しようとします.
454         </para>
455         <para>
456             <literal>xdot0</literal> が <emphasis>initial</emphasis>引数で指定された場合,
457             <literal>g(t,x0,xdot0)=0</literal>を満たす互換性がある微係数または近似値の
458             どちらかとなります.
459             後者の場合,
460             <literal>
461                 <link linkend="daeoptions">%DAEOPTIONS</link>(7)
462             </literal>
463             を 1 に設定する必要があります.
464         </para>
465         <para>Scilab および C でコーディングされた外部ルーチンを使用する詳細な例については,
466             <literal>modules/differential_equations/tests/unit_tests/dassldasrt.tst</literal>
467             および
468             <literal>modules/differential_equations/tests/unit_tests/daskr.tst</literal>
469             で提供されています.
470         </para>
471     </refsection>
472     <refsection>
473         <title>例</title>
474         <para>
475             例 #1: dassl (求解なし)
476         </para>
477         <programlisting role="example"><![CDATA[
478 // Scilabコードを使用する例
479 function [r,ires]=chemres(t,y,yd)
480     r(1) = -0.04*y(1) + 1d4*y(2)*y(3) - yd(1);
481     r(2) =  0.04*y(1) - 1d4*y(2)*y(3) - 3d7*y(2)*y(2) - yd(2);
482     r(3) =       y(1) +     y(2)      + y(3)-1;
483     ires =  0;
484 endfunction
485 function pd=chemjac(x,y,yd,cj)
486     pd=[-0.04-cj , 1d4*y(3)               , 1d4*y(2);
487          0.04    ,-1d4*y(3)-2*3d7*y(2)-cj ,-1d4*y(2);
488          1       , 1                      , 1       ]
489 endfunction
490
491 x0 = [1; 0; 0];
492 xd0 = [-0.04; 0.04; 0];
493 t = [1.d-5:0.02:.4, 0.41:.1:4, 40, 400, 4000, 40000, 4d5, 4d6, 4d7, 4d8, 4d9, 4d10];
494
495 y = dae([x0,xd0],0,t,chemres);// 指定した観測時刻での点を返す
496
497 %DAEOPTIONS = list([],1,[],[],[],0,0); // dae メッシュ点を返すかどうかを確認
498 y = dae([x0,xd0],0,4d10,chemres); // ヤコビアン指定なし
499 y = dae([x0,xd0],0,4d10,chemres,chemjac); // ヤコビアン指定あり
500      ]]></programlisting>
501         <para>
502             例 #2: dasrt ("root")
503         </para>
504         <programlisting role="example"><![CDATA[
505 // C コードの例 (Cコンパイラが必要) --------------------------------------------------
506 bOK = haveacompiler();
507 if bOK <> %t
508     [btn] = messagebox(["You need a C compiler for this example."; "Execution of this example is canceled."], "Software problem", 'info');
509     return
510 end
511
512 //-1- Cコードを TMPDIR に作成 - Vanderpol 方程式, 陰形式
513 code=['#include <math.h>'
514       'void res22(double *t,double *y,double *yd,double *res,int *ires,double *rpar,int *ipar)'
515       '{res[0] = yd[0] - y[1];'
516       ' res[1] = yd[1] - (100.0*(1.0 - y[0]*y[0])*y[1] - y[0]);}'
517       ' '
518       'void jac22(double *t,double *y,double *yd,double *pd,double *cj,double *rpar,int *ipar)'
519       '{pd[0]=*cj - 0.0;'
520       ' pd[1]=    - (-200.0*y[0]*y[1] - 1.0);'
521       ' pd[2]=    - 1.0;'
522       ' pd[3]=*cj - (100.0*(1.0 - y[0]*y[0]));}'
523       ' '
524       'void gr22(int *neq, double *t, double *y, int *ng, double *groot, double *rpar, int *ipar)'
525       '{ groot[0] = y[0];}']
526 mputl(code,TMPDIR+'/t22.c')
527 //-2- コンパイルの後,ロード
528 ilib_for_link(['res22' 'jac22' 'gr22'],'t22.c',[],'c',[],TMPDIR+'/t22loader.sce');
529 exec(TMPDIR+'/t22loader.sce')
530 //-3- 実行
531 rtol = [1.d-6;1.d-6];atol=[1.d-6;1.d-4];
532 t0 = 0; y0 = [2;0]; y0d = [0;-2]; t = [20:20:200]; ng = 1;
533 //簡単なシミュレーション
534 t = 0:0.003:300;
535 yy = dae([y0,y0d],t0,t,atol,rtol,'res22','jac22');
536 clf(); plot(yy(1,:),yy(2,:))
537 //find first point where yy(1)=0
538 [yy,nn,hotd] = dae("root",[y0,y0d],t0,300,atol,rtol,'res22','jac22',ng,'gr22');
539 plot(yy(1,1),yy(2,1),'r+')
540 xstring(yy(1,1)+0.1,yy(2,1),string(nn(1)))
541
542 //次の点までホットリスタート
543 t01 = nn(1);
544 [pp,qq] = size(yy);
545 y01 = yy(2:3,qq); y0d1 = yy(3:4,qq);
546 [yy,nn,hotd] = dae("root",[y01,y0d1],t01,300,atol,rtol,'res22','jac22',ng,'gr22',hotd);
547 plot(yy(1,1),yy(2,1),'r+')
548 xstring(yy(1,1)+0.1,yy(2,1),string(nn(1)));
549 cd(previous_dir);
550      ]]></programlisting>
551         <scilab:image localized="false"><![CDATA[
552 code = ['#include <math.h>'
553       'void res22(double *t, double *y, double *yd, double *res, int *ires, double *rpar, int *ipar)'
554       '{res[0] = yd[0] - y[1];'
555       ' res[1] = yd[1] - (100.0*(1.0 - y[0]*y[0])*y[1] - y[0]);}'
556       ' '
557       'void jac22(double *t, double *y, double *yd, double *pd, double *cj, double *rpar, int *ipar)'
558       '{pd[0] = *cj - 0.0;'
559       ' pd[1] =     - (-200.0*y[0]*y[1] - 1.0);'
560       ' pd[2] =     - 1.0;'
561       ' pd[3] = *cj - (100.0*(1.0 - y[0]*y[0]));}'
562       ' '
563       'void gr22(int *neq, double *t, double *y, int *ng, double *groot, double *rpar, int *ipar)'
564       '{ groot[0] = y[0];}']
565 previous_dir = pwd();
566 cd TMPDIR;
567 mputl(code, 't22.c')
568 ilib_for_link(['res22' 'jac22' 'gr22'], 't22.c', [], 'c', [], 't22loader.sce');
569 exec('t22loader.sce')
570 rtol = [1.d-6; 1.d-6];
571 atol = [1.d-6; 1.d-4];
572 t0 = 0; t = [20:20:200];
573 y0 = [2; 0]; y0d = [0; -2];
574 ng = 1;
575 t = 0:0.003:300;
576 yy = dae([y0, y0d], t0, t, atol, rtol, 'res22', 'jac22');
577 clf(); plot(yy(1, :), yy(2, :))
578 [yy, nn, hotd] = dae("root", [y0, y0d], t0, 300, atol, rtol, 'res22', 'jac22', ng, 'gr22');
579 plot(yy(1, 1), yy(2, 1), 'r+')
580 xstring(yy(1, 1)+0.1, yy(2, 1), string(nn(1)));
581 t01 = nn(1);
582 [pp, qq] = size(yy);
583 y01 = yy(2:3, qq);
584 y0d1 = yy(3:4, qq);
585 [yy, nn, hotd] = dae("root", [y01, y0d1], t01, 300, atol, rtol, 'res22', 'jac22', ng, 'gr22', hotd);
586 plot(yy(1, 1), yy(2, 1), 'r+')
587 xstring(yy(1, 1)+0.1, yy(2, 1), string(nn(1)));
588 cd(previous_dir);
589  ]]></scilab:image>
590         <para>
591             例 #3: daskr ("root2"), デフォルトの 'psol' および 'pjac' ルーチンを使用
592         </para>
593         <programlisting role="example"><![CDATA[
594 // Example with C code (C compiler needed)
595 //--------------------------------------------------
596 bOK = haveacompiler();
597 if bOK <> %t
598     [btn] = messagebox(["You need a C compiler for this example."; "Execution of this example is canceled."], "Software problem", 'info');
599     return
600 end
601
602 //-1- Create the C codes in TMPDIR - Vanderpol equation, implicit form
603 code = ['#include <math.h>'
604       'void res22(double *t, double *y, double *yd, double *res, int *ires, double *rpar, int *ipar)'
605       '{res[0] = yd[0] - y[1];'
606       ' res[1] = yd[1] - (100.0*(1.0 - y[0]*y[0])*y[1] - y[0]);}'
607       ' '
608       'void jac22(double *t, double *y, double *yd, double *pd, double *cj, double *rpar, int *ipar)'
609       '{pd[0] = *cj - 0.0;'
610       ' pd[1] =     - (-200.0*y[0]*y[1] - 1.0);'
611       ' pd[2] =     - 1.0;'
612       ' pd[3] = *cj - (100.0*(1.0 - y[0]*y[0]));}'
613       ' '
614       'void gr22(int *neq, double *t, double *y, int *ng, double *groot, double *rpar, int *ipar)'
615       '{ groot[0] = y[0];}']
616 previous_dir = pwd();
617 cd TMPDIR;
618 mputl(code, 't22.c')
619
620 //-2- Compile and load them
621 ilib_for_link(['res22' 'jac22' 'gr22'], 't22.c', [], 'c', [], 't22loader.sce');
622 exec('t22loader.sce')
623
624 //-3- Run
625 rtol = [1.d-6; 1.d-6];
626 atol = [1.d-6; 1.d-4];
627 t0 = 0; t = [20:20:200];
628 y0 = [2; 0]; y0d = [0; -2];
629 ng = 1;
630
631 // Simple simulation
632 t = 0:0.003:300;
633 yy = dae([y0, y0d], t0, t, atol, rtol, 'res22', 'jac22');
634 clf(); plot(yy(1, :), yy(2, :))
635 // Find first point where yy(1) = 0
636 %DAEOPTIONS = list([] , 0, [], [], [], 0, [], 1, [], 0, 1, [], [], 1);
637 [yy, nn, hotd] = dae("root2", [y0, y0d], t0, 300, atol, rtol, 'res22', 'jac22', ng, 'gr22', 'psol1', 'pjac1');
638 plot(yy(1, 1), yy(2, 1), 'r+')
639 xstring(yy(1, 1)+0.1, yy(2, 1), string(nn(1)));
640
641 // Hot restart for next point
642 t01 = nn(1);
643 [pp, qq] = size(yy);
644 y01 = yy(2:3, qq); y0d1 = yy(3:4, qq);
645 [yy, nn, hotd] = dae("root2", [y01, y0d1], t01, 300, atol, rtol, 'res22', 'jac22', ng, 'gr22', 'psol1', 'pjac1', hotd);
646 plot(yy(1, 1), yy(2, 1), 'r+')
647 xstring(yy(1, 1)+0.1, yy(2, 1), string(nn(1)));
648 cd(previous_dir);
649      ]]></programlisting>
650         <scilab:image><![CDATA[
651 code = ['#include <math.h>'
652       'void res22(double *t, double *y, double *yd, double *res, int *ires, double *rpar, int *ipar)'
653       '{res[0] = yd[0] - y[1];'
654       ' res[1] = yd[1] - (100.0*(1.0 - y[0]*y[0])*y[1] - y[0]);}'
655       ' '
656       'void jac22(double *t, double *y, double *yd, double *pd, double *cj, double *rpar, int *ipar)'
657       '{pd[0] = *cj - 0.0;'
658       ' pd[1] =     - (-200.0*y[0]*y[1] - 1.0);'
659       ' pd[2] =     - 1.0;'
660       ' pd[3] = *cj - (100.0*(1.0 - y[0]*y[0]));}'
661       ' '
662       'void gr22(int *neq, double *t, double *y, int *ng, double *groot, double *rpar, int *ipar)'
663       '{ groot[0] = y[0];}']
664 previous_dir = pwd();
665 cd TMPDIR;
666 mputl(code, 't22.c')
667 ilib_for_link(['res22' 'jac22' 'gr22'], 't22.c', [], 'c', [], 't22loader.sce');
668 exec('t22loader.sce')
669 rtol = [1.d-6; 1.d-6];
670 atol = [1.d-6; 1.d-4];
671 t0 = 0; t = [20:20:200];
672 y0 = [2; 0]; y0d = [0; -2];
673 ng = 1;
674 t = 0:0.003:300;
675 yy = dae([y0, y0d], t0, t, atol, rtol, 'res22', 'jac22');
676 clf(); plot(yy(1, :), yy(2, :))
677 %DAEOPTIONS = list([], 0, [], [], [], 0, [], 1, [], 0, 1, [], [], 1);
678 [yy, nn, hotd] = dae("root2", [y0, y0d], t0, 300, atol, rtol, 'res22', 'jac22', ng, 'gr22', 'psol1', 'pjac1');
679 plot(yy(1, 1), yy(2, 1), 'r+')
680 xstring(yy(1, 1)+0.1, yy(2, 1), string(nn(1)));
681 t01 = nn(1);
682 [pp, qq] = size(yy);
683 y01 = yy(2:3, qq);
684 y0d1 = yy(3:4, qq);
685 [yy, nn, hotd] = dae("root2", [y01, y0d1], t01, 300, atol, rtol, 'res22', 'jac22', ng, 'gr22', 'psol1', 'pjac1', hotd);
686 plot(yy(1, 1), yy(2, 1), 'r+')
687 xstring(yy(1, 1)+0.1, yy(2, 1), string(nn(1)));
688 cd(previous_dir);
689  ]]></scilab:image>
690     </refsection>
691     <refsection role="see also">
692         <title>参照</title>
693         <simplelist type="inline">
694             <member>
695                 <link linkend="ode">ode</link>
696             </member>
697             <member>
698                 <link linkend="daeoptions">daeoptions</link>
699             </member>
700             <member>
701                 <link linkend="dassl">dassl</link>
702             </member>
703             <member>
704                 <link linkend="dasrt">dasrt</link>
705             </member>
706             <member>
707                 <link linkend="daskr">daskr</link>
708             </member>
709             <member>
710                 <link linkend="impl">impl</link>
711             </member>
712             <member>
713                 <link linkend="call">call</link>
714             </member>
715             <member>
716                 <link linkend="link">link</link>
717             </member>
718             <member>
719                 <link linkend="external">external</link>
720             </member>
721         </simplelist>
722     </refsection>
723 </refentry>