hi allen,

thank you for the answer.

i've figured this out. there were two issues

1) analytical jacobian needs to be supplied to solve succesfully.

2) two sign typos in the paper need to be corrected (in second and fourth

equation of the third example).

so the final working code exactly matching the numerical examples from

section V of the paper looks like

# exam1

exam1x = [1.0; 0.0];

exam1xdot = [1.0; 0.0];

exam1xav = [0; 1];

function f = exam1f (x, xdot, t)

f(1) = - t*cos(t) + x(1) - (1.0 + t) * x(2) + xdot(1);

f(2) = sin(t) - x(2);

endfunction

function j = exam1jac (x, xdot, t, c)

j(1,1) = 1.0 + c;

j(1,2) = -(1.0+t);

j(2,1) = 0.0;

j(2,2) = -1.0;

endfunction

# exam2

exam2x = [1.0; 1.0];

exam2xdot = [1.0; 0.0];

exam2xav = [0; 1];

function f = exam2f (x, xdot, t)

f(1) = - x(2) + xdot(1);

f(2) = - x(1)^2.0 + x(2)^3.0;

endfunction

function j = exam2jac (x, xdot, t, c)

j(1,1) = c;

j(1,2) = -1.0;

j(2,1) = -2.0*x(1);

j(2,2) = +3.0*x(2)^2.0;

endfunction

# exam3

exam3x = [5.0; 1.0; -1.0; 0.0];

exam3xdot = [1.0; 0.0; 0.0; 0.0];

exam3xav = [0;0;1;1];

function f = exam3f (x, xdot, t)

f(1) = xdot(1) + t*x(2) + (1.0 + t)*x(3);

f(2) = xdot(2) - t*x(1) + (1.0 + t)*x(4);

f(3) = (x(1) - x(4)) / 5.0 - cos(t^2.0/2.0);

f(4) = (x(2) + x(3)) / 5.0 - sin(t^2.0/2.0);

endfunction

function j = exam3jac (x, xdot, t, c)

j(1,1) = c;

j(1,2) = t;

j(1,3) = 1.0 + t;

j(1,4) = 0.0;

j(2,1) = -t;

j(2,2) = c;

j(2,3) = 0.0;

j(2,4) = 1.0 + t;

j(3,1) = 1.0/5.0;

j(3,2) = 0.0;

j(3,3) = 0.0;

j(3,4) = -1.0/5.0;

j(4,1) = 0.0;

j(4,2) = 1.0/5.0;

j(4,3) = 1.0/5.0;

j(4,4) = 0.0;

endfunction

# solve

function solve (hdr, fcn, x0, xdot0, xav, t)

daspk_options("algebraic variables", xav);

daspk_options("exclude algebraic variables from error test", 1);

daspk_options("initial step size", 0.00001);

[x, xdot] = daspk (fcn, x0, xdot0, t);

printf("%s\n", hdr);

for i = 1:length(t)

printf("%10.4f ", t(i));

for j = 1:size(x,2)

printf("%10.4f ", x(i,j));

endfor

printf("\n");

endfor

endfunction

# solve

t = linspace(0,10,6);

solve("Exam 1", {"exam1f"; "exam1jac"}, exam1x, exam1xdot, exam1xav, t);

solve("Exam 2", {"exam2f"; "exam2jac"}, exam2x, exam2xdot, exam2xav, t);

solve("Exam 3", {"exam3f"; "exam3jac"}, exam3x, exam3xdot, exam3xav, t);

hope this could help someone, best regards.

vit mach-zizka

--

Sent from:

https://octave.1599824.n4.nabble.com/Octave-General-f1599825.html