Solving DAE using daspk function

classic Classic list List threaded Threaded
4 messages Options
vmz
Reply | Threaded
Open this post in threaded view
|

Solving DAE using daspk function

vmz
hello,

i would like to solve 3 DAE examples from  this paper
<https://pdfs.semanticscholar.org/c750/77a37def9367e373a50e1570d69ee97fdc8f.pdf>  
using daspk function. first two of them seems to work ok. the last one
doesn't. here is my code

# exam1
exam1x = [1.0; 0.0];
exam1xdot = [1.0; 0.0];
exam1xav = [0,1];
function f = exam1f (x, xdot, t)
        f(1) = xdot(1) - t*cos(t) + x(1) - (1 + t) * x(2);
        f(2) = sin(t) - x(2);
endfunction

# exam2
exam2x = [1.0; 1.0];
exam2xdot = [1.0; 0.0];
exam2xav = [0,1];
function f = exam2f (x, xdot, t)
        f(1) = xdot(1) - x(2);
        f(2) = x(2)^3.0 - x(1)^2;
endfunction

# exam3
exam3x = [5.0; 1.0; -1.0; 0.0];
exam3xdot = [1.0; 0.0; 0.0; 1.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

# solve
function solve (f, x0, xdot0, xav, t)
        daspk_options("algebraic variables", xav);
        daspk_options("initial step size", 0.1);
        [x, xdot] = daspk (f, x0, xdot0, t);
        printf("%10s %10s %10s\n", "t", "x0", "x1");
        for i = 1:length(t)
                printf("%10.4f %10.4f %10.4f\n", t(i), x(i,1), x(i,2));
        endfor
endfunction

# solve
t = linspace(0,10,6);
solve("exam1f", exam1x, exam1xdot, exam1xav, t);
solve("exam2f", exam2x, exam2xdot, exam2xav, t);
solve("exam3f", exam3x, exam3xdot, exam3xav, t);

i'd appreciate any help / advice. thanks in advance.
vit mach-zizka






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


Reply | Threaded
Open this post in threaded view
|

Re: Solving DAE using daspk function

Windhorn, Allen E [ACIM/LSA/MKT]
vmz,

I thank you for your interesting question of the 20th, because I need to learn
to solve differential-algebraic equations.  I'm not particularly familiar with
daspk yet, so I don't know why your third example isn't working, but you can
solve the third and fourth equations (the algebraic ones) for x(3) and x(4) by
hand, and just substitute those values back into the first two equations,
turning the system into a pure differential equation (please check my work to
see if I got the equations correct):

function f = exam4f(x, xdot, t)
  f(1) = xdot(1)+t*x(2)+(1.0+t)*(x(2)-5*sin(t^2/2));
  f(2) = xdot(2)+t*x(1)+(1.0+t)*5*(x(1)-5*cos(t^2/2));
endfunction

Still using daspk (because it was handy), I tried this function and got a solution.
The solution gets really truly enormous very quickly, so I don't know whether
to believe it or not, but at least it doesn't crash.

Regards,
Allen
--
Allen Windhorn P.E. (Mn), CEng| Senior Principal Engineer
Leroy-Somer Americas | Kato Engineering, Inc.
2075 Howard Dr. West | North Mankato, MN 56003 | USA
T +1 507-345-2782 | F +1 507-345-2798
[hidden email] | [hidden email]



vmz
Reply | Threaded
Open this post in threaded view
|

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


Reply | Threaded
Open this post in threaded view
|

RE: Solving DAE using daspk function

Windhorn, Allen E [ACIM/LSA/MKT]
vmz,

> -----Original Message-----
> From: Help-octave <help-octave-bounces+allen.windhorn=[hidden email]> On Behalf Of vmz
>
> 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 ...

I'm glad you got it working -- thanks for posting the result, as I have an
interest in this.

Regards,
Allen