find numerical solution of differential-algebraic initial value problems - Maple Help

Online Help

All Products    Maple    MapleSim


Home : Support : Online Help : Mathematics : Factorization and Solving Equations : Numerical Solutions : dsolve/numeric/DAE

dsolve/numeric/DAE - find numerical solution of differential-algebraic initial value problems

Calling Sequence

dsolve(daesys, numeric, vars, options)

Parameters

daesys

-

set or list; ordinary differential equation(s), algebraic equation(s) and initial conditions

numeric

-

name; instruct dsolve to find a numerical solution

vars

-

(optional) dependent variable or a set or list of dependent variables for daesys

options

-

(optional) equations of the form keyword = value

Description

• 

The dsolve command with the numeric or type=numeric option and a real-valued differential-algebraic initial value problem (DAE IVP) finds a numerical solution for the DAE IVP. If the optional equation method=numericmethod is provided (where numericmethod is one of rkf45_dae, ck45_dae, rosenbrock_dae, or mebdfi), dsolve uses that method to obtain the numerical solution.

  

In most cases dsolve is able to detect if the problem is a DAE system, as opposed to an ODE system, namely the cases in which pure algebraic equations in the dependent variables are present. If the input is a DAE system containing no purely algebraic equations, the method must be included to specify that the system is a DAE system.

• 

Constrained mechanical systems often give rise to DAE problems. (See the pendulum example below.)

• 

The return value of dsolve and the following high level or common options are discussed in dsolve[numeric] and dsolve[Error_Control].

'output'

=

keyword or array

'range'

=

numeric..numeric

'abserr'

=

numeric or list

'relerr'

=

numeric

'stiff'

=

boolean

'known'

=

name or list of names

'optimize'

=

boolean

  

The exception is that for the DAE solvers the absolute error tolerance can be specified as a per-component list. For expected behavior, the variables of the problem must also be specified as a list, and the entries of 'abserr' must correspond 1-1 to the variables in the list, or to the variables in the system converted to first order using the order of the variables in the list. For more information, see the ?dsolve[numeric,IVP] help page.

• 

The default DAE IVP method is a modified Runge-Kutta Fehlberg method, which uses a base order 4-5 method, but has been modified to find solutions for DAE problems. The default stiff method is a Rosenbrock method, which uses a base order 3-4 method. For a description of the modifications done to these methods in extending them to DAE solution, see the ?dae_extension help page. The other method available for DAE IVP is the dsolve[mebdfi] method, which is short for Modified Extended Backward-Differentiation Formula Implicit method.

• 

In general, the DAE IVP solvers are very similar to the standard differential IVP solvers, so this page is primarily concerned with outlining the differences between them.

• 

The DAE solvers are currently restricted to finding solutions for real-valued problems.

• 

For use of any of the methods, the specified initial conditions must satisfy all hidden constraints of the system (that is, they must be a consistent set of initial conditions with respect to the DAE). In the event that they do not, an error results, and information is provided on the unsatisfied condition.

  

In some cases, it may be necessary to use fsolve to compute consistent initial conditions for the problem.

• 

If the ?dae_extension methods are in use, the differential option is set to true, and the system is sufficiently linear in the algebraic variables (i.e., variables which have no derivatives appearing in the input system), it is possible to skip initial conditions for those variables. If the initial conditions are skipped when they are required, an error will be produced.

• 

The following options are also available for some or all of the DAE methods:

'minstep'

=

numeric

'maxstep'

=

numeric

'initstep'

=

numeric

'startinit'

=

boolean

'events'

=

list

'event_pre'

=

keyword

'event_maxiter'

=

integer

'projection'

=

boolean

'differential'

=

boolean

'implicit'

=

boolean

'parameters'

=

list of names

  

'minstep', 'maxstep', and 'initstep'

  

These options are discussed in dsolve[Error_Control], and are only available for the dsolve[mebdfi] method.

  

 

  

'startinit'= boolean

  

This option controls the behavior of the numerical integration with respect to consecutive calculations. This option is described in dsolve[numeric,IVP].

  

For the methods, rkf45_dae, ck45_dae, and rosenbrock_dae, when a 'range' has been specified, the solution is recomputed only when the new value of the independent variable is not between the initial point and any other previously requested solution point. This has the effect of never reversing the direction of integration, and making evaluation of the solution for an already computed interval quite inexpensive (the solution values are obtained by interpolation). The storage of the solution can also be enabled by using the storage argument. The startinit parameter also forces these methods to recompute the solution whenever a solution value is desired.

  

 

  

'events'= list

  

'event_pre'= keyword

  

'event_maxiter'= integer

  

