Mercurial > hg > octave-lojdl
annotate src/DLD-FUNCTIONS/lsode.cc @ 8920:eb63fbe60fab
update copyright notices
| author | John W. Eaton <jwe@octave.org> |
|---|---|
| date | Sat, 07 Mar 2009 10:41:27 -0500 |
| parents | 5b23faa8113c |
| children | 7c02ec148a3c |
| rev | line source |
|---|---|
| 2928 | 1 /* |
| 2 | |
| 7017 | 3 Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2005, |
| 8920 | 4 2006, 2007, 2008, 2009 John W. Eaton |
| 2928 | 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. | |
| 2928 | 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/>. | |
| 2928 | 21 |
| 22 */ | |
| 23 | |
| 24 #ifdef HAVE_CONFIG_H | |
| 25 #include <config.h> | |
| 26 #endif | |
| 27 | |
| 28 #include <string> | |
| 29 | |
| 3567 | 30 #include <iomanip> |
| 3523 | 31 #include <iostream> |
| 2928 | 32 |
| 33 #include "LSODE.h" | |
| 34 #include "lo-mappers.h" | |
| 35 | |
| 36 #include "defun-dld.h" | |
| 37 #include "error.h" | |
| 38 #include "gripes.h" | |
| 39 #include "oct-obj.h" | |
| 2968 | 40 #include "ov-fcn.h" |
| 5729 | 41 #include "ov-cell.h" |
| 2928 | 42 #include "pager.h" |
| 3952 | 43 #include "pr-output.h" |
| 3243 | 44 #include "unwind-prot.h" |
| 2928 | 45 #include "utils.h" |
| 46 #include "variables.h" | |
| 47 | |
| 3998 | 48 #include "LSODE-opts.cc" |
| 49 | |
| 2928 | 50 // Global pointer for user defined function required by lsode. |
| 2968 | 51 static octave_function *lsode_fcn; |
| 2928 | 52 |
| 53 // Global pointer for optional user defined jacobian function used by lsode. | |
| 2968 | 54 static octave_function *lsode_jac; |
| 2928 | 55 |
| 4140 | 56 // Have we warned about imaginary values returned from user function? |
| 57 static bool warned_fcn_imaginary = false; | |
| 58 static bool warned_jac_imaginary = false; | |
| 59 | |
| 3243 | 60 // Is this a recursive call? |
| 61 static int call_depth = 0; | |
| 62 | |
| 2928 | 63 ColumnVector |
| 64 lsode_user_function (const ColumnVector& x, double t) | |
| 65 { | |
| 66 ColumnVector retval; | |
| 67 | |
| 68 octave_value_list args; | |
| 69 args(1) = t; | |
| 4628 | 70 args(0) = x; |
| 2928 | 71 |
| 72 if (lsode_fcn) | |
| 73 { | |
| 3544 | 74 octave_value_list tmp = lsode_fcn->do_multi_index_op (1, args); |
| 2928 | 75 |
| 76 if (error_state) | |
| 77 { | |
| 78 gripe_user_supplied_eval ("lsode"); | |
| 79 return retval; | |
| 80 } | |
| 81 | |
| 82 if (tmp.length () > 0 && tmp(0).is_defined ()) | |
| 83 { | |
| 4140 | 84 if (! warned_fcn_imaginary && tmp(0).is_complex_type ()) |
| 85 { | |
| 86 warning ("lsode: ignoring imaginary part returned from user-supplied function"); | |
| 87 warned_fcn_imaginary = true; | |
| 88 } | |
| 89 | |
| 3419 | 90 retval = ColumnVector (tmp(0).vector_value ()); |
| 2928 | 91 |
| 92 if (error_state || retval.length () == 0) | |
| 93 gripe_user_supplied_eval ("lsode"); | |
| 94 } | |
| 95 else | |
| 96 gripe_user_supplied_eval ("lsode"); | |
| 97 } | |
| 98 | |
| 99 return retval; | |
| 100 } | |
| 101 | |
| 102 Matrix | |
| 103 lsode_user_jacobian (const ColumnVector& x, double t) | |
| 104 { | |
| 105 Matrix retval; | |
| 106 | |
| 107 octave_value_list args; | |
| 108 args(1) = t; | |
| 4628 | 109 args(0) = x; |
| 2928 | 110 |
| 111 if (lsode_jac) | |
| 112 { | |
| 3544 | 113 octave_value_list tmp = lsode_jac->do_multi_index_op (1, args); |
| 2928 | 114 |
| 115 if (error_state) | |
| 116 { | |
| 117 gripe_user_supplied_eval ("lsode"); | |
| 118 return retval; | |
| 119 } | |
| 120 | |
| 121 if (tmp.length () > 0 && tmp(0).is_defined ()) | |
| 122 { | |
| 4140 | 123 if (! warned_jac_imaginary && tmp(0).is_complex_type ()) |
| 124 { | |
| 125 warning ("lsode: ignoring imaginary part returned from user-supplied jacobian function"); | |
| 126 warned_jac_imaginary = true; | |
| 127 } | |
| 128 | |
| 2928 | 129 retval = tmp(0).matrix_value (); |
| 130 | |
| 131 if (error_state || retval.length () == 0) | |
| 132 gripe_user_supplied_eval ("lsode"); | |
| 133 } | |
| 134 else | |
| 135 gripe_user_supplied_eval ("lsode"); | |
| 136 } | |
| 137 | |
| 138 return retval; | |
| 139 } | |
| 140 | |
| 3323 | 141 #define LSODE_ABORT() \ |
| 142 do \ | |
| 143 { \ | |
| 144 unwind_protect::run_frame ("Flsode"); \ | |
| 145 return retval; \ | |
| 146 } \ | |
| 147 while (0) | |
| 148 | |
| 149 #define LSODE_ABORT1(msg) \ | |
| 150 do \ | |
| 151 { \ | |
| 3747 | 152 ::error ("lsode: " msg); \ |
| 3323 | 153 LSODE_ABORT (); \ |
| 154 } \ | |
| 155 while (0) | |
| 156 | |
| 157 #define LSODE_ABORT2(fmt, arg) \ | |
| 158 do \ | |
| 159 { \ | |
| 3747 | 160 ::error ("lsode: " fmt, arg); \ |
| 3323 | 161 LSODE_ABORT (); \ |
| 162 } \ | |
| 163 while (0) | |
| 164 | |
| 2928 | 165 DEFUN_DLD (lsode, args, nargout, |
| 3373 | 166 "-*- texinfo -*-\n\ |
| 6755 | 167 @deftypefn {Loadable Function} {[@var{x}, @var{istate}, @var{msg}] =} lsode (@var{fcn}, @var{x_0}, @var{t}, @var{t_crit})\n\ |
| 4115 | 168 Solve the set of differential equations\n\ |
| 169 @tex\n\ | |
| 170 $$ {dx \\over dt} = f (x, t) $$\n\ | |
| 171 with\n\ | |
| 172 $$ x(t_0) = x_0 $$\n\ | |
| 173 @end tex\n\ | |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7562
diff
changeset
|
174 @ifnottex\n\ |
| 4115 | 175 \n\ |
| 176 @example\n\ | |
| 177 dx\n\ | |
| 178 -- = f(x, t)\n\ | |
| 4117 | 179 dt\n\ |
| 4115 | 180 @end example\n\ |
| 181 \n\ | |
| 182 with\n\ | |
| 183 \n\ | |
| 184 @example\n\ | |
| 185 x(t_0) = x_0\n\ | |
| 186 @end example\n\ | |
| 187 \n\ | |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7562
diff
changeset
|
188 @end ifnottex\n\ |
| 4115 | 189 The solution is returned in the matrix @var{x}, with each row\n\ |
| 190 corresponding to an element of the vector @var{t}. The first element\n\ | |
| 191 of @var{t} should be @math{t_0} and should correspond to the initial\n\ | |
| 192 state of the system @var{x_0}, so that the first row of the output\n\ | |
| 193 is @var{x_0}.\n\ | |
| 3373 | 194 \n\ |
| 8787 | 195 The first argument, @var{fcn}, is a string, inline, or function handle\n\ |
| 196 that names the function @math{f} to call to compute the vector of right\n\ | |
| 197 hand sides for the set of equations. The function must have the form\n\ | |
| 2928 | 198 \n\ |
| 3373 | 199 @example\n\ |
| 200 @var{xdot} = f (@var{x}, @var{t})\n\ | |
| 201 @end example\n\ | |
| 2928 | 202 \n\ |
| 3373 | 203 @noindent\n\ |
| 4115 | 204 in which @var{xdot} and @var{x} are vectors and @var{t} is a scalar.\n\ |
| 205 \n\ | |
| 8787 | 206 If @var{fcn} is a two-element string array or a two-element cell array\n\ |
| 207 of strings, inlines, or function handles, the first element names\n\ | |
| 208 the function @math{f} described above, and the second element names a\n\ | |
| 209 function to compute the Jacobian of @math{f}. The Jacobian function\n\ | |
| 210 must have the form\n\ | |
| 4115 | 211 \n\ |
| 212 @example\n\ | |
| 213 @var{jac} = j (@var{x}, @var{t})\n\ | |
| 214 @end example\n\ | |
| 215 \n\ | |
| 216 in which @var{jac} is the matrix of partial derivatives\n\ | |
| 217 @tex\n\ | |
| 218 $$ J = {\\partial f_i \\over \\partial x_j} = \\left[\\matrix{\n\ | |
| 219 {\\partial f_1 \\over \\partial x_1}\n\ | |
| 220 & {\\partial f_1 \\over \\partial x_2}\n\ | |
| 221 & \\cdots\n\ | |
| 222 & {\\partial f_1 \\over \\partial x_N} \\cr\n\ | |
| 223 {\\partial f_2 \\over \\partial x_1}\n\ | |
| 224 & {\\partial f_2 \\over \\partial x_2}\n\ | |
| 225 & \\cdots\n\ | |
| 226 & {\\partial f_2 \\over \\partial x_N} \\cr\n\ | |
| 227 \\vdots & \\vdots & \\ddots & \\vdots \\cr\n\ | |
| 228 {\\partial f_3 \\over \\partial x_1}\n\ | |
| 229 & {\\partial f_3 \\over \\partial x_2}\n\ | |
| 230 & \\cdots\n\ | |
| 231 & {\\partial f_3 \\over \\partial x_N} \\cr}\\right]$$\n\ | |
| 232 @end tex\n\ | |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7562
diff
changeset
|
233 @ifnottex\n\ |
| 4115 | 234 \n\ |
| 235 @example\n\ | |
| 236 | df_1 df_1 df_1 |\n\ | |
| 237 | ---- ---- ... ---- |\n\ | |
| 238 | dx_1 dx_2 dx_N |\n\ | |
| 239 | |\n\ | |
| 240 | df_2 df_2 df_2 |\n\ | |
| 241 | ---- ---- ... ---- |\n\ | |
| 242 df_i | dx_1 dx_2 dx_N |\n\ | |
| 243 jac = ---- = | |\n\ | |
| 244 dx_j | . . . . |\n\ | |
| 245 | . . . . |\n\ | |
| 246 | . . . . |\n\ | |
| 247 | |\n\ | |
| 248 | df_N df_N df_N |\n\ | |
| 249 | ---- ---- ... ---- |\n\ | |
| 250 | dx_1 dx_2 dx_N |\n\ | |
| 251 @end example\n\ | |
| 252 \n\ | |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7562
diff
changeset
|
253 @end ifnottex\n\ |
| 4115 | 254 \n\ |
| 7001 | 255 The second and third arguments specify the initial state of the system,\n\ |
| 4115 | 256 @math{x_0}, and the initial value of the independent variable @math{t_0}.\n\ |
| 2928 | 257 \n\ |
| 3373 | 258 The fourth argument is optional, and may be used to specify a set of\n\ |
| 259 times that the ODE solver should not integrate past. It is useful for\n\ | |
| 260 avoiding difficulties with singularities and points where there is a\n\ | |
| 261 discontinuity in the derivative.\n\ | |
| 3964 | 262 \n\ |
| 4115 | 263 After a successful computation, the value of @var{istate} will be 2\n\ |
| 264 (consistent with the Fortran version of @sc{Lsode}).\n\ | |
| 265 \n\ | |
| 266 If the computation is not successful, @var{istate} will be something\n\ | |
| 267 other than 2 and @var{msg} will contain additional information.\n\ | |
| 268 \n\ | |
| 3964 | 269 You can use the function @code{lsode_options} to set optional\n\ |
| 270 parameters for @code{lsode}.\n\ | |
| 5694 | 271 @seealso{daspk, dassl, dasrt}\n\ |
| 5646 | 272 @end deftypefn") |
| 2928 | 273 { |
| 274 octave_value_list retval; | |
| 275 | |
| 4140 | 276 warned_fcn_imaginary = false; |
| 277 warned_jac_imaginary = false; | |
| 278 | |
| 3243 | 279 unwind_protect::begin_frame ("Flsode"); |
| 2928 | 280 |
| 3243 | 281 unwind_protect_int (call_depth); |
| 282 call_depth++; | |
| 283 | |
| 284 if (call_depth > 1) | |
| 3323 | 285 LSODE_ABORT1 ("invalid recursive call"); |
| 2928 | 286 |
| 3243 | 287 int nargin = args.length (); |
| 2928 | 288 |
| 3959 | 289 if (nargin > 2 && nargin < 5 && nargout < 4) |
| 2928 | 290 { |
| 5729 | 291 std::string fcn_name, fname, jac_name, jname; |
| 3991 | 292 lsode_fcn = 0; |
| 293 lsode_jac = 0; | |
| 294 | |
| 3243 | 295 octave_value f_arg = args(0); |
| 2928 | 296 |
| 5729 | 297 if (f_arg.is_cell ()) |
| 298 { | |
| 299 Cell c = f_arg.cell_value (); | |
| 300 if (c.length() == 1) | |
| 301 f_arg = c(0); | |
| 302 else if (c.length() == 2) | |
| 303 { | |
| 304 if (c(0).is_function_handle () || c(0).is_inline_function ()) | |
| 305 lsode_fcn = c(0).function_value (); | |
| 306 else | |
| 307 { | |
| 308 fcn_name = unique_symbol_name ("__lsode_fcn__"); | |
| 309 fname = "function y = "; | |
| 310 fname.append (fcn_name); | |
| 311 fname.append (" (x, t) y = "); | |
| 312 lsode_fcn = extract_function | |
| 313 (c(0), "lsode", fcn_name, fname, "; endfunction"); | |
| 314 } | |
| 315 | |
| 316 if (lsode_fcn) | |
| 317 { | |
| 318 if (c(1).is_function_handle () || c(1).is_inline_function ()) | |
| 319 lsode_jac = c(1).function_value (); | |
| 320 else | |
| 321 { | |
| 322 jac_name = unique_symbol_name ("__lsode_jac__"); | |
| 323 jname = "function jac = "; | |
| 324 jname.append(jac_name); | |
| 325 jname.append (" (x, t) jac = "); | |
| 326 lsode_jac = extract_function | |
| 327 (c(1), "lsode", jac_name, jname, "; endfunction"); | |
| 2928 | 328 |
| 5729 | 329 if (!lsode_jac) |
| 330 { | |
| 331 if (fcn_name.length()) | |
| 332 clear_function (fcn_name); | |
| 333 lsode_fcn = 0; | |
| 334 } | |
| 335 } | |
| 336 } | |
| 337 } | |
| 338 else | |
| 339 LSODE_ABORT1 ("incorrect number of elements in cell array"); | |
| 340 } | |
| 2928 | 341 |
| 5729 | 342 if (!lsode_fcn && ! f_arg.is_cell()) |
| 343 { | |
| 344 if (f_arg.is_function_handle () || f_arg.is_inline_function ()) | |
| 345 lsode_fcn = f_arg.function_value (); | |
| 346 else | |
| 347 { | |
| 348 switch (f_arg.rows ()) | |
| 349 { | |
| 350 case 1: | |
| 351 do | |
| 352 { | |
| 353 fcn_name = unique_symbol_name ("__lsode_fcn__"); | |
| 354 fname = "function y = "; | |
| 355 fname.append (fcn_name); | |
| 356 fname.append (" (x, t) y = "); | |
| 357 lsode_fcn = extract_function | |
| 358 (f_arg, "lsode", fcn_name, fname, "; endfunction"); | |
| 359 } | |
| 360 while (0); | |
| 361 break; | |
| 2928 | 362 |
| 5729 | 363 case 2: |
| 3243 | 364 { |
| 5729 | 365 string_vector tmp = f_arg.all_strings (); |
| 366 | |
| 367 if (! error_state) | |
| 368 { | |
| 369 fcn_name = unique_symbol_name ("__lsode_fcn__"); | |
| 370 fname = "function y = "; | |
| 371 fname.append (fcn_name); | |
| 372 fname.append (" (x, t) y = "); | |
| 373 lsode_fcn = extract_function | |
| 374 (tmp(0), "lsode", fcn_name, fname, "; endfunction"); | |
| 3243 | 375 |
| 5729 | 376 if (lsode_fcn) |
| 377 { | |
| 378 jac_name = unique_symbol_name ("__lsode_jac__"); | |
| 379 jname = "function jac = "; | |
| 380 jname.append(jac_name); | |
| 381 jname.append (" (x, t) jac = "); | |
| 382 lsode_jac = extract_function | |
| 383 (tmp(1), "lsode", jac_name, jname, | |
| 384 "; endfunction"); | |
| 385 | |
| 386 if (!lsode_jac) | |
| 387 { | |
| 388 if (fcn_name.length()) | |
| 389 clear_function (fcn_name); | |
| 390 lsode_fcn = 0; | |
| 391 } | |
| 392 } | |
| 393 } | |
| 3243 | 394 } |
| 5729 | 395 break; |
| 2928 | 396 |
| 5729 | 397 default: |
| 398 LSODE_ABORT1 | |
| 399 ("first arg should be a string or 2-element string array"); | |
| 400 } | |
| 401 } | |
| 3243 | 402 } |
| 2928 | 403 |
| 3243 | 404 if (error_state || ! lsode_fcn) |
| 3323 | 405 LSODE_ABORT (); |
| 2928 | 406 |
| 3419 | 407 ColumnVector state (args(1).vector_value ()); |
| 2928 | 408 |
| 409 if (error_state) | |
| 3323 | 410 LSODE_ABORT1 ("expecting state vector as second argument"); |
| 3243 | 411 |
| 3419 | 412 ColumnVector out_times (args(2).vector_value ()); |
| 3243 | 413 |
| 414 if (error_state) | |
| 3323 | 415 LSODE_ABORT1 ("expecting output time vector as third argument"); |
| 2928 | 416 |
| 3243 | 417 ColumnVector crit_times; |
| 2928 | 418 |
| 3243 | 419 int crit_times_set = 0; |
| 420 if (nargin > 3) | |
| 421 { | |
| 3419 | 422 crit_times = ColumnVector (args(3).vector_value ()); |
| 2928 | 423 |
| 3243 | 424 if (error_state) |
| 3323 | 425 LSODE_ABORT1 ("expecting critical time vector as fourth argument"); |
| 2928 | 426 |
| 3243 | 427 crit_times_set = 1; |
| 428 } | |
| 2928 | 429 |
| 3243 | 430 double tzero = out_times (0); |
| 431 | |
| 432 ODEFunc func (lsode_user_function); | |
| 433 if (lsode_jac) | |
| 434 func.set_jacobian_function (lsode_user_jacobian); | |
| 2928 | 435 |
| 3243 | 436 LSODE ode (state, tzero, func); |
| 437 | |
| 4122 | 438 ode.set_options (lsode_opts); |
| 3243 | 439 |
| 3859 | 440 Matrix output; |
| 3243 | 441 if (crit_times_set) |
| 442 output = ode.integrate (out_times, crit_times); | |
| 443 else | |
| 444 output = ode.integrate (out_times); | |
| 445 | |
| 5729 | 446 if (fcn_name.length()) |
| 447 clear_function (fcn_name); | |
| 448 if (jac_name.length()) | |
| 449 clear_function (jac_name); | |
| 450 | |
| 3243 | 451 if (! error_state) |
| 452 { | |
| 3997 | 453 std::string msg = ode.error_message (); |
| 454 | |
| 455 retval(2) = msg; | |
| 3959 | 456 retval(1) = static_cast<double> (ode.integration_state ()); |
| 457 | |
| 458 if (ode.integration_ok ()) | |
| 3971 | 459 retval(0) = output; |
| 3959 | 460 else |
| 461 { | |
| 3971 | 462 retval(0) = Matrix (); |
| 3959 | 463 |
| 464 if (nargout < 2) | |
| 3997 | 465 error ("lsode: %s", msg.c_str ()); |
| 3959 | 466 } |
| 3243 | 467 } |
| 468 } | |
| 2928 | 469 else |
| 5823 | 470 print_usage (); |
| 2928 | 471 |
| 3243 | 472 unwind_protect::run_frame ("Flsode"); |
| 2928 | 473 |
| 474 return retval; | |
| 475 } | |
| 476 | |
| 477 /* | |
|
7562
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
478 |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
479 %% dassl-1.m |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
480 %% |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
481 %% Test lsode() function |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
482 %% |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
483 %% Author: David Billinghurst (David.Billinghurst@riotinto.com.au) |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
484 %% Comalco Research and Technology |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
485 %% 20 May 1998 |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
486 %% |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
487 %% Problem |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
488 %% |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
489 %% y1' = -y2, y1(0) = 1 |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
490 %% y2' = y1, y2(0) = 0 |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
491 %% |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
492 %% Solution |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
493 %% |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
494 %% y1(t) = cos(t) |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
495 %% y2(t) = sin(t) |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
496 %!function xdot = f (x, t) |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
497 %! xdot = [-x(2); x(1)]; |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
498 %!test |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
499 %! |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
500 %! x0 = [1; 0]; |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
501 %! xdot0 = [0; 1]; |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
502 %! t = (0:1:10)'; |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
503 %! |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
504 %! tol = 500 * lsode_options ("relative tolerance"); |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
505 %! |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
506 %! |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
507 %! x = lsode ("f", x0, t); |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
508 %! |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
509 %! y = [cos(t), sin(t)]; |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
510 %! |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
511 %! assert(all (all (abs (x - y) < tol))); |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
512 |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
513 %!function xdotdot = f (x, t) |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
514 %! xdotdot = [x(2); -x(1)]; |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
515 %!test |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
516 %! |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
517 %! x0 = [1; 0]; |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
518 %! t = [0; 2*pi]; |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
519 %! tol = 100 * dassl_options ("relative tolerance"); |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
520 %! |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
521 %! x = lsode ("f", x0, t); |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
522 %! |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
523 %! y = [1, 0; 1, 0]; |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
524 %! |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
525 %! assert(all (all (abs (x - y) < tol))); |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
526 |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
527 %!function xdot = f (x, t) |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
528 %! xdot = x; |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
529 %!test |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
530 %! |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
531 %! x0 = 1; |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
532 %! t = [0; 1]; |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
533 %! tol = 100 * dassl_options ("relative tolerance"); |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
534 %! |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
535 %! x = lsode ("f", x0, t); |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
536 %! |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
537 %! y = [1; e]; |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
538 %! |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
539 %! assert(all (all (abs (x - y) < tol))); |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
540 |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
541 %!test |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
542 %! lsode_options ("absolute tolerance", eps); |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
543 %! assert(lsode_options ("absolute tolerance") == eps); |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
544 |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
545 %!error <Invalid call to lsode_options.*> lsode_options ("foo", 1, 2); |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
546 |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
547 */ |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
548 |
|
c827f5673321
move tests to individual source files
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
549 /* |
| 2928 | 550 ;;; Local Variables: *** |
| 551 ;;; mode: C++ *** | |
| 552 ;;; End: *** | |
| 553 */ |