annotate src/DLD-FUNCTIONS/dasrt.cc @ 4140:303b28a7a7e4

[project @ 2002-11-01 02:53:13 by jwe]
author jwe
date Fri, 01 Nov 2002 02:53:14 +0000
parents 19a1626b8d57
children c0121a3b9cbe
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
3990
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
1 /*
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
2
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
3 Copyright (C) 2002 John W. Eaton
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
4
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
5 This file is part of Octave.
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
6
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
7 Octave is free software; you can redistribute it and/or modify it
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
8 under the terms of the GNU General Public License as published by the
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
9 Free Software Foundation; either version 2, or (at your option) any
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
10 later version.
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
11
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
12 Octave is distributed in the hope that it will be useful, but WITHOUT
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
15 for more details.
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
16
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
17 You should have received a copy of the GNU General Public License
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
18 along with Octave; see the file COPYING. If not, write to the Free
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
19 Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
20
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
21 */
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
22
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
23 #ifdef HAVE_CONFIG_H
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
24 #include <config.h>
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
25 #endif
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
26
4052
c2562b2593e2 [project @ 2002-08-17 22:15:13 by jwe]
jwe
parents: 4035
diff changeset
27 #include <iostream>
3990
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
28 #include <string>
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
29
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
30 #include "DASRT.h"
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
31 #include "lo-mappers.h"
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
32
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
33 #include "defun-dld.h"
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
34 #include "error.h"
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
35 #include "gripes.h"
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
36 #include "oct-obj.h"
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
37 #include "ov-fcn.h"
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
38 #include "pager.h"
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
39 #include "parse.h"
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
40 #include "unwind-prot.h"
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
41 #include "utils.h"
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
42 #include "variables.h"
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
43
3998
f6df65db67f9 [project @ 2002-07-24 18:10:39 by jwe]
jwe
parents: 3997
diff changeset
44 #include "DASRT-opts.cc"
f6df65db67f9 [project @ 2002-07-24 18:10:39 by jwe]
jwe
parents: 3997
diff changeset
45
4115
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
46 // Global pointers for user defined function required by dasrt.
3990
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
47 static octave_function *dasrt_f;
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
48 static octave_function *dasrt_j;
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
49 static octave_function *dasrt_cf;
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
50
4140
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
51 // Have we warned about imaginary values returned from user function?
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
52 static bool warned_fcn_imaginary = false;
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
53 static bool warned_jac_imaginary = false;
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
54 static bool warned_cf_imaginary = false;
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
55
3990
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
56 // Is this a recursive call?
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
57 static int call_depth = 0;
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
58
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
59 static ColumnVector
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
60 dasrt_user_f (const ColumnVector& x, const ColumnVector& xprime,
3993
f23bc69132cc [project @ 2002-07-16 20:18:56 by jwe]
jwe
parents: 3992
diff changeset
61 double t, int& ires)
3990
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
62 {
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
63 ColumnVector retval;
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
64
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
65 octave_value_list args;
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
66
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
67 int n = x.length ();
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
68
3993
f23bc69132cc [project @ 2002-07-16 20:18:56 by jwe]
jwe
parents: 3992
diff changeset
69 args(2) = t;
f23bc69132cc [project @ 2002-07-16 20:18:56 by jwe]
jwe
parents: 3992
diff changeset
70
3990
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
71 if (n > 1)
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
72 {
3993
f23bc69132cc [project @ 2002-07-16 20:18:56 by jwe]
jwe
parents: 3992
diff changeset
73 args(1) = xprime;
3990
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
74 args(0) = x;
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
75 }
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
76 else if (n == 1)
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
77 {
3993
f23bc69132cc [project @ 2002-07-16 20:18:56 by jwe]
jwe
parents: 3992
diff changeset
78 args(1) = xprime(0);
3990
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
79 args(0) = x(0);
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
80 }
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
81 else
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
82 {
3993
f23bc69132cc [project @ 2002-07-16 20:18:56 by jwe]
jwe
parents: 3992
diff changeset
83 args(1) = Matrix ();
3990
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
84 args(0) = Matrix ();
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
85 }
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
86
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
87 if (dasrt_f)
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
88 {
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
89 octave_value_list tmp = dasrt_f->do_multi_index_op (1, args);
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
90
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
91 if (error_state)
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
92 {
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
93 gripe_user_supplied_eval ("dasrt");
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
94 return retval;
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
95 }
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
96
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
97 if (tmp.length () > 0 && tmp(0).is_defined ())
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
98 {
4140
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
99 if (! warned_fcn_imaginary && tmp(0).is_complex_type ())
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
100 {
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
101 warning ("dasrt: ignoring imaginary part returned from user-supplied function");
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
102 warned_fcn_imaginary = true;
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
103 }
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
104
3990
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
105 retval = ColumnVector (tmp(0).vector_value ());
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
106
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
107 if (error_state || retval.length () == 0)
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
108 gripe_user_supplied_eval ("dasrt");
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
109 }
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
110 else
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
111 gripe_user_supplied_eval ("dasrt");
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
112 }
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
113
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
114 return retval;
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
115 }
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
116
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
117 static ColumnVector
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
118 dasrt_user_cf (const ColumnVector& x, double t)
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
119 {
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
120 ColumnVector retval;
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
121
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
122 octave_value_list args;
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
123
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
124 int n = x.length ();
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
125
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
126 if (n > 1)
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
127 args(0) = x;
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
128 else if (n == 1)
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
129 args(0) = x(0);
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
130 else
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
131 args(0) = Matrix ();
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
132
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
133 args(1) = t;
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
134
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
135 if (dasrt_cf)
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
136 {
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
137 octave_value_list tmp = dasrt_cf->do_multi_index_op (1, args);
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
138
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
139 if (error_state)
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
140 {
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
141 gripe_user_supplied_eval ("dasrt");
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
142 return retval;
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
143 }
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
144
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
145 if (tmp.length () > 0 && tmp(0).is_defined ())
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
146 {
4140
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
147 if (! warned_cf_imaginary && tmp(0).is_complex_type ())
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
148 {
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
149 warning ("dasrt: ignoring imaginary part returned from user-supplied constraint function");
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
150 warned_cf_imaginary = true;
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
151 }
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
152
3990
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
153 retval = ColumnVector (tmp(0).vector_value ());
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
154
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
155 if (error_state || retval.length () == 0)
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
156 gripe_user_supplied_eval ("dasrt");
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
157 }
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
158 else
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
159 gripe_user_supplied_eval ("dasrt");
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
160 }
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
161
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
162 return retval;
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
163 }
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
164
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
165 static Matrix
3993
f23bc69132cc [project @ 2002-07-16 20:18:56 by jwe]
jwe
parents: 3992
diff changeset
166 dasrt_user_j (const ColumnVector& x, const ColumnVector& xdot,
f23bc69132cc [project @ 2002-07-16 20:18:56 by jwe]
jwe
parents: 3992
diff changeset
167 double t, double cj)
3990
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
168 {
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
169 Matrix retval;
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
170
3993
f23bc69132cc [project @ 2002-07-16 20:18:56 by jwe]
jwe
parents: 3992
diff changeset
171 int nstates = x.capacity ();
f23bc69132cc [project @ 2002-07-16 20:18:56 by jwe]
jwe
parents: 3992
diff changeset
172
f23bc69132cc [project @ 2002-07-16 20:18:56 by jwe]
jwe
parents: 3992
diff changeset
173 assert (nstates == xdot.capacity ());
3990
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
174
3993
f23bc69132cc [project @ 2002-07-16 20:18:56 by jwe]
jwe
parents: 3992
diff changeset
175 octave_value_list args;
f23bc69132cc [project @ 2002-07-16 20:18:56 by jwe]
jwe
parents: 3992
diff changeset
176
f23bc69132cc [project @ 2002-07-16 20:18:56 by jwe]
jwe
parents: 3992
diff changeset
177 args(3) = cj;
f23bc69132cc [project @ 2002-07-16 20:18:56 by jwe]
jwe
parents: 3992
diff changeset
178 args(2) = t;
3990
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
179
3993
f23bc69132cc [project @ 2002-07-16 20:18:56 by jwe]
jwe
parents: 3992
diff changeset
180 if (nstates > 1)
f23bc69132cc [project @ 2002-07-16 20:18:56 by jwe]
jwe
parents: 3992
diff changeset
181 {
f23bc69132cc [project @ 2002-07-16 20:18:56 by jwe]
jwe
parents: 3992
diff changeset
182 Matrix m1 (nstates, 1);
f23bc69132cc [project @ 2002-07-16 20:18:56 by jwe]
jwe
parents: 3992
diff changeset
183 Matrix m2 (nstates, 1);
f23bc69132cc [project @ 2002-07-16 20:18:56 by jwe]
jwe
parents: 3992
diff changeset
184 for (int i = 0; i < nstates; i++)
f23bc69132cc [project @ 2002-07-16 20:18:56 by jwe]
jwe
parents: 3992
diff changeset
185 {
f23bc69132cc [project @ 2002-07-16 20:18:56 by jwe]
jwe
parents: 3992
diff changeset
186 m1 (i, 0) = x (i);
f23bc69132cc [project @ 2002-07-16 20:18:56 by jwe]
jwe
parents: 3992
diff changeset
187 m2 (i, 0) = xdot (i);
f23bc69132cc [project @ 2002-07-16 20:18:56 by jwe]
jwe
parents: 3992
diff changeset
188 }
f23bc69132cc [project @ 2002-07-16 20:18:56 by jwe]
jwe
parents: 3992
diff changeset
189 octave_value state (m1);
f23bc69132cc [project @ 2002-07-16 20:18:56 by jwe]
jwe
parents: 3992
diff changeset
190 octave_value deriv (m2);
f23bc69132cc [project @ 2002-07-16 20:18:56 by jwe]
jwe
parents: 3992
diff changeset
191 args(1) = deriv;
f23bc69132cc [project @ 2002-07-16 20:18:56 by jwe]
jwe
parents: 3992
diff changeset
192 args(0) = state;
f23bc69132cc [project @ 2002-07-16 20:18:56 by jwe]
jwe
parents: 3992
diff changeset
193 }
f23bc69132cc [project @ 2002-07-16 20:18:56 by jwe]
jwe
parents: 3992
diff changeset
194 else
f23bc69132cc [project @ 2002-07-16 20:18:56 by jwe]
jwe
parents: 3992
diff changeset
195 {
f23bc69132cc [project @ 2002-07-16 20:18:56 by jwe]
jwe
parents: 3992
diff changeset
196 double d1 = x (0);
f23bc69132cc [project @ 2002-07-16 20:18:56 by jwe]
jwe
parents: 3992
diff changeset
197 double d2 = xdot (0);
f23bc69132cc [project @ 2002-07-16 20:18:56 by jwe]
jwe
parents: 3992
diff changeset
198 octave_value state (d1);
f23bc69132cc [project @ 2002-07-16 20:18:56 by jwe]
jwe
parents: 3992
diff changeset
199 octave_value deriv (d2);
f23bc69132cc [project @ 2002-07-16 20:18:56 by jwe]
jwe
parents: 3992
diff changeset
200 args(1) = deriv;
f23bc69132cc [project @ 2002-07-16 20:18:56 by jwe]
jwe
parents: 3992
diff changeset
201 args(0) = state;
f23bc69132cc [project @ 2002-07-16 20:18:56 by jwe]
jwe
parents: 3992
diff changeset
202 }
3990
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
203
3993
f23bc69132cc [project @ 2002-07-16 20:18:56 by jwe]
jwe
parents: 3992
diff changeset
204 if (dasrt_j)
f23bc69132cc [project @ 2002-07-16 20:18:56 by jwe]
jwe
parents: 3992
diff changeset
205 {
f23bc69132cc [project @ 2002-07-16 20:18:56 by jwe]
jwe
parents: 3992
diff changeset
206 octave_value_list tmp = dasrt_j->do_multi_index_op (1, args);
3990
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
207
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
208 if (error_state)
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
209 {
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
210 gripe_user_supplied_eval ("dasrt");
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
211 return retval;
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
212 }
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
213
3993
f23bc69132cc [project @ 2002-07-16 20:18:56 by jwe]
jwe
parents: 3992
diff changeset
214 int tlen = tmp.length ();
f23bc69132cc [project @ 2002-07-16 20:18:56 by jwe]
jwe
parents: 3992
diff changeset
215 if (tlen > 0 && tmp(0).is_defined ())
3990
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
216 {
4140
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
217 if (! warned_jac_imaginary && tmp(0).is_complex_type ())
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
218 {
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
219 warning ("dasrt: ignoring imaginary part returned from user-supplied jacobian function");
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
220 warned_jac_imaginary = true;
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
221 }
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
222
3990
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
223 retval = tmp(0).matrix_value ();
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
224
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
225 if (error_state || retval.length () == 0)
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
226 gripe_user_supplied_eval ("dasrt");
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
227 }
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
228 else
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
229 gripe_user_supplied_eval ("dasrt");
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
230 }
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
231
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
232 return retval;
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
233 }
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
234
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
235 #define DASRT_ABORT \
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
236 do \
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
237 { \
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
238 unwind_protect::run_frame ("Fdasrt"); \
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
239 return retval; \
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
240 } \
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
241 while (0)
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
242
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
243 #define DASRT_ABORT1(msg) \
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
244 do \
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
245 { \
4035
92776b806c55 [project @ 2002-08-10 08:04:07 by jwe]
jwe
parents: 3998
diff changeset
246 ::error ("dasrt: " msg); \
3990
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
247 DASRT_ABORT; \
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
248 } \
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
249 while (0)
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
250
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
251 #define DASRT_ABORT2(fmt, arg) \
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
252 do \
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
253 { \
4035
92776b806c55 [project @ 2002-08-10 08:04:07 by jwe]
jwe
parents: 3998
diff changeset
254 ::error ("dasrt: " fmt, arg); \
3990
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
255 DASRT_ABORT; \
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
256 } \
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
257 while (0)
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
258
3997
d4091aff6468 [project @ 2002-07-17 18:00:06 by jwe]
jwe
parents: 3994
diff changeset
259 DEFUN_DLD (dasrt, args, nargout,
3990
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
260 "-*- texinfo -*-\n\
4115
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
261 @deftypefn {Loadable Function} {[@var{x}, @var{xdot}, @var{t_out}, @var{istat}, @var{msg}] =} dasrt (@var{fcn} [, @var{g}], @var{x_0}, @var{xdot_0}, @var{t} [, @var{t_crit}])\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
262 Solve the set of differential-algebraic equations\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
263 @tex\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
264 $$ 0 = f (\\dot{x}, x, t) $$\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
265 with\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
266 $$ x(t_0) = x_0, \\dot{x}(t_0) = \\dot{x}_0 $$\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
267 @end tex\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
268 @ifinfo\n\
3990
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
269 \n\
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
270 @example\n\
4115
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
271 0 = f (xdot, x, t)\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
272 @end example\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
273 \n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
274 with\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
275 \n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
276 @example\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
277 x(t_0) = x_0, xdot(t_0) = xdot_0\n\
3990
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
278 @end example\n\
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
279 \n\
4115
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
280 @end ifinfo\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
281 with functional stopping criteria (root solving).\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
282 \n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
283 The solution is returned in the matrices @var{x} and @var{xdot},\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
284 with each row in the result matrices corresponding to one of the\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
285 elements in the vector @var{t_out}. The first element of @var{t}\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
286 should be @math{t_0} and correspond to the initial state of the\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
287 system @var{x_0} and its derivative @var{xdot_0}, so that the first\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
288 row of the output @var{x} is @var{x_0} and the first row\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
289 of the output @var{xdot} is @var{xdot_0}.\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
290 \n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
291 The vector @var{t} provides an upper limit on the length of the\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
292 integration. If the stopping condition is met, the vector\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
293 @var{t_out} will be shorter than @var{t}, and the final element of\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
294 @var{t_out} will be the point at which the stopping condition was met,\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
295 and may not correspond to any element of the vector @var{t}.\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
296 \n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
297 The first argument, @var{fcn}, is a string that names the function to\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
298 call to compute the vector of residuals for the set of equations.\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
299 It must have the form\n\
3990
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
300 \n\
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
301 @example\n\
4115
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
302 @var{res} = f (@var{x}, @var{xdot}, @var{t})\n\
3990
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
303 @end example\n\
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
304 \n\
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
305 @noindent\n\
4115
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
306 in which @var{x}, @var{xdot}, and @var{res} are vectors, and @var{t} is a\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
307 scalar.\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
308 \n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
309 If @var{fcn} is a two-element string array, the first element names\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
310 the function @math{f} described above, and the second element names\n\
4117
944b276d8856 [project @ 2002-10-23 01:06:46 by jwe]
jwe
parents: 4115
diff changeset
311 a function to compute the modified Jacobian\n\
3990
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
312 \n\
4115
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
313 @tex\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
314 $$\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
315 J = {\\partial f \\over \\partial x}\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
316 + c {\\partial f \\over \\partial \\dot{x}}\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
317 $$\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
318 @end tex\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
319 @ifinfo\n\
3990
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
320 \n\
4115
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
321 @example\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
322 df df\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
323 jac = -- + c ------\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
324 dx d xdot\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
325 @end example\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
326 \n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
327 @end ifinfo\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
328 \n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
329 The modified Jacobian function must have the form\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
330 \n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
331 @example\n\
3990
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
332 \n\
4115
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
333 @var{jac} = j (@var{x}, @var{xdot}, @var{t}, @var{c})\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
334 \n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
335 @end example\n\
3990
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
336 \n\
4115
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
337 The optional second argument names a function that defines the\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
338 constraint functions whose roots are desired during the integration.\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
339 This function must have the form\n\
3990
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
340 \n\
4115
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
341 @example\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
342 @var{g_out} = g (@var{x}, @var{t})\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
343 @end example\n\
3990
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
344 \n\
4115
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
345 and return a vector of the constraint function values.\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
346 If the value of any of the constraint functions changes sign, @sc{Dasrt}\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
347 will attempt to stop the integration at the point of the sign change.\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
348 \n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
349 If the name of the constraint function is omitted, @code{dasrt} solves\n\
4117
944b276d8856 [project @ 2002-10-23 01:06:46 by jwe]
jwe
parents: 4115
diff changeset
350 the saem problem as @code{daspk} or @code{dassl}.\n\
3990
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
351 \n\
4115
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
352 Note that because of numerical errors in the constraint functions\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
353 due to roundoff and integration error, @sc{Dasrt} may return false\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
354 roots, or return the same root at two or more nearly equal values of\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
355 @var{T}. If such false roots are suspected, the user should consider\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
356 smaller error tolerances or higher precision in the evaluation of the\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
357 constraint functions.\n\
3990
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
358 \n\
4115
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
359 If a root of some constraint function defines the end of the problem,\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
360 the input to @sc{Dasrt} should nevertheless allow integration to a\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
361 point slightly past that root, so that @sc{Dasrt} can locate the root\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
362 by interpolation.\n\
3990
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
363 \n\
4115
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
364 The third and fourth arguments to @code{dasrt} specify the initial\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
365 condition of the states and their derivatives, and the fourth argument\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
366 specifies a vector of output times at which the solution is desired,\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
367 including the time corresponding to the initial condition.\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
368 \n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
369 The set of initial states and derivatives are not strictly required to\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
370 be consistent. In practice, however, @sc{Dassl} is not very good at\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
371 determining a consistent set for you, so it is best if you ensure that\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
372 the initial values result in the function evaluating to zero.\n\
3990
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
373 \n\
4115
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
374 The sixth argument is optional, and may be used to specify a set of\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
375 times that the DAE solver should not integrate past. It is useful for\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
376 avoiding difficulties with singularities and points where there is a\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
377 discontinuity in the derivative.\n\
3990
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
378 \n\
4115
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
379 After a successful computation, the value of @var{istate} will be\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
380 greater than zero (consistent with the Fortran version of @sc{Dassl}).\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
381 \n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
382 If the computation is not successful, the value of @var{istate} will be\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
383 less than zero and @var{msg} will contain additional information.\n\
3990
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
384 \n\
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
385 You can use the function @code{dasrt_options} to set optional\n\
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
386 parameters for @code{dasrt}.\n\
4115
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
387 @end deftypefn\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 4052
diff changeset
388 @seealso{daspk, dasrt, lsode, odessa}")
3990
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
389 {
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
390 octave_value_list retval;
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
391
4140
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
392 warned_fcn_imaginary = false;
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
393 warned_jac_imaginary = false;
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
394 warned_cf_imaginary = false;
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
395
3990
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
396 unwind_protect::begin_frame ("Fdasrt");
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
397
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
398 unwind_protect_int (call_depth);
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
399 call_depth++;
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
400
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
401 if (call_depth > 1)
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
402 DASRT_ABORT1 ("invalid recursive call");
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
403
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
404 int argp = 0;
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
405
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
406 int nargin = args.length ();
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
407
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
408 if (nargin < 4 || nargin > 6)
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
409 {
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
410 print_usage ("dasrt");
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
411 unwind_protect::run_frame ("Fdasrt");
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
412 return retval;
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
413 }
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
414
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
415 dasrt_f = 0;
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
416 dasrt_j = 0;
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
417 dasrt_cf = 0;
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
418
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
419 // Check all the arguments. Are they the right animals?
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
420
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
421 // Here's where I take care of f and j in one shot:
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
422
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
423 octave_value f_arg = args(0);
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
424
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
425 switch (f_arg.rows ())
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
426 {
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
427 case 1:
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
428 dasrt_f = extract_function
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
429 (args(0), "dasrt", "__dasrt_fcn__",
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
430 "function res = __dasrt_fcn__ (x, xdot, t) res = ",
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
431 "; endfunction");
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
432 break;
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
433
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
434 case 2:
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
435 {
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
436 string_vector tmp = args(0).all_strings ();
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
437
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
438 if (! error_state)
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
439 {
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
440 dasrt_f = extract_function
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
441 (tmp(0), "dasrt", "__dasrt_fcn__",
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
442 "function res = __dasrt_fcn__ (x, xdot, t) res = ",
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
443 "; endfunction");
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
444
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
445 if (dasrt_f)
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
446 {
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
447 dasrt_j = extract_function
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
448 (tmp(1), "dasrt", "__dasrt_jac__",
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
449 "function jac = __dasrt_jac__ (x, xdot, t, cj) jac = ",
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
450 "; endfunction");
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
451
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
452 if (! dasrt_j)
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
453 dasrt_f = 0;
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
454 }
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
455 }
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
456 }
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
457 break;
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
458
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
459 default:
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
460 DASRT_ABORT1
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
461 ("first arg should be a string or 2-element string array");
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
462 }
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
463
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
464 if (error_state || (! dasrt_f))
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
465 DASRT_ABORT;
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
466
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
467 DAERTFunc func (dasrt_user_f);
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
468
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
469 argp++;
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
470
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
471 if (args(1).is_string ())
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
472 {
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
473 dasrt_cf = is_valid_function (args(1), "dasrt", true);
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
474
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
475 if (! dasrt_cf)
3992
53b4eab68976 [project @ 2002-07-16 19:36:52 by jwe]
jwe
parents: 3990
diff changeset
476 DASRT_ABORT1 ("expecting function name as argument 2");
3990
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
477
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
478 argp++;
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
479
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
480 func.set_constraint_function (dasrt_user_cf);
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
481 }
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
482
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
483 ColumnVector state (args(argp++).vector_value ());
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
484
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
485 if (error_state)
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
486 DASRT_ABORT2 ("expecting state vector as argument %d", argp);
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
487
3992
53b4eab68976 [project @ 2002-07-16 19:36:52 by jwe]
jwe
parents: 3990
diff changeset
488 ColumnVector stateprime (args(argp++).vector_value ());
3990
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
489
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
490 if (error_state)
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
491 DASRT_ABORT2
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
492 ("expecting time derivative of state vector as argument %d", argp);
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
493
3994
a41827ec5677 [project @ 2002-07-16 23:57:09 by jwe]
jwe
parents: 3993
diff changeset
494 ColumnVector out_times (args(argp++).vector_value ());
3990
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
495
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
496 if (error_state)
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
497 DASRT_ABORT2
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
498 ("expecting output time vector as %s argument %d", argp);
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
499
3994
a41827ec5677 [project @ 2002-07-16 23:57:09 by jwe]
jwe
parents: 3993
diff changeset
500 double tzero = out_times (0);
3990
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
501
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
502 ColumnVector crit_times;
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
503
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
504 bool crit_times_set = false;
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
505
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
506 if (argp < nargin)
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
507 {
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
508 crit_times = ColumnVector (args(argp++).vector_value ());
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
509
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
510 if (error_state)
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
511 DASRT_ABORT2
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
512 ("expecting critical time vector as argument %d", argp);
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
513
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
514 crit_times_set = true;
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
515 }
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
516
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
517 if (dasrt_j)
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
518 func.set_jacobian_function (dasrt_user_j);
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
519
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
520 DASRT_result output;
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
521
3992
53b4eab68976 [project @ 2002-07-16 19:36:52 by jwe]
jwe
parents: 3990
diff changeset
522 DASRT dae = DASRT (state, stateprime, tzero, func);
3990
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
523
4122
19a1626b8d57 [project @ 2002-10-23 22:10:53 by jwe]
jwe
parents: 4117
diff changeset
524 dae.set_options (dasrt_opts);
3990
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
525
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
526 if (crit_times_set)
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
527 output = dae.integrate (out_times, crit_times);
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
528 else
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
529 output = dae.integrate (out_times);
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
530
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
531 if (! error_state)
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
532 {
3997
d4091aff6468 [project @ 2002-07-17 18:00:06 by jwe]
jwe
parents: 3994
diff changeset
533 std::string msg = dae.error_message ();
d4091aff6468 [project @ 2002-07-17 18:00:06 by jwe]
jwe
parents: 3994
diff changeset
534
d4091aff6468 [project @ 2002-07-17 18:00:06 by jwe]
jwe
parents: 3994
diff changeset
535 retval(4) = msg;
d4091aff6468 [project @ 2002-07-17 18:00:06 by jwe]
jwe
parents: 3994
diff changeset
536 retval(3) = static_cast<double> (dae.integration_state ());
d4091aff6468 [project @ 2002-07-17 18:00:06 by jwe]
jwe
parents: 3994
diff changeset
537
d4091aff6468 [project @ 2002-07-17 18:00:06 by jwe]
jwe
parents: 3994
diff changeset
538 if (dae.integration_ok ())
d4091aff6468 [project @ 2002-07-17 18:00:06 by jwe]
jwe
parents: 3994
diff changeset
539 {
d4091aff6468 [project @ 2002-07-17 18:00:06 by jwe]
jwe
parents: 3994
diff changeset
540 retval(2) = output.times ();
d4091aff6468 [project @ 2002-07-17 18:00:06 by jwe]
jwe
parents: 3994
diff changeset
541 retval(1) = output.deriv ();
d4091aff6468 [project @ 2002-07-17 18:00:06 by jwe]
jwe
parents: 3994
diff changeset
542 retval(0) = output.state ();
d4091aff6468 [project @ 2002-07-17 18:00:06 by jwe]
jwe
parents: 3994
diff changeset
543 }
d4091aff6468 [project @ 2002-07-17 18:00:06 by jwe]
jwe
parents: 3994
diff changeset
544 else
d4091aff6468 [project @ 2002-07-17 18:00:06 by jwe]
jwe
parents: 3994
diff changeset
545 {
d4091aff6468 [project @ 2002-07-17 18:00:06 by jwe]
jwe
parents: 3994
diff changeset
546 retval(2) = Matrix ();
d4091aff6468 [project @ 2002-07-17 18:00:06 by jwe]
jwe
parents: 3994
diff changeset
547 retval(1) = Matrix ();
d4091aff6468 [project @ 2002-07-17 18:00:06 by jwe]
jwe
parents: 3994
diff changeset
548 retval(0) = Matrix ();
d4091aff6468 [project @ 2002-07-17 18:00:06 by jwe]
jwe
parents: 3994
diff changeset
549
d4091aff6468 [project @ 2002-07-17 18:00:06 by jwe]
jwe
parents: 3994
diff changeset
550 if (nargout < 4)
d4091aff6468 [project @ 2002-07-17 18:00:06 by jwe]
jwe
parents: 3994
diff changeset
551 error ("dasrt: %s", msg.c_str ());
d4091aff6468 [project @ 2002-07-17 18:00:06 by jwe]
jwe
parents: 3994
diff changeset
552 }
3990
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
553 }
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
554
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
555 unwind_protect::run_frame ("Fdasrt");
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
556
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
557 return retval;
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
558 }
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
559
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
560 /*
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
561 ;;; Local Variables: ***
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
562 ;;; mode: C++ ***
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
563 ;;; End: ***
46388d6a4e44 [project @ 2002-07-16 06:20:39 by jwe]
jwe
parents:
diff changeset
564 */