These options are available only for the rkf45_dae, ck45_dae, and rosenbrock_dae methods, and are used to specify and control the behavior of events. These are the same as for standard IVP problems. For a description, see the ?dsolve[numeric,IVP] and ?dsolve[Events] help pages.

  

 

  

'projection', 'differential', and 'implicit'

  

The 'projection', 'differential' and 'implicit' options are specific to the extension methods, so they are discussed there.

  

 

• 

In addition to the computation of solution values for the given input problem, the procedure returns (that is, output=procedurelist, listprocedure or operator) provide additional interactive features, including the ability to specify parameters. Information on these features is provided on the dsolve[numeric,interactive] page. The exception is that the mebdfi method cannot work with parametrized problems.

Examples

As a first example, we consider the problem of modeling the dynamics of a mass on a string of unit length in 2-D Cartesian coordinates (the pendulum problem). We let r be the position of the mass on the string, and v the velocity:

VectorCalculus[SetCoordinates](cartesian):

VectorCalculus[BasisFormat](false):

r := VectorCalculus[Vector]([x(t),y(t)]);

r:=xtyt

(1)

v := VectorCalculus[diff](r,t);

v:=ⅆⅆtxtⅆⅆtyt

(2)

Since the mass is on a string of fixed length l, we have the constraint:

con := VectorCalculus[DotProduct](r,r)-l^2;

con:=xt2+yt2l2

(3)

Now we want to construct the DAE system using the Euler-Lagrange formulation, so we compute the kinetic energy T and the potential energy V as:

T := m/2*VectorCalculus[DotProduct](v,v);

T:=12mⅆⅆtxt2+ⅆⅆtyt2

(4)

V := m*g*y(t);

V:=mgyt

(5)

where m is the mass, and g is the gravitational constant (we will use 9.8). The Lagrangian and modified Lagrangian are given by:

L := T-V;

L:=12mⅆⅆtxt2+ⅆⅆtyt2mgyt

(6)

ML := L - lambda(t)*con;

ML:=12mⅆⅆtxt2+ⅆⅆtyt2mgytλtxt2+yt2l2

(7)

We can then construct the Euler-Lagrange formulation via:

Q := remove(has,VariationalCalculus:-EulerLagrange(
                 ML,t,[x(t),y(t),lambda(t)]),K);

Q:=2λtxtmⅆ2ⅆt2xt,l2xt2yt2,mg2ytλtmⅆ2ⅆt2yt

(8)

EQx := isolate(select(has,Q,diff(x(t),t,t))[1],diff(x(t),t,t));

EQx:=ⅆ2ⅆt2xt=2λtxtm

(9)

EQy := isolate(select(has,Q,diff(y(t),t,t))[1],diff(y(t),t,t));

EQy:=ⅆ2ⅆt2yt=mg+2ytλtm

(10)

EQc := remove(has,Q,diff)[1];

EQc:=l2xt2yt2

(11)

Now we have the equations of motion for the pendulum. Next, we need to determine consistent initial conditions. To do so, we must identify any hidden constraints of the system. These are easy to find, as we have only one constraint.

Dcon := diff(EQc,t);

Dcon:=2xtⅆⅆtxt2ytⅆⅆtyt

(12)

DDcon := eval(diff(Dcon,t),{EQx,EQy});

DDcon:=2ⅆⅆtxt2+4xt2λtm2ⅆⅆtyt2+2ytmg+2ytλtm

(13)

Our initial conditions must satisfy EQc, Dcon, and DDcon at the initial point, leaving only 2 degrees of freedom for the conditions. So for a pendulum starting at the minimum value of y0=l having an initial horizontal velocity of Dx0=vx, we get:

sys := {y(0)=-l,D(x)(0)=vx} union
       eval(convert({EQc,Dcon,DDcon},D),t=0);

sys:=2x0Dx02y0Dy0,l2x02y02,2Dx02+4x02λ0m2Dy02+2y0mg+2y0λ0m,y0=l,Dx0=vx

(14)

ini := [solve(sys,{x(0),y(0),D(x)(0),D(y)(0),lambda(0)})][1];

ini:=λ0=12mgl+vx2l2,x0=0,y0=l,Dx0=vx,Dy0=0

(15)

So we consider the above with a pendulum of unit length l=1 having unit mass m=1 and an initial horizontal velocity of vx=110, giving us the DAE system and initial conditions:

dsys := eval({EQx,EQy,EQc},{l=1,m=1,g=9.8,vx=1/10});

dsys:=1xt2yt2,ⅆ2ⅆt2xt=2λtxt,ⅆ2ⅆt2yt=9.82ytλt

(16)

