help-octave
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: Solving DAE using daspk function


From: vmz
Subject: Re: Solving DAE using daspk function
Date: Thu, 5 Dec 2019 11:06:27 -0600 (CST)

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



reply via email to

[Prev in Thread] Current Thread [Next in Thread]