From bug-request at octave dot org Sun Feb  6 05:18:34 2005
Subject: is daspk working?
From: Claude =?iso-8859-1?Q?Lacoursi=E8re?= <claude at hpc2n dot umu dot se>
To: bug-octave at bevo dot che dot wisc dot edu
Date: Sun, 6 Feb 2005 12:22:56 +0100

Hi. 

>From the documentation, it would appear that the daspk
interface is identical to that of dassl which it supersedes.
However, though I can get daspk to integrate a simple ODE, I
cannot sucessfully integrate a simple DAE which dassl can
process properly. 

Anyone aware of issues with the daspk solver?  I cannot find
a sample script and there is nothing relevant in the list
archives. 

Currently using version 2.1.64

When I run the script below using daspk, I get: 
DASPK--  AT T (=R1) AND STEPSIZE H (=R2) THE
      In above,  R1 =  0.1000000000000E-01   R2 =
      0.9094947017729E-17
DASPK--  ERROR TEST FAILED REPEATEDLY OR WITH ABS(H)=HMIN

but I get a correct answer using dassl, whether or not I
provide my own Jacobian function. 

Best regards, 
Claude.

===============simple script follows===================

function X = bead ()

  ## usage:  X = bead (h, n)
  ##
  ## integrate 2D particle falling under gravity but constrained
  ## to move on a wire i.e., the line : x + y = 0. 
  ## Mass is set to 1, so is the constrant downward gravity force. 
  
  h = 0.01; 			# time step for solution vector
  n = 100; 			# number of sample points

  func = char("b", "bjac");	# use explicit Jacobian evaluation
  ##func = char("b");		# use daspk's numerical Jacobian evaluation

  ## start everything at rest at the origin.
  x0 = zeros(5, 1); v0 = zeros(5, 1); 
  T = h*(1:n)';
  [X, V, IS, MSG ] = daspk (func, x0, v0, T); 
  disp(IS);
  if ( IS  < 0 )
    disp(MSG);
  endif
  ## build plottable solution vector
  X = [X, T];
endfunction

## Evaluate residual.  The equations of motion are second order so we do
## a mapping to a first order equation by introducing the velocities as
## extra variables. 
function r = b (x, v, t)
  y1 = [0; 1];			# unit gravity pointing downwards
  y2 = [1; 1];			# the constraint is a line at 45 degrees
  r(1:2, 1) = v(1:2, 1)  - x(3:4, 1); 
  r(3:4, 1) = v(3:4, 1)  + y1 + x(5, 1)*y2; 
  r(5  , 1) = y2'*x(1:2, 1) ; 
endfunction


## Evaluate the Jacobian. 
function J = bjac(x, v, t, c)
  J = zeros(5, 5);
  y2 = [1; 1];
  J(1:2, 3:4) = -eye(2, 2); 
  J(3:4, 5) = y2; 
  J(5, 1:2) = y2';
  J2 = eye(5, 5);  J2(5, 5)  = 0 ; 
  J  = J + c*J2; 
endfunction

==================script end==========================



-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------