Re: Solving DAE using daspk function

[vmz]

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