dini := eval(ini,{l=1,m=1,g=9.8,vx=1/10});

dini:=λ0=4.905000000,x0=0,y0=1,Dx0=110,Dy0=0

(17)

We can then obtain the solution as:

dsol1 := dsolve(dsys union dini, numeric);

dsol1:=procx_rkf45_dae...end proc

(18)

dsol1(1/2);

t=0.500000000000000,λt=4.89750023467855,xt=0.0319392699415791,ⅆⅆtxt=0.000564539406239432,yt=0.999489811490763,ⅆⅆtyt=0.0000180384378593043

(19)

dsol1(1);

t=1.,λt=4.90499905723651,xt=0.000360891517794992,ⅆⅆtxt=0.0999937504894365,yt=0.999999934755133,ⅆⅆtyt=0.0000360794101890703

(20)

Solution with rosenbrock_dae:

dsol2 := dsolve(dsys union dini, numeric, method=rosenbrock_dae);

dsol2:=procx_rosenbrock_dae...end proc

(21)

dsol2(1/2);

t=0.500000000000000,λt=4.89750019118586,xt=0.0319392364763694,ⅆⅆtxt=0.000564626004686083,yt=0.999489813543334,ⅆⅆtyt=0.0000180482728717227

(22)

dsol2(1);

t=1.,λt=4.90499906207811,xt=0.000360968463263150,ⅆⅆtxt=0.0999936779452163,yt=0.999999936691318,ⅆⅆtyt=0.0000360897908606990

(23)

Solution with mebdfi:

dsol3 := dsolve(dsys union dini, numeric, method=mebdfi);

dsol3:=procx_mebdfi...end proc

(24)

dsol3(1/2);

t=0.500000000000000,λt=4.89749908804820855,xt=0.0319389748677520596,ⅆⅆtxt=0.000565199414553933343,yt=0.999489820801930273,ⅆⅆtyt=0.0000180581765179766362

(25)

dsol3(1);

t=1.,λt=4.90499581413915475,xt=0.000361181747066772996,ⅆⅆtxt=0.0999926673290293111,yt=0.99999993019115307,ⅆⅆtyt=0.0000360957473722881540

(26)

Now consider a similar problem as above, but in addition add a second mass supported from the first by another string, this one of length 1/2 (the double pendulum). The system can be obtained and solved as:

dsysd := {
   diff(x1(t),t,t)    +2*lambda1(t)*x1(t)+2*lambda2(t)*(x1(t)-x2(t)),
   diff(y1(t),t,t)+9.8+2*lambda1(t)*y1(t)+2*lambda2(t)*(y1(t)-y2(t)),
   diff(x2(t),t,t)                       -2*lambda2(t)*(x1(t)-x2(t)),
   diff(y2(t),t,t)+9.8                   -2*lambda2(t)*(y1(t)-y2(t)),
   x1(t)^2+y1(t)^2-1,
   (x1(t)-x2(t))^2+(y1(t)-y2(t))^2-1/4
};

dsysd:=ⅆ2ⅆt2x2t2λ2tx1tx2t,x1tx2t2+y1ty2t214,x1t2+y1t21,ⅆ2ⅆt2x1t+2λ1tx1t+2λ2tx1tx2t,ⅆ2ⅆt2y2t+9.82λ2ty1ty2t,ⅆ2ⅆt2y1t+9.8+2λ1ty1t+2λ2ty1ty2t

(27)

ics := {x1(0)=0,D(x1)(0)=-3,y1(0)=-1,D(y1)(0)=0,
        x2(0)=0,D(x2)(0)=4,y2(0)=-3/2,D(y2)(0)=0};

ics:=x10=0,x20=0,y10=1,y20=32,Dx10=3,Dx20=4,Dy10=0,Dy20=0

(28)

dsold := dsolve(dsysd union ics,numeric);

dsold:=procx_rkf45_dae...end proc

(29)

The trajectory of the second mass can be plotted via:

plots[odeplot](dsold,[x2(t),y2(t)],0..5,numpoints=1000);

Note that we did not specify initial conditions for lambda1t and lambda2t as we were using an extension method, and the system was sufficiently linear in the lambda.

See Also

dsolve/dae_extension, dsolve/Error_Control, dsolve[ck45], dsolve[Events], dsolve[maxfun], dsolve[mebdfi], dsolve[numeric,interactive], dsolve[numeric,IVP], dsolve[numeric], dsolve[rkf45], dsolve[rosenbrock], dsolve[Stiffness], plots[odeplot]


Download Help Document

Was this information helpful?



Please add your Comment (Optional)
E-mail Address (Optional)
What is ? This question helps us to combat spam