Mercurial > octave-nkf

comparison src/dassl.cc @ 1:78fd87e624cb

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1993-08-08 01:13:40 by jwe] Initial revision
author jwe
date Sun, 08 Aug 1993 01:13:40 +0000
parents
children d68036bcad4c
comparison
equal deleted inserted replaced
0:22412e3a4641 1:78fd87e624cb
1 // tc-dassl.cc -*- C++ -*-
2 /*
3
4 Copyright (C) 1993 John W. Eaton
5
6 This file is part of Octave.
7
8 Octave is free software; you can redistribute it and/or modify it
9 under the terms of the GNU General Public License as published by the
10 Free Software Foundation; either version 2, or (at your option) any
11 later version.
12
13 Octave is distributed in the hope that it will be useful, but WITHOUT
14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
16 for more details.
17
18 You should have received a copy of the GNU General Public License
19 along with Octave; see the file COPYING. If not, write to the Free
20 Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
21
22 */
23
24 #ifdef __GNUG__
25 #pragma implementation
26 #endif
27
28 #include "DAE.h"
29
30 #include "tree-const.h"
31 #include "variables.h"
32 #include "gripes.h"
33 #include "error.h"
34 #include "utils.h"
35
36 // Global pointer for user defined function required by dassl.
37 static tree *dassl_fcn;
38
39 #ifdef WITH_DLD
40 tree_constant *
41 builtin_dassl_2 (tree_constant *args, int nargin, int nargout)
42 {
43 return dassl (args, nargin, nargout);
44 }
45 #endif
46
47 ColumnVector
48 dassl_user_function (const ColumnVector& x, const ColumnVector& xdot, double t)
49 {
50 ColumnVector retval;
51
52 int nstates = x.capacity ();
53
54 assert (nstates == xdot.capacity ());
55
56 // tree_constant name (dassl_fcn->name ());
57 tree_constant *args = new tree_constant [4];
58 // args[0] = name;
59 args[3] = tree_constant (t);
60
61 if (nstates > 1)
62 {
63 Matrix m1 (nstates, 1);
64 Matrix m2 (nstates, 1);
65 for (int i = 0; i < nstates; i++)
66 {
67 m1 (i, 0) = x.elem (i);
68 m2 (i, 0) = xdot.elem (i);
69 }
70 tree_constant state (m1);
71 tree_constant deriv (m2);
72 args[1] = state;
73 args[2] = deriv;
74 }
75 else
76 {
77 double d1 = x.elem (0);
78 double d2 = xdot.elem (0);
79 tree_constant state (d1);
80 tree_constant deriv (d2);
81 args[1] = state;
82 args[2] = deriv;
83 }
84
85 if (dassl_fcn != NULL_TREE)
86 {
87 tree_constant *tmp = dassl_fcn->eval (args, 4, 1, 0);
88 delete [] args;
89 if (tmp != NULL_TREE_CONST && tmp[0].is_defined ())
90 {
91 retval = tmp[0].to_vector ();
92 delete [] tmp;
93 }
94 else
95 {
96 delete [] tmp;
97 gripe_user_supplied_eval ("dassl");
98 jump_to_top_level ();
99 }
100 }
101
102 return retval;
103 }
104
105 tree_constant *
106 dassl (tree_constant *args, int nargin, int nargout)
107 {
108 // Assumes that we have been given the correct number of arguments.
109
110 tree_constant *retval = NULL_TREE_CONST;
111
112 dassl_fcn = is_valid_function (args[1], "dassl", 1);
113 if (dassl_fcn == NULL_TREE
114 || takes_correct_nargs (dassl_fcn, 4, "dassl", 1) != 1)
115 return retval;
116
117 ColumnVector state = args[2].to_vector ();
118 ColumnVector deriv = args[3].to_vector ();
119 ColumnVector out_times = args[4].to_vector ();
120 ColumnVector crit_times;
121 int crit_times_set = 0;
122 if (nargin > 5)
123 {
124 crit_times = args[5].to_vector ();
125 crit_times_set = 1;
126 }
127
128 if (state.capacity () != deriv.capacity ())
129 {
130 message ("dassl", "x and xdot must have the same size");
131 return retval;
132 }
133
134 double tzero = out_times.elem (0);
135
136 DAEFunc func (dassl_user_function);
137 DAE dae (state, deriv, tzero, func);
138
139 Matrix output;
140 Matrix deriv_output;
141
142 if (crit_times_set)
143 output = dae.integrate (out_times, deriv_output, crit_times);
144 else
145 output = dae.integrate (out_times, deriv_output);
146
147 retval = new tree_constant [3];
148 retval[0] = tree_constant (output);
149 retval[1] = tree_constant (deriv_output);
150 return retval;
151 }
152
153 /*
154 ;;; Local Variables: ***
155 ;;; mode: C++ ***
156 ;;; page-delimiter: "^/\\*" ***
157 ;;; End: ***
158 */
159