dae_root

Arguments

y0

a real vector or matrix (initial conditions).

t0

a real scalar (initial time).

t

a real vector (times at which the solution is computed).

f

an external i.e. function or character string or list, computes the value of f(t, y, ydot).

rtol, atol

a real scalar or a column vector of the same size as y0.

jac

an external i.e. function or character string or list, computes the value of dg/dx+cj*dg/dxdot for a given value of parameter cj.

Syntax of jac: r = jac(t, y, ydot, cj)

ng

an integer.

g

an external i.e. function or character string or list, computes the value of the column vector g(t, x) with ng components. Each component defines a surface.

psol

an external i.e. function or character string or list, solves a linear system P*x = b, with P being the factored preconditioner that routine pjac computed beforehand and stored in wp and iwp.

Syntax of psol: [r, ier] = psol(wp, iwp, b)

pjac

an external i.e. function or character string or list, computes the value of dg/dy + cj*dg/dydot for a given value of parameter cj and LU-factorizes it in two arrays, real and integer.

Syntax of pjac: [wp, iwp, ires] = pjac(neq, t, y, ydot, h, cj, rewt, savr)

hd

real vector which allows to store the dae context and to resume integration.

y

a real vector or matrix. The solution.

nn

a real vector with two entries [times num], times is the value of the time at which the surface is crossed, num is the number of the crossed surface.

Description

With this syntax (first argument equal to "root" or "root2") dae solves the implicit differential equation:

            g(t,y,ydot) = 0
            y(t0) = y0  and   ydot(t0) = ydot0
            

Returns the surface crossing instants and the number of the surface reached in nn.

Other arguments and other options are the same as for dae, see the dae help.

Examples

Example #1: use "root"

// dy/dt = ((2*log(y)+8)/t -5)*y,  y(1) = 1,  1<=t<=6
// g1 = ((2*log(y)+8)/t - 5)*y
// g2 = log(y) - 2.2491
y0 = 1; t = 2:6; t0 = 1; y0d = 3;
atol = 1.d-6; rtol = 0; ng = 2;

deff('[delta,ires] = res1(t,y,ydot)', 'ires=0; delta=ydot-((2*log(y)+8)/t-5)*y')
deff('rts = gr1(t,y)', 'rts=[((2*log(y)+8)/t-5)*y;log(y)-2.2491]')

[yy,nn] = dae("root", [y0,y0d],t0,t,rtol,atol,res1,ng,gr1);
//(Should return nn=[2.4698972 2])

Example #2: use "root2"

// dy1/dt = y2
// dy2/dt = 100 * (1 - y1^2) * y2 - y1
// g = y1
t0 = 0;
y0 = [2;0];
y0d = [0; -2];
t = [20:20:200];
ng = 1;
rtol = [1.d-6; 1.d-6];
atol = [1.d-6; 1.d-4];

deff("[delta,ires]=res2(t,y,ydot)",...
"ires=0;y1=y(1),y2=y(2),delta=[ydot-[y2;100*(1-y1*y1)*y2-y1]]")

deff("s=gr2(t,y,yd)","s=y(1)")

[yy, nn]=dae("root2", [y0, y0d], t0, t, rtol, atol, res2, ng, gr2);
nn

// Should return nn = [81.163512 1]

See also

• dae — Differential algebraic equations solver