annotate src/DLD-FUNCTIONS/dassl.cc @ 4628:c0121a3b9cbe

[project @ 2003-11-17 20:19:07 by jwe]
author jwe
date Mon, 17 Nov 2003 20:19:07 +0000
parents 303b28a7a7e4
children 737fd5c0ae06
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
1 /*
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
2
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
3 Copyright (C) 1996, 1997 John W. Eaton
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
4
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
5 This file is part of Octave.
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
6
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
7 Octave is free software; you can redistribute it and/or modify it
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
8 under the terms of the GNU General Public License as published by the
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
9 Free Software Foundation; either version 2, or (at your option) any
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
10 later version.
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
11
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
12 Octave is distributed in the hope that it will be useful, but WITHOUT
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
15 for more details.
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
16
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
17 You should have received a copy of the GNU General Public License
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
18 along with Octave; see the file COPYING. If not, write to the Free
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
19 Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
20
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
21 */
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
22
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
23 #ifdef HAVE_CONFIG_H
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
24 #include <config.h>
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
25 #endif
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
26
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
27 #include <string>
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
28
3567
dd8998699aa8 [project @ 2000-02-03 09:54:28 by jwe]
jwe
parents: 3544
diff changeset
29 #include <iomanip>
3523
b80bbb43a1a9 [project @ 2000-02-02 10:25:52 by jwe]
jwe
parents: 3419
diff changeset
30 #include <iostream>
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
31
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
32 #include "DASSL.h"
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
33
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
34 #include "defun-dld.h"
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
35 #include "error.h"
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
36 #include "gripes.h"
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
37 #include "oct-obj.h"
2968
e330cb788508 [project @ 1997-05-15 17:40:29 by jwe]
jwe
parents: 2928
diff changeset
38 #include "ov-fcn.h"
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
39 #include "pager.h"
3243
dd00769643ae [project @ 1999-05-28 04:19:00 by jwe]
jwe
parents: 3131
diff changeset
40 #include "unwind-prot.h"
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
41 #include "utils.h"
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
42 #include "variables.h"
295f037b4b3e [project @ 1997-05-05 05:32:33 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 "DASSL-opts.cc"
f6df65db67f9 [project @ 2002-07-24 18:10:39 by jwe]
jwe
parents: 3997
diff changeset
45
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
46 // Global pointer for user defined function required by dassl.
2968
e330cb788508 [project @ 1997-05-15 17:40:29 by jwe]
jwe
parents: 2928
diff changeset
47 static octave_function *dassl_fcn;
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
48
3991
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
49 // Global pointer for optional user defined jacobian function.
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
50 static octave_function *dassl_jac;
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
51
4140
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
52 // 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
53 static bool warned_fcn_imaginary = false;
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
54 static bool warned_jac_imaginary = false;
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
55
3243
dd00769643ae [project @ 1999-05-28 04:19:00 by jwe]
jwe
parents: 3131
diff changeset
56 // Is this a recursive call?
dd00769643ae [project @ 1999-05-28 04:19:00 by jwe]
jwe
parents: 3131
diff changeset
57 static int call_depth = 0;
dd00769643ae [project @ 1999-05-28 04:19:00 by jwe]
jwe
parents: 3131
diff changeset
58
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
59 ColumnVector
3849
5266e351a19c [project @ 2001-11-02 04:50:09 by jwe]
jwe
parents: 3747
diff changeset
60 dassl_user_function (const ColumnVector& x, const ColumnVector& xdot,
5266e351a19c [project @ 2001-11-02 04:50:09 by jwe]
jwe
parents: 3747
diff changeset
61 double t, int& ires)
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
62 {
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
63 ColumnVector retval;
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
64
4628
c0121a3b9cbe [project @ 2003-11-17 20:19:07 by jwe]
jwe
parents: 4140
diff changeset
65 assert (x.capacity () == xdot.capacity ());
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
66
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
67 octave_value_list args;
4628
c0121a3b9cbe [project @ 2003-11-17 20:19:07 by jwe]
jwe
parents: 4140
diff changeset
68
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
69 args(2) = t;
4628
c0121a3b9cbe [project @ 2003-11-17 20:19:07 by jwe]
jwe
parents: 4140
diff changeset
70 args(1) = xdot;
c0121a3b9cbe [project @ 2003-11-17 20:19:07 by jwe]
jwe
parents: 4140
diff changeset
71 args(0) = x;
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
72
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
73 if (dassl_fcn)
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
74 {
3544
71bd2d124119 [project @ 2000-02-02 21:02:31 by jwe]
jwe
parents: 3523
diff changeset
75 octave_value_list tmp = dassl_fcn->do_multi_index_op (1, args);
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
76
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
77 if (error_state)
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
78 {
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
79 gripe_user_supplied_eval ("dassl");
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
80 return retval;
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
81 }
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
82
3849
5266e351a19c [project @ 2001-11-02 04:50:09 by jwe]
jwe
parents: 3747
diff changeset
83 int tlen = tmp.length ();
5266e351a19c [project @ 2001-11-02 04:50:09 by jwe]
jwe
parents: 3747
diff changeset
84 if (tlen > 0 && tmp(0).is_defined ())
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
85 {
4140
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
86 if (! warned_fcn_imaginary && tmp(0).is_complex_type ())
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
87 {
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
88 warning ("dassl: ignoring imaginary part returned from user-supplied function");
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
89 warned_fcn_imaginary = true;
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
90 }
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
91
3419
e71b3d1dd327 [project @ 2000-01-12 05:23:34 by jwe]
jwe
parents: 3373
diff changeset
92 retval = ColumnVector (tmp(0).vector_value ());
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
93
3849
5266e351a19c [project @ 2001-11-02 04:50:09 by jwe]
jwe
parents: 3747
diff changeset
94 if (tlen > 1)
5266e351a19c [project @ 2001-11-02 04:50:09 by jwe]
jwe
parents: 3747
diff changeset
95 ires = tmp(1).int_value ();
5266e351a19c [project @ 2001-11-02 04:50:09 by jwe]
jwe
parents: 3747
diff changeset
96
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
97 if (error_state || retval.length () == 0)
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
98 gripe_user_supplied_eval ("dassl");
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
99 }
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
100 else
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
101 gripe_user_supplied_eval ("dassl");
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
102 }
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
103
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
104 return retval;
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
105 }
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
106
3991
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
107 Matrix
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
108 dassl_user_jacobian (const ColumnVector& x, const ColumnVector& xdot,
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
109 double t, double cj)
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
110 {
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
111 Matrix retval;
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
112
4628
c0121a3b9cbe [project @ 2003-11-17 20:19:07 by jwe]
jwe
parents: 4140
diff changeset
113 assert (x.capacity () == xdot.capacity ());
3991
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
114
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
115 octave_value_list args;
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
116
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
117 args(3) = cj;
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
118 args(2) = t;
4628
c0121a3b9cbe [project @ 2003-11-17 20:19:07 by jwe]
jwe
parents: 4140
diff changeset
119 args(1) = xdot;
c0121a3b9cbe [project @ 2003-11-17 20:19:07 by jwe]
jwe
parents: 4140
diff changeset
120 args(0) = x;
3991
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
121
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
122 if (dassl_jac)
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
123 {
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
124 octave_value_list tmp = dassl_jac->do_multi_index_op (1, args);
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
125
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
126 if (error_state)
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
127 {
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
128 gripe_user_supplied_eval ("dassl");
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
129 return retval;
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
130 }
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
131
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
132 int tlen = tmp.length ();
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
133 if (tlen > 0 && tmp(0).is_defined ())
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
134 {
4140
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
135 if (! warned_jac_imaginary && tmp(0).is_complex_type ())
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
136 {
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
137 warning ("dassl: ignoring imaginary part returned from user-supplied jacobian function");
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
138 warned_jac_imaginary = true;
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
139 }
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
140
3991
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
141 retval = tmp(0).matrix_value ();
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
142
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
143 if (error_state || retval.length () == 0)
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
144 gripe_user_supplied_eval ("dassl");
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
145 }
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
146 else
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
147 gripe_user_supplied_eval ("dassl");
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
148 }
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
149
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
150 return retval;
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
151 }
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
152
3323
511caaa5e98e [project @ 1999-10-27 02:21:46 by jwe]
jwe
parents: 3321
diff changeset
153 #define DASSL_ABORT() \
511caaa5e98e [project @ 1999-10-27 02:21:46 by jwe]
jwe
parents: 3321
diff changeset
154 do \
511caaa5e98e [project @ 1999-10-27 02:21:46 by jwe]
jwe
parents: 3321
diff changeset
155 { \
511caaa5e98e [project @ 1999-10-27 02:21:46 by jwe]
jwe
parents: 3321
diff changeset
156 unwind_protect::run_frame ("Fdassl"); \
511caaa5e98e [project @ 1999-10-27 02:21:46 by jwe]
jwe
parents: 3321
diff changeset
157 return retval; \
511caaa5e98e [project @ 1999-10-27 02:21:46 by jwe]
jwe
parents: 3321
diff changeset
158 } \
511caaa5e98e [project @ 1999-10-27 02:21:46 by jwe]
jwe
parents: 3321
diff changeset
159 while (0)
511caaa5e98e [project @ 1999-10-27 02:21:46 by jwe]
jwe
parents: 3321
diff changeset
160
511caaa5e98e [project @ 1999-10-27 02:21:46 by jwe]
jwe
parents: 3321
diff changeset
161 #define DASSL_ABORT1(msg) \
511caaa5e98e [project @ 1999-10-27 02:21:46 by jwe]
jwe
parents: 3321
diff changeset
162 do \
511caaa5e98e [project @ 1999-10-27 02:21:46 by jwe]
jwe
parents: 3321
diff changeset
163 { \
3747
50f30e40abca [project @ 2000-12-07 05:47:06 by jwe]
jwe
parents: 3568
diff changeset
164 ::error ("dassl: " msg); \
3323
511caaa5e98e [project @ 1999-10-27 02:21:46 by jwe]
jwe
parents: 3321
diff changeset
165 DASSL_ABORT (); \
511caaa5e98e [project @ 1999-10-27 02:21:46 by jwe]
jwe
parents: 3321
diff changeset
166 } \
511caaa5e98e [project @ 1999-10-27 02:21:46 by jwe]
jwe
parents: 3321
diff changeset
167 while (0)
511caaa5e98e [project @ 1999-10-27 02:21:46 by jwe]
jwe
parents: 3321
diff changeset
168
511caaa5e98e [project @ 1999-10-27 02:21:46 by jwe]
jwe
parents: 3321
diff changeset
169 #define DASSL_ABORT2(fmt, arg) \
511caaa5e98e [project @ 1999-10-27 02:21:46 by jwe]
jwe
parents: 3321
diff changeset
170 do \
511caaa5e98e [project @ 1999-10-27 02:21:46 by jwe]
jwe
parents: 3321
diff changeset
171 { \
3747
50f30e40abca [project @ 2000-12-07 05:47:06 by jwe]
jwe
parents: 3568
diff changeset
172 ::error ("dassl: " fmt, arg); \
3323
511caaa5e98e [project @ 1999-10-27 02:21:46 by jwe]
jwe
parents: 3321
diff changeset
173 DASSL_ABORT (); \
511caaa5e98e [project @ 1999-10-27 02:21:46 by jwe]
jwe
parents: 3321
diff changeset
174 } \
511caaa5e98e [project @ 1999-10-27 02:21:46 by jwe]
jwe
parents: 3321
diff changeset
175 while (0)
511caaa5e98e [project @ 1999-10-27 02:21:46 by jwe]
jwe
parents: 3321
diff changeset
176
3991
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
177 DEFUN_DLD (dassl, args, nargout,
3373
36405da8e173 [project @ 1999-11-23 20:54:17 by jwe]
jwe
parents: 3325
diff changeset
178 "-*- texinfo -*-\n\
4115
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 3998
diff changeset
179 @deftypefn {Loadable Function} {[@var{x}, @var{xdot}, @var{istate}, @var{msg}] =} dassl (@var{fcn}, @var{x_0}, @var{xdot_0}, @var{t}, @var{t_crit})\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 3998
diff changeset
180 Solve the set of differential-algebraic equations\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 3998
diff changeset
181 @tex\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 3998
diff changeset
182 $$ 0 = f (\\dot{x}, x, t) $$\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 3998
diff changeset
183 with\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 3998
diff changeset
184 $$ 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: 3998
diff changeset
185 @end tex\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 3998
diff changeset
186 @ifinfo\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 3998
diff changeset
187 \n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 3998
diff changeset
188 @example\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 3998
diff changeset
189 0 = f (xdot, x, t)\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 3998
diff changeset
190 @end example\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 3998
diff changeset
191 \n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 3998
diff changeset
192 with\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 3998
diff changeset
193 \n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 3998
diff changeset
194 @example\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 3998
diff changeset
195 x(t_0) = x_0, xdot(t_0) = xdot_0\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 3998
diff changeset
196 @end example\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 3998
diff changeset
197 \n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 3998
diff changeset
198 @end ifinfo\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 3998
diff changeset
199 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: 3998
diff changeset
200 with each row in the result matrices corresponding to one of the\n\
3373
36405da8e173 [project @ 1999-11-23 20:54:17 by jwe]
jwe
parents: 3325
diff changeset
201 elements in the vector @var{t}. The first element of @var{t}\n\
4115
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 3998
diff changeset
202 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: 3998
diff changeset
203 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: 3998
diff changeset
204 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: 3998
diff changeset
205 of the output @var{xdot} is @var{xdot_0}.\n\
3373
36405da8e173 [project @ 1999-11-23 20:54:17 by jwe]
jwe
parents: 3325
diff changeset
206 \n\
36405da8e173 [project @ 1999-11-23 20:54:17 by jwe]
jwe
parents: 3325
diff changeset
207 The first argument, @var{fcn}, is a string that names the function to\n\
36405da8e173 [project @ 1999-11-23 20:54:17 by jwe]
jwe
parents: 3325
diff changeset
208 call to compute the vector of residuals for the set of equations.\n\
36405da8e173 [project @ 1999-11-23 20:54:17 by jwe]
jwe
parents: 3325
diff changeset
209 It must have the form\n\
36405da8e173 [project @ 1999-11-23 20:54:17 by jwe]
jwe
parents: 3325
diff changeset
210 \n\
36405da8e173 [project @ 1999-11-23 20:54:17 by jwe]
jwe
parents: 3325
diff changeset
211 @example\n\
36405da8e173 [project @ 1999-11-23 20:54:17 by jwe]
jwe
parents: 3325
diff changeset
212 @var{res} = f (@var{x}, @var{xdot}, @var{t})\n\
36405da8e173 [project @ 1999-11-23 20:54:17 by jwe]
jwe
parents: 3325
diff changeset
213 @end example\n\
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
214 \n\
3373
36405da8e173 [project @ 1999-11-23 20:54:17 by jwe]
jwe
parents: 3325
diff changeset
215 @noindent\n\
4115
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 3998
diff changeset
216 in which @var{x}, @var{xdot}, and @var{res} are vectors, and @var{t} is a\n\
3373
36405da8e173 [project @ 1999-11-23 20:54:17 by jwe]
jwe
parents: 3325
diff changeset
217 scalar.\n\
36405da8e173 [project @ 1999-11-23 20:54:17 by jwe]
jwe
parents: 3325
diff changeset
218 \n\
4115
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 3998
diff changeset
219 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: 3998
diff changeset
220 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
221 a function to compute the modified Jacobian\n\
4115
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 3998
diff changeset
222 \n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 3998
diff changeset
223 @tex\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 3998
diff changeset
224 $$\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 3998
diff changeset
225 J = {\\partial f \\over \\partial x}\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 3998
diff changeset
226 + c {\\partial f \\over \\partial \\dot{x}}\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 3998
diff changeset
227 $$\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 3998
diff changeset
228 @end tex\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 3998
diff changeset
229 @ifinfo\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 3998
diff changeset
230 df df\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 3998
diff changeset
231 jac = -- + c ------\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 3998
diff changeset
232 dx d xdot\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 3998
diff changeset
233 @example\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 3998
diff changeset
234 @end example\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 3998
diff changeset
235 \n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 3998
diff changeset
236 @end ifinfo\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 3998
diff changeset
237 \n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 3998
diff changeset
238 The modified Jacobian function must have the form\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 3998
diff changeset
239 \n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 3998
diff changeset
240 @example\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 3998
diff changeset
241 \n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 3998
diff changeset
242 @var{jac} = j (@var{x}, @var{xdot}, @var{t}, @var{c})\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 3998
diff changeset
243 \n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 3998
diff changeset
244 @end example\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 3998
diff changeset
245 \n\
3373
36405da8e173 [project @ 1999-11-23 20:54:17 by jwe]
jwe
parents: 3325
diff changeset
246 The second and third arguments to @code{dassl} specify the initial\n\
36405da8e173 [project @ 1999-11-23 20:54:17 by jwe]
jwe
parents: 3325
diff changeset
247 condition of the states and their derivatives, and the fourth argument\n\
4115
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 3998
diff changeset
248 specifies a vector of output times at which the solution is desired,\n\
3373
36405da8e173 [project @ 1999-11-23 20:54:17 by jwe]
jwe
parents: 3325
diff changeset
249 including the time corresponding to the initial condition.\n\
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
250 \n\
3373
36405da8e173 [project @ 1999-11-23 20:54:17 by jwe]
jwe
parents: 3325
diff changeset
251 The set of initial states and derivatives are not strictly required to\n\
36405da8e173 [project @ 1999-11-23 20:54:17 by jwe]
jwe
parents: 3325
diff changeset
252 be consistent. In practice, however, @sc{Dassl} is not very good at\n\
36405da8e173 [project @ 1999-11-23 20:54:17 by jwe]
jwe
parents: 3325
diff changeset
253 determining a consistent set for you, so it is best if you ensure that\n\
36405da8e173 [project @ 1999-11-23 20:54:17 by jwe]
jwe
parents: 3325
diff changeset
254 the initial values result in the function evaluating to zero.\n\
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
255 \n\
3373
36405da8e173 [project @ 1999-11-23 20:54:17 by jwe]
jwe
parents: 3325
diff changeset
256 The fifth argument is optional, and may be used to specify a set of\n\
36405da8e173 [project @ 1999-11-23 20:54:17 by jwe]
jwe
parents: 3325
diff changeset
257 times that the DAE solver should not integrate past. It is useful for\n\
36405da8e173 [project @ 1999-11-23 20:54:17 by jwe]
jwe
parents: 3325
diff changeset
258 avoiding difficulties with singularities and points where there is a\n\
36405da8e173 [project @ 1999-11-23 20:54:17 by jwe]
jwe
parents: 3325
diff changeset
259 discontinuity in the derivative.\n\
3964
3030cb6cb559 [project @ 2002-06-28 17:17:55 by jwe]
jwe
parents: 3849
diff changeset
260 \n\
4115
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 3998
diff changeset
261 After a successful computation, the value of @var{istate} will be\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 3998
diff changeset
262 greater than zero (consistent with the Fortran version of @sc{Dassl}).\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 3998
diff changeset
263 \n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 3998
diff changeset
264 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: 3998
diff changeset
265 less than zero and @var{msg} will contain additional information.\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 3998
diff changeset
266 \n\
3964
3030cb6cb559 [project @ 2002-06-28 17:17:55 by jwe]
jwe
parents: 3849
diff changeset
267 You can use the function @code{dassl_options} to set optional\n\
3030cb6cb559 [project @ 2002-06-28 17:17:55 by jwe]
jwe
parents: 3849
diff changeset
268 parameters for @code{dassl}.\n\
4115
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 3998
diff changeset
269 @end deftypefn\n\
fc2048d4cd21 [project @ 2002-10-22 21:28:42 by jwe]
jwe
parents: 3998
diff changeset
270 @seealso{daspk, dasrt, lsode, odessa}")
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
271 {
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
272 octave_value_list retval;
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
273
4140
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
274 warned_fcn_imaginary = false;
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
275 warned_jac_imaginary = false;
303b28a7a7e4 [project @ 2002-11-01 02:53:13 by jwe]
jwe
parents: 4122
diff changeset
276
3243
dd00769643ae [project @ 1999-05-28 04:19:00 by jwe]
jwe
parents: 3131
diff changeset
277 unwind_protect::begin_frame ("Fdassl");
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
278
3243
dd00769643ae [project @ 1999-05-28 04:19:00 by jwe]
jwe
parents: 3131
diff changeset
279 unwind_protect_int (call_depth);
dd00769643ae [project @ 1999-05-28 04:19:00 by jwe]
jwe
parents: 3131
diff changeset
280 call_depth++;
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
281
3243
dd00769643ae [project @ 1999-05-28 04:19:00 by jwe]
jwe
parents: 3131
diff changeset
282 if (call_depth > 1)
3323
511caaa5e98e [project @ 1999-10-27 02:21:46 by jwe]
jwe
parents: 3321
diff changeset
283 DASSL_ABORT1 ("invalid recursive call");
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
284
3243
dd00769643ae [project @ 1999-05-28 04:19:00 by jwe]
jwe
parents: 3131
diff changeset
285 int nargin = args.length ();
dd00769643ae [project @ 1999-05-28 04:19:00 by jwe]
jwe
parents: 3131
diff changeset
286
3991
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
287 if (nargin > 3 && nargin < 6 && nargout < 5)
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
288 {
3991
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
289 dassl_fcn = 0;
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
290 dassl_jac = 0;
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
291
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
292 octave_value f_arg = args(0);
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
293
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
294 switch (f_arg.rows ())
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
295 {
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
296 case 1:
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
297 dassl_fcn = extract_function
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
298 (f_arg, "dassl", "__dassl_fcn__",
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
299 "function res = __dassl_fcn__ (x, xdot, t) res = ",
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
300 "; endfunction");
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
301 break;
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
302
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
303 case 2:
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
304 {
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
305 string_vector tmp = f_arg.all_strings ();
3243
dd00769643ae [project @ 1999-05-28 04:19:00 by jwe]
jwe
parents: 3131
diff changeset
306
3991
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
307 if (! error_state)
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
308 {
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
309 dassl_fcn = extract_function
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
310 (tmp(0), "dassl", "__dassl_fcn__",
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
311 "function res = __dassl_fcn__ (x, xdot, t) res = ",
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
312 "; endfunction");
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
313
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
314 if (dassl_fcn)
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
315 {
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
316 dassl_jac = extract_function
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
317 (tmp(1), "dassl", "__dassl_jac__",
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
318 "function jac = __dassl_jac__ (x, xdot, t, cj) jac = ",
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
319 "; endfunction");
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
320
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
321 if (! dassl_jac)
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
322 dassl_fcn = 0;
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
323 }
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
324 }
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
325 }
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
326 }
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
327
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
328 if (error_state || ! dassl_fcn)
3323
511caaa5e98e [project @ 1999-10-27 02:21:46 by jwe]
jwe
parents: 3321
diff changeset
329 DASSL_ABORT ();
3243
dd00769643ae [project @ 1999-05-28 04:19:00 by jwe]
jwe
parents: 3131
diff changeset
330
3419
e71b3d1dd327 [project @ 2000-01-12 05:23:34 by jwe]
jwe
parents: 3373
diff changeset
331 ColumnVector state = ColumnVector (args(1).vector_value ());
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
332
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
333 if (error_state)
3323
511caaa5e98e [project @ 1999-10-27 02:21:46 by jwe]
jwe
parents: 3321
diff changeset
334 DASSL_ABORT1 ("expecting state vector as second argument");
3243
dd00769643ae [project @ 1999-05-28 04:19:00 by jwe]
jwe
parents: 3131
diff changeset
335
3419
e71b3d1dd327 [project @ 2000-01-12 05:23:34 by jwe]
jwe
parents: 3373
diff changeset
336 ColumnVector deriv (args(2).vector_value ());
3243
dd00769643ae [project @ 1999-05-28 04:19:00 by jwe]
jwe
parents: 3131
diff changeset
337
dd00769643ae [project @ 1999-05-28 04:19:00 by jwe]
jwe
parents: 3131
diff changeset
338 if (error_state)
3323
511caaa5e98e [project @ 1999-10-27 02:21:46 by jwe]
jwe
parents: 3321
diff changeset
339 DASSL_ABORT1 ("expecting derivative vector as third argument");
3243
dd00769643ae [project @ 1999-05-28 04:19:00 by jwe]
jwe
parents: 3131
diff changeset
340
3419
e71b3d1dd327 [project @ 2000-01-12 05:23:34 by jwe]
jwe
parents: 3373
diff changeset
341 ColumnVector out_times (args(3).vector_value ());
3243
dd00769643ae [project @ 1999-05-28 04:19:00 by jwe]
jwe
parents: 3131
diff changeset
342
dd00769643ae [project @ 1999-05-28 04:19:00 by jwe]
jwe
parents: 3131
diff changeset
343 if (error_state)
3323
511caaa5e98e [project @ 1999-10-27 02:21:46 by jwe]
jwe
parents: 3321
diff changeset
344 DASSL_ABORT1 ("expecting output time vector as fourth argument");
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
345
3243
dd00769643ae [project @ 1999-05-28 04:19:00 by jwe]
jwe
parents: 3131
diff changeset
346 ColumnVector crit_times;
dd00769643ae [project @ 1999-05-28 04:19:00 by jwe]
jwe
parents: 3131
diff changeset
347 int crit_times_set = 0;
dd00769643ae [project @ 1999-05-28 04:19:00 by jwe]
jwe
parents: 3131
diff changeset
348 if (nargin > 4)
dd00769643ae [project @ 1999-05-28 04:19:00 by jwe]
jwe
parents: 3131
diff changeset
349 {
3419
e71b3d1dd327 [project @ 2000-01-12 05:23:34 by jwe]
jwe
parents: 3373
diff changeset
350 crit_times = ColumnVector (args(4).vector_value ());
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
351
3243
dd00769643ae [project @ 1999-05-28 04:19:00 by jwe]
jwe
parents: 3131
diff changeset
352 if (error_state)
3323
511caaa5e98e [project @ 1999-10-27 02:21:46 by jwe]
jwe
parents: 3321
diff changeset
353 DASSL_ABORT1 ("expecting critical time vector as fifth argument");
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
354
3243
dd00769643ae [project @ 1999-05-28 04:19:00 by jwe]
jwe
parents: 3131
diff changeset
355 crit_times_set = 1;
dd00769643ae [project @ 1999-05-28 04:19:00 by jwe]
jwe
parents: 3131
diff changeset
356 }
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
357
3243
dd00769643ae [project @ 1999-05-28 04:19:00 by jwe]
jwe
parents: 3131
diff changeset
358 if (state.capacity () != deriv.capacity ())
3323
511caaa5e98e [project @ 1999-10-27 02:21:46 by jwe]
jwe
parents: 3321
diff changeset
359 DASSL_ABORT1 ("x and xdot must have the same size");
3243
dd00769643ae [project @ 1999-05-28 04:19:00 by jwe]
jwe
parents: 3131
diff changeset
360
dd00769643ae [project @ 1999-05-28 04:19:00 by jwe]
jwe
parents: 3131
diff changeset
361 double tzero = out_times (0);
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
362
3243
dd00769643ae [project @ 1999-05-28 04:19:00 by jwe]
jwe
parents: 3131
diff changeset
363 DAEFunc func (dassl_user_function);
3991
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
364 if (dassl_jac)
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
365 func.set_jacobian_function (dassl_user_jacobian);
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
366
3243
dd00769643ae [project @ 1999-05-28 04:19:00 by jwe]
jwe
parents: 3131
diff changeset
367 DASSL dae (state, deriv, tzero, func);
3991
48d2bc4a3729 [project @ 2002-07-16 17:46:50 by jwe]
jwe
parents: 3964
diff changeset
368
4122
19a1626b8d57 [project @ 2002-10-23 22:10:53 by jwe]
jwe
parents: 4117
diff changeset
369 dae.set_options (dassl_opts);
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
370
3243
dd00769643ae [project @ 1999-05-28 04:19:00 by jwe]
jwe
parents: 3131
diff changeset
371 Matrix output;
dd00769643ae [project @ 1999-05-28 04:19:00 by jwe]
jwe
parents: 3131
diff changeset
372 Matrix deriv_output;
dd00769643ae [project @ 1999-05-28 04:19:00 by jwe]
jwe
parents: 3131
diff changeset
373
dd00769643ae [project @ 1999-05-28 04:19:00 by jwe]
jwe
parents: 3131
diff changeset
374 if (crit_times_set)
dd00769643ae [project @ 1999-05-28 04:19:00 by jwe]
jwe
parents: 3131
diff changeset
375 output = dae.integrate (out_times, deriv_output, crit_times);
dd00769643ae [project @ 1999-05-28 04:19:00 by jwe]
jwe
parents: 3131
diff changeset
376 else
dd00769643ae [project @ 1999-05-28 04:19:00 by jwe]
jwe
parents: 3131
diff changeset
377 output = dae.integrate (out_times, deriv_output);
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
378
3243
dd00769643ae [project @ 1999-05-28 04:19:00 by jwe]
jwe
parents: 3131
diff changeset
379 if (! error_state)
dd00769643ae [project @ 1999-05-28 04:19:00 by jwe]
jwe
parents: 3131
diff changeset
380 {
3997
d4091aff6468 [project @ 2002-07-17 18:00:06 by jwe]
jwe
parents: 3991
diff changeset
381 std::string msg = dae.error_message ();
d4091aff6468 [project @ 2002-07-17 18:00:06 by jwe]
jwe
parents: 3991
diff changeset
382
d4091aff6468 [project @ 2002-07-17 18:00:06 by jwe]
jwe
parents: 3991
diff changeset
383 retval(3) = msg;
d4091aff6468 [project @ 2002-07-17 18:00:06 by jwe]
jwe
parents: 3991
diff changeset
384 retval(2) = static_cast<double> (dae.integration_state ());
d4091aff6468 [project @ 2002-07-17 18:00:06 by jwe]
jwe
parents: 3991
diff changeset
385
d4091aff6468 [project @ 2002-07-17 18:00:06 by jwe]
jwe
parents: 3991
diff changeset
386 if (dae.integration_ok ())
d4091aff6468 [project @ 2002-07-17 18:00:06 by jwe]
jwe
parents: 3991
diff changeset
387 {
d4091aff6468 [project @ 2002-07-17 18:00:06 by jwe]
jwe
parents: 3991
diff changeset
388 retval(1) = deriv_output;
d4091aff6468 [project @ 2002-07-17 18:00:06 by jwe]
jwe
parents: 3991
diff changeset
389 retval(0) = output;
d4091aff6468 [project @ 2002-07-17 18:00:06 by jwe]
jwe
parents: 3991
diff changeset
390 }
d4091aff6468 [project @ 2002-07-17 18:00:06 by jwe]
jwe
parents: 3991
diff changeset
391 else
d4091aff6468 [project @ 2002-07-17 18:00:06 by jwe]
jwe
parents: 3991
diff changeset
392 {
d4091aff6468 [project @ 2002-07-17 18:00:06 by jwe]
jwe
parents: 3991
diff changeset
393 retval(1) = Matrix ();
d4091aff6468 [project @ 2002-07-17 18:00:06 by jwe]
jwe
parents: 3991
diff changeset
394 retval(0) = Matrix ();
d4091aff6468 [project @ 2002-07-17 18:00:06 by jwe]
jwe
parents: 3991
diff changeset
395
d4091aff6468 [project @ 2002-07-17 18:00:06 by jwe]
jwe
parents: 3991
diff changeset
396 if (nargout < 3)
d4091aff6468 [project @ 2002-07-17 18:00:06 by jwe]
jwe
parents: 3991
diff changeset
397 error ("dassl: %s", msg.c_str ());
d4091aff6468 [project @ 2002-07-17 18:00:06 by jwe]
jwe
parents: 3991
diff changeset
398 }
3243
dd00769643ae [project @ 1999-05-28 04:19:00 by jwe]
jwe
parents: 3131
diff changeset
399 }
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
400 }
3243
dd00769643ae [project @ 1999-05-28 04:19:00 by jwe]
jwe
parents: 3131
diff changeset
401 else
dd00769643ae [project @ 1999-05-28 04:19:00 by jwe]
jwe
parents: 3131
diff changeset
402 print_usage ("dassl");
dd00769643ae [project @ 1999-05-28 04:19:00 by jwe]
jwe
parents: 3131
diff changeset
403
dd00769643ae [project @ 1999-05-28 04:19:00 by jwe]
jwe
parents: 3131
diff changeset
404 unwind_protect::run_frame ("Fdassl");
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
405
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
406 return retval;
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
407 }
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
408
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
409 /*
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
410 ;;; Local Variables: ***
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
411 ;;; mode: C++ ***
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
412 ;;; End: ***
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
413 */