Mercurial > octave-nkf
annotate src/DLD-FUNCTIONS/daspk.cc @ 9067:8970b4b10e9f
Cleanup documentation for quad.texi and diffeq.texi
Grammarcheck input .txi files
Spellcheck .texi files
| author | Rik <rdrider0-list@yahoo.com> |
|---|---|
| date | Sat, 28 Mar 2009 22:18:51 -0700 |
| parents | 7c02ec148a3c |
| children | 610bf90fce2a |
| rev | line source |
|---|---|
| 3912 | 1 /* |
| 2 | |
| 8920 | 3 Copyright (C) 1996, 1997, 2002, 2003, 2004, 2005, 2006, 2007, 2009 |
| 7017 | 4 John W. Eaton |
| 3912 | 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 | |
| 7016 | 10 Free Software Foundation; either version 3 of the License, or (at your |
| 11 option) any later version. | |
| 3912 | 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 | |
| 7016 | 19 along with Octave; see the file COPYING. If not, see |
| 20 <http://www.gnu.org/licenses/>. | |
| 3912 | 21 |
| 22 */ | |
| 23 | |
| 24 #ifdef HAVE_CONFIG_H | |
| 25 #include <config.h> | |
| 26 #endif | |
| 27 | |
| 28 #include <string> | |
| 29 | |
| 30 #include <iomanip> | |
| 31 #include <iostream> | |
| 32 | |
| 33 #include "DASPK.h" | |
| 34 | |
| 35 #include "defun-dld.h" | |
| 36 #include "error.h" | |
| 37 #include "gripes.h" | |
| 38 #include "oct-obj.h" | |
| 39 #include "ov-fcn.h" | |
| 5729 | 40 #include "ov-cell.h" |
| 3912 | 41 #include "pager.h" |
| 42 #include "unwind-prot.h" | |
| 43 #include "utils.h" | |
| 44 #include "variables.h" | |
| 45 | |
| 3998 | 46 #include "DASPK-opts.cc" |
| 47 | |
| 3912 | 48 // Global pointer for user defined function required by daspk. |
| 49 static octave_function *daspk_fcn; | |
| 50 | |
| 4140 | 51 // Global pointer for optional user defined jacobian function. |
| 52 static octave_function *daspk_jac; | |
| 53 | |
| 54 // Have we warned about imaginary values returned from user function? | |
| 55 static bool warned_fcn_imaginary = false; | |
| 56 static bool warned_jac_imaginary = false; | |
| 57 | |
| 3912 | 58 // Is this a recursive call? |
| 59 static int call_depth = 0; | |
| 60 | |
| 61 ColumnVector | |
| 62 daspk_user_function (const ColumnVector& x, const ColumnVector& xdot, | |
| 5275 | 63 double t, octave_idx_type& ires) |
| 3912 | 64 { |
| 65 ColumnVector retval; | |
| 66 | |
| 4628 | 67 assert (x.capacity () == xdot.capacity ()); |
| 3912 | 68 |
| 69 octave_value_list args; | |
| 4628 | 70 |
| 3912 | 71 args(2) = t; |
| 4628 | 72 args(1) = xdot; |
| 73 args(0) = x; | |
| 3912 | 74 |
| 75 if (daspk_fcn) | |
| 76 { | |
| 77 octave_value_list tmp = daspk_fcn->do_multi_index_op (1, args); | |
| 78 | |
| 79 if (error_state) | |
| 80 { | |
| 81 gripe_user_supplied_eval ("daspk"); | |
| 82 return retval; | |
| 83 } | |
| 84 | |
| 85 int tlen = tmp.length (); | |
| 86 if (tlen > 0 && tmp(0).is_defined ()) | |
| 87 { | |
| 4140 | 88 if (! warned_fcn_imaginary && tmp(0).is_complex_type ()) |
| 89 { | |
| 90 warning ("daspk: ignoring imaginary part returned from user-supplied function"); | |
| 91 warned_fcn_imaginary = true; | |
| 92 } | |
| 93 | |
| 3912 | 94 retval = ColumnVector (tmp(0).vector_value ()); |
| 95 | |
| 96 if (tlen > 1) | |
| 97 ires = tmp(1).int_value (); | |
| 98 | |
| 99 if (error_state || retval.length () == 0) | |
| 100 gripe_user_supplied_eval ("daspk"); | |
| 101 } | |
| 102 else | |
| 103 gripe_user_supplied_eval ("daspk"); | |
| 104 } | |
| 105 | |
| 106 return retval; | |
| 107 } | |
| 108 | |
| 4140 | 109 Matrix |
| 110 daspk_user_jacobian (const ColumnVector& x, const ColumnVector& xdot, | |
| 111 double t, double cj) | |
| 112 { | |
| 113 Matrix retval; | |
| 114 | |
| 4628 | 115 assert (x.capacity () == xdot.capacity ()); |
| 4140 | 116 |
| 117 octave_value_list args; | |
| 118 | |
| 119 args(3) = cj; | |
| 120 args(2) = t; | |
| 4628 | 121 args(1) = xdot; |
| 122 args(0) = x; | |
| 4140 | 123 |
| 124 if (daspk_jac) | |
| 125 { | |
| 126 octave_value_list tmp = daspk_jac->do_multi_index_op (1, args); | |
| 127 | |
| 128 if (error_state) | |
| 129 { | |
| 130 gripe_user_supplied_eval ("daspk"); | |
| 131 return retval; | |
| 132 } | |
| 133 | |
| 134 int tlen = tmp.length (); | |
| 135 if (tlen > 0 && tmp(0).is_defined ()) | |
| 136 { | |
| 137 if (! warned_jac_imaginary && tmp(0).is_complex_type ()) | |
| 138 { | |
| 139 warning ("daspk: ignoring imaginary part returned from user-supplied jacobian function"); | |
| 140 warned_jac_imaginary = true; | |
| 141 } | |
| 142 | |
| 143 retval = tmp(0).matrix_value (); | |
| 144 | |
| 145 if (error_state || retval.length () == 0) | |
| 146 gripe_user_supplied_eval ("daspk"); | |
| 147 } | |
| 148 else | |
| 149 gripe_user_supplied_eval ("daspk"); | |
| 150 } | |
| 151 | |
| 152 return retval; | |
| 153 } | |
| 154 | |
| 3912 | 155 #define DASPK_ABORT() \ |
| 156 do \ | |
| 157 { \ | |
| 158 unwind_protect::run_frame ("Fdaspk"); \ | |
| 159 return retval; \ | |
| 160 } \ | |
| 161 while (0) | |
| 162 | |
| 163 #define DASPK_ABORT1(msg) \ | |
| 164 do \ | |
| 165 { \ | |
| 166 ::error ("daspk: " msg); \ | |
| 167 DASPK_ABORT (); \ | |
| 168 } \ | |
| 169 while (0) | |
| 170 | |
| 171 #define DASPK_ABORT2(fmt, arg) \ | |
| 172 do \ | |
| 173 { \ | |
| 174 ::error ("daspk: " fmt, arg); \ | |
| 175 DASPK_ABORT (); \ | |
| 176 } \ | |
| 177 while (0) | |
| 178 | |
| 3997 | 179 DEFUN_DLD (daspk, args, nargout, |
| 3912 | 180 "-*- texinfo -*-\n\ |
| 4115 | 181 @deftypefn {Loadable Function} {[@var{x}, @var{xdot}, @var{istate}, @var{msg}] =} daspk (@var{fcn}, @var{x_0}, @var{xdot_0}, @var{t}, @var{t_crit})\n\ |
| 182 Solve the set of differential-algebraic equations\n\ | |
| 183 @tex\n\ | |
| 4852 | 184 $$ 0 = f (x, \\dot{x}, t) $$\n\ |
| 4115 | 185 with\n\ |
| 186 $$ x(t_0) = x_0, \\dot{x}(t_0) = \\dot{x}_0 $$\n\ | |
| 187 @end tex\n\ | |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7017
diff
changeset
|
188 @ifnottex\n\ |
| 4115 | 189 \n\ |
| 190 @example\n\ | |
| 4698 | 191 0 = f (x, xdot, t)\n\ |
| 4115 | 192 @end example\n\ |
| 193 \n\ | |
| 194 with\n\ | |
| 195 \n\ | |
| 196 @example\n\ | |
| 197 x(t_0) = x_0, xdot(t_0) = xdot_0\n\ | |
| 198 @end example\n\ | |
| 199 \n\ | |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7017
diff
changeset
|
200 @end ifnottex\n\ |
| 4115 | 201 The solution is returned in the matrices @var{x} and @var{xdot},\n\ |
| 202 with each row in the result matrices corresponding to one of the\n\ | |
| 3912 | 203 elements in the vector @var{t}. The first element of @var{t}\n\ |
| 4115 | 204 should be @math{t_0} and correspond to the initial state of the\n\ |
| 205 system @var{x_0} and its derivative @var{xdot_0}, so that the first\n\ | |
| 206 row of the output @var{x} is @var{x_0} and the first row\n\ | |
| 207 of the output @var{xdot} is @var{xdot_0}.\n\ | |
| 3912 | 208 \n\ |
| 8787 | 209 The first argument, @var{fcn}, is a string, inline, or function handle\n\ |
| 210 that names the function @math{f} to call to compute the vector of\n\ | |
| 211 residuals for the set of equations. It must have the form\n\ | |
| 3912 | 212 \n\ |
| 213 @example\n\ | |
| 214 @var{res} = f (@var{x}, @var{xdot}, @var{t})\n\ | |
| 215 @end example\n\ | |
| 216 \n\ | |
| 217 @noindent\n\ | |
| 4115 | 218 in which @var{x}, @var{xdot}, and @var{res} are vectors, and @var{t} is a\n\ |
| 3912 | 219 scalar.\n\ |
| 220 \n\ | |
| 8787 | 221 If @var{fcn} is a two-element string array or a two-element cell array\n\ |
|
9067
8970b4b10e9f
Cleanup documentation for quad.texi and diffeq.texi
Rik <rdrider0-list@yahoo.com>
parents:
9064
diff
changeset
|
222 of strings, inline functions, or function handles, the first element names\n\ |
| 8787 | 223 the function @math{f} described above, and the second element names a\n\ |
| 224 function to compute the modified Jacobian\n\ | |
| 4115 | 225 @tex\n\ |
| 226 $$\n\ | |
| 227 J = {\\partial f \\over \\partial x}\n\ | |
| 228 + c {\\partial f \\over \\partial \\dot{x}}\n\ | |
| 229 $$\n\ | |
| 230 @end tex\n\ | |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7017
diff
changeset
|
231 @ifnottex\n\ |
| 4144 | 232 \n\ |
| 233 @example\n\ | |
|
9064
7c02ec148a3c
Check grammar on all .cc files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
234 @group\n\ |
| 4115 | 235 df df\n\ |
| 236 jac = -- + c ------\n\ | |
| 237 dx d xdot\n\ | |
|
9064
7c02ec148a3c
Check grammar on all .cc files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
238 @end group\n\ |
| 4115 | 239 @end example\n\ |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7017
diff
changeset
|
240 @end ifnottex\n\ |
| 4115 | 241 \n\ |
| 242 The modified Jacobian function must have the form\n\ | |
| 243 \n\ | |
| 244 @example\n\ | |
|
9064
7c02ec148a3c
Check grammar on all .cc files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
245 @group\n\ |
| 4115 | 246 \n\ |
| 247 @var{jac} = j (@var{x}, @var{xdot}, @var{t}, @var{c})\n\ | |
| 248 \n\ | |
|
9064
7c02ec148a3c
Check grammar on all .cc files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
249 @end group\n\ |
| 4115 | 250 @end example\n\ |
| 251 \n\ | |
| 3912 | 252 The second and third arguments to @code{daspk} specify the initial\n\ |
| 253 condition of the states and their derivatives, and the fourth argument\n\ | |
| 4115 | 254 specifies a vector of output times at which the solution is desired,\n\ |
| 3912 | 255 including the time corresponding to the initial condition.\n\ |
| 256 \n\ | |
| 257 The set of initial states and derivatives are not strictly required to\n\ | |
| 4115 | 258 be consistent. If they are not consistent, you must use the\n\ |
| 259 @code{daspk_options} function to provide additional information so\n\ | |
| 260 that @code{daspk} can compute a consistent starting point.\n\ | |
| 3912 | 261 \n\ |
| 262 The fifth argument is optional, and may be used to specify a set of\n\ | |
| 263 times that the DAE solver should not integrate past. It is useful for\n\ | |
| 264 avoiding difficulties with singularities and points where there is a\n\ | |
| 265 discontinuity in the derivative.\n\ | |
| 3964 | 266 \n\ |
| 4115 | 267 After a successful computation, the value of @var{istate} will be\n\ |
| 268 greater than zero (consistent with the Fortran version of @sc{Daspk}).\n\ | |
| 269 \n\ | |
| 270 If the computation is not successful, the value of @var{istate} will be\n\ | |
| 271 less than zero and @var{msg} will contain additional information.\n\ | |
| 272 \n\ | |
| 3964 | 273 You can use the function @code{daspk_options} to set optional\n\ |
| 274 parameters for @code{daspk}.\n\ | |
| 5642 | 275 @seealso{dassl}\n\ |
| 276 @end deftypefn") | |
| 3912 | 277 { |
| 278 octave_value_list retval; | |
| 279 | |
| 4140 | 280 warned_fcn_imaginary = false; |
| 281 warned_jac_imaginary = false; | |
| 282 | |
| 3912 | 283 unwind_protect::begin_frame ("Fdaspk"); |
| 284 | |
| 285 unwind_protect_int (call_depth); | |
| 286 call_depth++; | |
| 287 | |
| 288 if (call_depth > 1) | |
| 289 DASPK_ABORT1 ("invalid recursive call"); | |
| 290 | |
| 291 int nargin = args.length (); | |
| 292 | |
| 293 if (nargin > 3 && nargin < 6) | |
| 294 { | |
| 5729 | 295 std::string fcn_name, fname, jac_name, jname; |
| 4140 | 296 daspk_fcn = 0; |
| 297 daspk_jac = 0; | |
| 298 | |
| 299 octave_value f_arg = args(0); | |
| 300 | |
| 5729 | 301 if (f_arg.is_cell ()) |
| 302 { | |
| 303 Cell c = f_arg.cell_value (); | |
| 304 if (c.length() == 1) | |
| 305 f_arg = c(0); | |
| 306 else if (c.length() == 2) | |
| 307 { | |
| 308 if (c(0).is_function_handle () || c(0).is_inline_function ()) | |
| 309 daspk_fcn = c(0).function_value (); | |
| 310 else | |
| 311 { | |
| 312 fcn_name = unique_symbol_name ("__daspk_fcn__"); | |
| 313 fname = "function y = "; | |
| 314 fname.append (fcn_name); | |
| 315 fname.append (" (x, xdot, t) y = "); | |
| 316 daspk_fcn = extract_function | |
| 317 (c(0), "daspk", fcn_name, fname, "; endfunction"); | |
| 318 } | |
| 319 | |
| 320 if (daspk_fcn) | |
| 321 { | |
| 322 if (c(1).is_function_handle () || c(1).is_inline_function ()) | |
| 323 daspk_jac = c(1).function_value (); | |
| 324 else | |
| 325 { | |
| 326 jac_name = unique_symbol_name ("__daspk_jac__"); | |
| 327 jname = "function jac = "; | |
| 328 jname.append(jac_name); | |
| 329 jname.append (" (x, xdot, t, cj) jac = "); | |
| 330 daspk_jac = extract_function | |
| 331 (c(1), "daspk", jac_name, jname, "; endfunction"); | |
| 4140 | 332 |
| 5729 | 333 if (!daspk_jac) |
| 334 { | |
| 335 if (fcn_name.length()) | |
| 336 clear_function (fcn_name); | |
| 337 daspk_fcn = 0; | |
| 338 } | |
| 339 } | |
| 340 } | |
| 341 } | |
| 342 else | |
| 343 DASPK_ABORT1 ("incorrect number of elements in cell array"); | |
| 344 } | |
| 3912 | 345 |
| 5729 | 346 if (!daspk_fcn && ! f_arg.is_cell()) |
| 347 { | |
| 348 if (f_arg.is_function_handle () || f_arg.is_inline_function ()) | |
| 349 daspk_fcn = f_arg.function_value (); | |
| 350 else | |
| 351 { | |
| 352 switch (f_arg.rows ()) | |
| 353 { | |
| 354 case 1: | |
| 355 do | |
| 356 { | |
| 357 fcn_name = unique_symbol_name ("__daspk_fcn__"); | |
| 358 fname = "function y = "; | |
| 359 fname.append (fcn_name); | |
| 360 fname.append (" (x, xdot, t) y = "); | |
| 361 daspk_fcn = extract_function | |
| 362 (f_arg, "daspk", fcn_name, fname, "; endfunction"); | |
| 363 } | |
| 364 while (0); | |
| 365 break; | |
| 4140 | 366 |
| 5729 | 367 case 2: |
| 4140 | 368 { |
| 5729 | 369 string_vector tmp = f_arg.all_strings (); |
| 370 | |
| 371 if (! error_state) | |
| 372 { | |
| 373 fcn_name = unique_symbol_name ("__daspk_fcn__"); | |
| 374 fname = "function y = "; | |
| 375 fname.append (fcn_name); | |
| 376 fname.append (" (x, xdot, t) y = "); | |
| 377 daspk_fcn = extract_function | |
| 378 (tmp(0), "daspk", fcn_name, fname, "; endfunction"); | |
| 4140 | 379 |
| 5729 | 380 if (daspk_fcn) |
| 381 { | |
| 382 jac_name = unique_symbol_name ("__daspk_jac__"); | |
| 383 jname = "function jac = "; | |
| 384 jname.append(jac_name); | |
| 385 jname.append (" (x, xdot, t, cj) jac = "); | |
| 386 daspk_jac = extract_function | |
| 387 (tmp(1), "daspk", jac_name, jname, | |
| 388 "; endfunction"); | |
| 389 | |
| 390 if (!daspk_jac) | |
| 391 { | |
| 392 if (fcn_name.length()) | |
| 393 clear_function (fcn_name); | |
| 394 daspk_fcn = 0; | |
| 395 } | |
| 396 } | |
| 397 } | |
| 4140 | 398 } |
| 5729 | 399 } |
| 400 } | |
| 4140 | 401 } |
| 402 | |
| 403 if (error_state || ! daspk_fcn) | |
| 3912 | 404 DASPK_ABORT (); |
| 405 | |
| 406 ColumnVector state = ColumnVector (args(1).vector_value ()); | |
| 407 | |
| 408 if (error_state) | |
| 409 DASPK_ABORT1 ("expecting state vector as second argument"); | |
| 410 | |
| 411 ColumnVector deriv (args(2).vector_value ()); | |
| 412 | |
| 413 if (error_state) | |
| 414 DASPK_ABORT1 ("expecting derivative vector as third argument"); | |
| 415 | |
| 416 ColumnVector out_times (args(3).vector_value ()); | |
| 417 | |
| 418 if (error_state) | |
| 419 DASPK_ABORT1 ("expecting output time vector as fourth argument"); | |
| 420 | |
| 421 ColumnVector crit_times; | |
| 422 int crit_times_set = 0; | |
| 423 if (nargin > 4) | |
| 424 { | |
| 425 crit_times = ColumnVector (args(4).vector_value ()); | |
| 426 | |
| 427 if (error_state) | |
| 428 DASPK_ABORT1 ("expecting critical time vector as fifth argument"); | |
| 429 | |
| 430 crit_times_set = 1; | |
| 431 } | |
| 432 | |
| 433 if (state.capacity () != deriv.capacity ()) | |
| 434 DASPK_ABORT1 ("x and xdot must have the same size"); | |
| 435 | |
| 436 double tzero = out_times (0); | |
| 437 | |
| 438 DAEFunc func (daspk_user_function); | |
| 4140 | 439 if (daspk_jac) |
| 440 func.set_jacobian_function (daspk_user_jacobian); | |
| 441 | |
| 3912 | 442 DASPK dae (state, deriv, tzero, func); |
| 4122 | 443 dae.set_options (daspk_opts); |
| 3912 | 444 |
| 445 Matrix output; | |
| 446 Matrix deriv_output; | |
| 447 | |
| 448 if (crit_times_set) | |
| 449 output = dae.integrate (out_times, deriv_output, crit_times); | |
| 450 else | |
| 451 output = dae.integrate (out_times, deriv_output); | |
| 452 | |
| 5729 | 453 if (fcn_name.length()) |
| 454 clear_function (fcn_name); | |
| 455 if (jac_name.length()) | |
| 456 clear_function (jac_name); | |
| 457 | |
| 3912 | 458 if (! error_state) |
| 459 { | |
| 3997 | 460 std::string msg = dae.error_message (); |
| 461 | |
| 462 retval(3) = msg; | |
| 463 retval(2) = static_cast<double> (dae.integration_state ()); | |
| 3912 | 464 |
| 3997 | 465 if (dae.integration_ok ()) |
| 466 { | |
| 467 retval(1) = deriv_output; | |
| 468 retval(0) = output; | |
| 469 } | |
| 470 else | |
| 471 { | |
| 472 retval(1) = Matrix (); | |
| 473 retval(0) = Matrix (); | |
| 474 | |
| 475 if (nargout < 3) | |
| 476 error ("daspk: %s", msg.c_str ()); | |
| 477 } | |
| 3912 | 478 } |
| 479 } | |
| 480 else | |
| 5823 | 481 print_usage (); |
| 3912 | 482 |
| 483 unwind_protect::run_frame ("Fdaspk"); | |
| 484 | |
| 485 return retval; | |
| 486 } | |
| 487 | |
| 488 /* | |
| 489 ;;; Local Variables: *** | |
| 490 ;;; mode: C++ *** | |
| 491 ;;; End: *** | |
| 492 */ |