dae_root
differential algebraic equation solver with roots finding
Syntax
[y, nn, [,hd]] = dae("root", y0, t0, t, f, ng, g) [y, nn, [,hd]] = dae("root", y0, t0, t [,rtol [,atol]], f [,jac], ng, g [,hd]) [y, nn, [,hd]] = dae("root2", y0, t0, t, f, ng, g) [y, nn, [,hd]] = dae("root2", y0, t0, t [,rtol [,atol]], f [,jac], ng, g [,psol, pjac] [,hd])
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
Report an issue | ||
<< dae | Дифференциальные уравнения | daeoptions >> |