Differential Equations - Maple Help
For the best experience, we recommend viewing online help using Google Chrome or Microsoft Edge.

Online Help

All Products    Maple    MapleSim


Home : Support : Online Help : System : Information : Updates : Maple 11 : Differential Equations

Updates to Differential Equation (DE) Solvers in Maple 11

 

Summary

Exact Solutions

Numeric Solutions and related Plotting

Summary

  

The differential equation (DE) theme for Maple 11 is the exploration of the geometrical structure underlying a DE. Thus for 2nd order nonlinear Ordinary Differential Equations (ODEs), new algorithms were developed to solve entire classes in terms of elliptic functions.

  

For 1st order ODEs, the nonlinear Abel AIR class, as well as those admitting conformal symmetries, are now all entirely solvable in terms of hypergeometric and elliptic functions.

  

For Partial Differential Equations (PDEs) the geometrical structure is also reflected in their symmetries, which is the subject of this release, together with the related group invariant solutions, representing the most important advancement of the last 5 years in the Maple libraries for the exact solution of PDEs.

Exact Solutions

Ordinary Differential Equations (ODEs)

New solutions in terms of elliptic functions for 2nd order nonlinear equations

  

New solving algorithms for various classes of nonlinear 2nd order ODEs admitting solutions expressible in terms of the WeierstrassP or any of the 12 JacobiPQ elliptic functions. These classes of equations admit only one point symmetry and the related reduction of order leads to 1st order nonlinear Abel ODEs, a problem for which Maple 11 also presents new algorithms.

  

Examples

PDEtools:-declare(y(x), prime=x);

yxwill now be displayed asy

derivatives with respect toxof functions of one variable will now be displayed with '

(1)

diff(y(x), x, x) = 4/x*diff(y(x), x) + 6/x^2*y(x)^2 - 6/x^2*y(x) - 1/2*x^2*a;

y''=4y'x+6y2x26yx2x2a2

(2)

dsolve((2));

y=WeierstrassPx+c__1,a,c__2x2

(3)

diff(y(x), x, x) = 2/x*a/(x+a)*diff(y(x), x) + 2*k^2/x^2/(x+a)^4*y(x)^3 + (-2*x^3*a+(-k^2-1-6*a^2)*x^2 - 6*a^3*x - 2*a^4)/x^2/(x+a)^4*y(x);

(4)

dsolve((4));

(5)
  

The algorithm underlying the output above can compute solutions free of integrals whenever they are linear in WeierstrassP or any of the 12 JacobiPQ functions.

New option singsol = <none, essential, all> and new elliptic function solutions for high degree 1st order nonlinear ODE

  

High-degree equations are nonlinear in the highest derivative. In addition to their general solution, these equations admit singular solutions automatically computed by dsolve.

  

The novelty in Maple 11 is twofold. You can now specify whether to compute all, none, or only the essential singular solutions by passing an extra argument singsol. Also, instead of the standard general solution composing fractional powers with RootOf of elliptic functions or uncomputed integrals, dsolve now computes these general solutions with no RootOfs and no uncomputed integrals but directly as explicit solutions linear in JacobiPQ, WeierstrassP, and WeierstrassPPrime elliptic functions when that is possible.

  

Examples

  

Two examples quadratic in y'. To see the general solution, pass only singsol=none (it also saves the time dsolve would spend splitting the problem into cases and computing the singular solutions).

(-2/x^3*y(x)+1/x^2*diff(y(x),x))^2-1-b^2/x^8*y(x)^4-(-b^2-1)/x^4*y(x)^2 = 0;

2yx3+y'x221b2y4x8b21y2x4=0

(6)

dsolve((6), singsol = none);

y=JacobiSNx+c__1&comma;bx2

(7)
  

The singular solutions for these high-degree ODEs are frequently the roots of some related high-degree polynomial.  To see this polynomial representing an implicit singular solution instead of their roots, use the optional argument implicit.

(-2/x^3*y(x)+1/x^2*diff(y(x),x))^2+1/x^2*y(x)*a-4/x^6*y(x)^3+b = 0;

2yx3+y'x22+yax24y3x6+b=0

(8)

dsolve((8), implicit);    # implicit avoids computing the roots of the cube

y3bx64yax44=0,y=WeierstrassPx+c__1&comma;a&comma;bx2

(9)
  

An example cubic in y' where the solution involves WeierstrassPPrime

-(x+a)^6*(y(x) + x*diff(y(x),x))^3 + 1/2*a^3/b^3 - 27/2*x^4*y(x)^4 - 27*b/c*x^2*y(x)^2 - 27/2*b^2/c^2 - 3/2*a/b*(x+a)^4*(y(x)+x*diff(y(x),x))^2 = 0;

x+a6y+y'x3+a32b327x4y4227bx2y2c27b22c23ax+a4y+y'x22b=0

(10)

[dsolve((10), implicit)];

y454b4x2y2c+a3c227b527x4b3c2=0&comma;y=WeierstrassPPrime1x+a+c__1&comma;ab&comma;bcx

(11)

map(odetest, (11), (10));  # verify these results for correctness

0&comma;0

(12)

New solutions in terms of hypergeometric functions for 1st order Abel equations of the AIR class

  

The new solving algorithms for first order Abel ODEs of the 3-parameter AIR class can solve the entire class of equations having rational coefficients by resolving an equivalence to the seed representative (see references) under changes of variables of the form x -> F(x), y(x) -> (P(x) y(x) + Q(x))/(R(x) y(x) + S(x)), for rational mappings F, P, Q, R, S. In addition, for restricted subclasses, the algorithms implemented can compute the solution also when these mappings, and the coefficients in the corresponding ODE, are non-rational. The Abel AIR class is presented as entirely solvable for the first time in the literature now with Maple 11.

  

Examples

  

Two examples involving non-rational coefficients

diff(y(x),x) = y(x)/x/(ln(x)*y(x)+(-ln(x)+b)*(-ln(x)+c));

y'=yxlnxy+lnx+blnx+c

(13)

dsolve((13));

c__1+lnxKummerMb&comma;1+cb&comma;ybKummerMb+1&comma;1+cb&comma;ylnxKummerUb&comma;1+cb&comma;y+bcKummerUb+1&comma;1+cb&comma;y=0

(14)
  

An example whose implicit solution can be obtained in terms of HeunC

diff(y(x),x) = y(x)*(y(x)-1)*exp(x)/(y(x)*a+exp(2*x)-exp(x)*c);

y'=yy1&ExponentialE;xya+&ExponentialE;2x&ExponentialE;xc

(15)

sol := dsolve((15));

solc__1+HeunCPrime0&comma;c&comma;−1&comma;0&comma;a+12&comma;yy1y&ExponentialE;2xy1HeunC0&comma;c&comma;−1&comma;0&comma;a+12&comma;yy1&ExponentialE;xy1cy1&ExponentialE;2x&ExponentialE;x+cHeunC0&comma;c&comma;−1&comma;0&comma;a+12&comma;yy1+yHeunCPrime0&comma;c&comma;−1&comma;0&comma;a+12&comma;yy1yc=0

(16)
  

To see this solution in terms of hypergeometric functions enter convert(sol, hypergeom);.  For more details, see convert/to_special_function. An example of this class where the solution is expressed directly in terms of hypergeometric functions

diff(y(x),x) = ((3+3*x^4+3*x^2-6*x^3+6*x)*y(x)^2 + (3*x^5+30*x^4-26*x^6-3*x^3-15*x^7+2*x^8)*y(x) - x^7*(29*x^3-31-27*x+x^5+21*x^2+5*x^4))/((x^4+x^2+1)*y(x) + x^3*(x^3+x-3*x^2+x^5-1-x^4))/x;

y'=3x46x3+3x2+6x+3y2+2x815x726x6+3x5+30x43x3yx7x5+5x4+29x3+21x227x31x4+x2+1y+x3x5x4+x33x2+x1x

(17)

dsolve((17));

(18)

New solutions for 1st order equations admitting conformal symmetries

  

In his original symmetry work, Sophus Lie was the first to observe that nonlinear first ODEs with conformal symmetries can have their solutions computed by quadratures. The infinitesimal generators [xi(x, y), eta(x, y)] of a conformal symmetry satisfy the conditions xi[x] = eta[y] and xi[y] = -eta[x], where [x] means differentiation with respect to x. In the framework of the Maple libraries, the key observation is that ODEs possessing this type of symmetry can systematically be transformed into equations admitting infinitesimals separable by product - a problem for which Maple implemented solving algorithms in Maple 6 (see references).

  

Examples

  

Kamke's example number 1.83

diff(y(x),x) = tan(x*y(x));

y'=tanxy

(19)

dsolve((19), implicit);

c__1+erf2x+Iy2erf2xIy2=0

(20)

Enhanced computation of solutions in terms of hypergeometric functions for 3rd and higher order linear equations

  

There are improvements in the solving algorithms for third and higher order linear ODEs to compute solutions expressible in terms of hypergeometric functions.

diff(y(x),x,x,x,x) = -15/2*x/(x^4+3)*diff(y(x),x)-45/2*x^2/(x^4+3)*diff(y(x),x,x)-10*x^3/(x^4+3)*diff(y(x),x,x,x)+15/16*y(x)/(x^4+3);

(21)

dsolve((21));

y=c__1hypergeom38&comma;58&comma;78&comma;98&comma;34&comma;54&comma;32&comma;3x4x32+c__2hypergeom18&comma;38&comma;58&comma;78&comma;12&comma;34&comma;54&comma;3x4x+c__3xhypergeom18&comma;18&comma;38&comma;58&comma;14&comma;12&comma;34&comma;3x4+c__4hypergeom58&comma;78&comma;98&comma;118&comma;54&comma;32&comma;74&comma;3x4x52

(22)

Partial Differential Equations (PDEs): separability, symmetry analysis and group invariant solutions

  

Nineteen new commands were added to PDEtools to perform the steps of the traditional symmetry analysis of PDE systems, as well as to automatically compute infinitesimal symmetry generators and related exact group invariant solutions departing directly from the given PDE system. The new set of commands includes one for reducing one PDE system taking into account another one.  Among the old PDEtools commands, separability has been extended to determining the separability by sum or by product of PDE systems.  Formerly it could determine the separability only of a single PDE.

Separability of PDE systems

  

Example

  

Given a nonlinear PDE system with two unknowns, u, and v, of three independent variables x, y, and z, invoke the declare facility to avoid redundancy in the display of the output.  Derivatives are displayed as indexed objects.

with(PDEtools);

CanonicalCoordinates&comma;ChangeSymmetry&comma;CharacteristicQ&comma;CharacteristicQInvariants&comma;ConservedCurrentTest&comma;ConservedCurrents&comma;ConsistencyTest&comma;D_Dx&comma;DeterminingPDE&comma;Eta_k&comma;Euler&comma;FirstIntegralSolver&comma;FromJet&comma;FunctionFieldSolutions&comma;InfinitesimalGenerator&comma;Infinitesimals&comma;IntegratingFactorTest&comma;IntegratingFactors&comma;InvariantEquation&comma;InvariantSolutions&comma;InvariantTransformation&comma;Invariants&comma;Laplace&comma;Library&comma;PDEplot&comma;PolynomialSolutions&comma;ReducedForm&comma;SimilaritySolutions&comma;SimilarityTransformation&comma;Solve&comma;SymmetryCommutator&comma;SymmetryGauge&comma;SymmetrySolutions&comma;SymmetryTest&comma;SymmetryTransformation&comma;TWSolutions&comma;ToJet&comma;ToMissingDependentVariable&comma;build&comma;casesplit&comma;charstrip&comma;dchange&comma;dcoeffs&comma;declare&comma;diff_table&comma;difforder&comma;dpolyform&comma;dsubs&comma;mapde&comma;separability&comma;splitstrip&comma;splitsys&comma;undeclare

(23)

declare((u,v)(x,y,z));

ux&comma;y&comma;zwill now be displayed asu

vx&comma;y&comma;zwill now be displayed asv

(24)

sys := [diff(u(x,y,z),x,y,z)*v(x,y,z) + diff(u(x,y,z),x,y)*diff(v(x,y,z),z) + diff(u(x,y,z),x,z)*diff(v(x,y,z),y) + diff(u(x,y,z),x)*diff(v(x,y,z),y,z) + diff(u(x,y,z),y,z)*diff(v(x,y,z),x) + diff(u(x,y,z),y)*diff(v(x,y,z),x,z) + diff(u(x,y,z),z)*diff(v(x,y,z),x,y) + u(x,y,z)*diff(v(x,y,z),x,y,z) = 0, diff(u(x,y,z),x,y,z) + diff(v(x,y,z),x,y,z) = 0];

sysux,y,zv+ux,yvz+ux,zvy+uxvy,z+uy,zvx+uyvx,z+uzvx,y+uvx,y,z=0&comma;ux,y,z+vx,y,z=0

(25)
  

As in the case of a single PDE, for PDE systems, unless indicated otherwise, the test is regarding separability by sum.

separability(sys);

0

(26)
  

Note that these separability test results are obtained only considering related integrability conditions, not actually separating the variables. Regarding separation by product,

remain := [separability(sys, `*`)];

remainuxvxuz2v2+vz2u2uy2v2+vy2u2uxuyuzv2+vxvyvzu2&comma;uyvyuz2v2+vz2u2ux2v2+vx2u2uxuyuzv2+vxvyvzu2&comma;uzvzuy2v2+vy2u2ux2v2+vx2u2uxuyuzv2+vxvyvzu2

(27)
  

This result means no separable solution exists unless the three equations above are zero, that is, there is no solution really involving the product of three functions respectively depending on x, y, z. That does not mean no separable solution exists where one or both of the unknowns {u(x, y, z), v(x, y, z)} depend on less variables (that is, is a constant with respect to one or more of x, y, z). For example, consider the following restriction:

restriction := [u(x,y,z) = u(x,z), v(x,y,z) = v(y,z)];

restrictionu=ux&comma;z&comma;v=vy&comma;z

(28)
  

Evaluate remain at this restriction

eval(remain, restriction);

0&comma;0&comma;0

(29)
  

So with this restriction the separability conditions are satisfied. In any case, the solutions separable by `+` and by `*` can be computed  with pdsolve using the HINT option.  For this particular system, only the solutions separable by `+` have u and v depending on the three variables.

pdsolve(sys, HINT = `+`);

u=f__1x+f__2y+f__3z&comma;v=f__4x+f__5y+f__6z

(30)

pdetest((30), sys);    # verify results - see ?pdetest

0&comma;0

(31)

pdsolve(sys, HINT = `*`);

u=f__1xf__2yc__1&comma;v=f__4xf__5yc__2

(32)

pdetest((32), sys);

0&comma;0

(33)

Symmetry analysis and group invariant solutions for PDE systems

  

Example

  

Using the new symmetry commands of PDEtools, the entire symmetry "analysis & solutions" cycle can be performed, automatically or step by step. Consider for instance

declare(u(x,t));

ux&comma;twill now be displayed asu

(34)

PDE := diff(u(x,t),x,x) - diff(u(x,t), t) = 0;

PDEux,xut=0

(35)
  

The determining PDE system whose solutions are the infinitesimals

DetSys := DeterminingPDE(PDE);

DetSys_&xi;tt,t,t=0&comma;_&eta;ut,u=_&xi;tt,t4&comma;_&eta;uu,u=0&comma;_&eta;ux,x=_&eta;ut&comma;_&xi;tu=0&comma;_&xi;tx=0&comma;_&xi;xt=2_&eta;uu,x&comma;_&xi;xu=0&comma;_&xi;xx=_&xi;tt2

(36)
  

You can now compute the infinitesimals of the generators of point symmetry transformations of PDE by solving this system with pdsolve.

pdsolve(DetSys);

_&eta;ux&comma;t&comma;u=x22tc__18+c__4x+c__5u+c__8&ExponentialE;_c1tc__6&ExponentialE;_c1x+c__8&ExponentialE;_c1tc__7&ExponentialE;_c1x&comma;_&xi;tx&comma;t&comma;u=12c__1t2+c__2t+c__3&comma;_&xi;xx&comma;t&comma;u=c__1x4c__4t2+c__2x2+c__9

(37)
  

Alternatively, you can call the new Infinitesimals command with PDE as the argument.

G := Infinitesimals(PDE);

G_&xi;xx&comma;t&comma;u=0&comma;_&xi;tx&comma;t&comma;u=1&comma;_&eta;ux&comma;t&comma;u=0,_&xi;xx&comma;t&comma;u=1&comma;_&xi;tx&comma;t&comma;u=0&comma;_&eta;ux&comma;t&comma;u=0,_&xi;xx&comma;t&comma;u=0&comma;_&xi;tx&comma;t&comma;u=0&comma;_&eta;ux&comma;t&comma;u=u,_&xi;xx&comma;t&comma;u=x2&comma;_&xi;tx&comma;t&comma;u=t&comma;_&eta;ux&comma;t&comma;u=0,_&xi;xx&comma;t&comma;u=2t&comma;_&xi;tx&comma;t&comma;u=0&comma;_&eta;ux&comma;t&comma;u=ux,_&xi;xx&comma;t&comma;u=0&comma;_&xi;tx&comma;t&comma;u=0&comma;_&eta;ux&comma;t&comma;u=&ExponentialE;_c1x&ExponentialE;_c1t,_&xi;xx&comma;t&comma;u=0&comma;_&xi;tx&comma;t&comma;u=0&comma;_&eta;ux&comma;t&comma;u=&ExponentialE;_c1t&ExponentialE;_c1x,_&xi;xx&comma;t&comma;u=tx2&comma;_&xi;tx&comma;t&comma;u=t22&comma;_&eta;ux&comma;t&comma;u=x2+2tu8

(38)
  

These infinitesimals can now be used to compute a transformation leaving PDE invariant.  Take for instance the first infinitesimal, the new variables are v(r, s).

SymmetryTransformation(G[1], u(x,t), v(r,s));

r=x&comma;s=t+_&epsilon;&comma;vr&comma;s=u

(39)
  

Solve this transformation for {x, t, u}

TR := solve((39), {x,t, u(x,t)});

TRt=s_&epsilon;&comma;x=r&comma;u=vr&comma;s

(40)
  

Applying now this transformation to PDE (see dchange) you obtain again PDE, that is, the meaning of "transformation leaving invariant the PDE."

PDE;

ux,xut=0

(41)

dchange(TR, (41), [r,s, v(r,s)]);

vr,rvs=0

(42)
  

These infinitesimals can also be used to compute a transformation reducing the number of independent variables, also called similarity transformation, a particular case of more general group invariant transformations.  For example, take again the first list of infinitesimals.

SimilarityTransformation(G[1], u(x,t), v(r,s), jetnotation = false);

r=x&comma;s=t&comma;vr=u,t=s&comma;x=r&comma;u=vr

(43)
  

In the above, the new variables are v(s). Change variables now in PDE to obtain a reduction in the number of independent variables

dchange((43)[2], PDE);

vr,r=0

(44)
  

Alternatively, the group invariant transformations (and their inverses) reducing the number of independent variables can be obtained for all the symmetries G computed above by calling the InvariantTransformation command with all the symmetries at once.

[InvariantTransformation([G], u(x,t), v(r,s))];

r=t&comma;vr&comma;s=u&comma;t=r&comma;u=vr&comma;s&comma;r=x&comma;vr&comma;s=u&comma;x=r&comma;u=vr&comma;s&comma;r=tx2&comma;s=t&comma;vr&comma;s=u&comma;t=s&comma;x=sr&comma;u=vr&comma;s&comma;r=t&comma;s=x&comma;vr&comma;s=u&ExponentialE;x24t&comma;t=r&comma;x=s&comma;u=vr&comma;s&ExponentialE;s24r&comma;r=tx&comma;s=t&comma;vr&comma;s=ux&ExponentialE;x24t&comma;t=s&comma;x=sr&comma;u=vr&comma;s&ExponentialE;s4r2sr

(45)
  

Change variables for instance using the transformation in the second list above. The transformation, new variables, and resulting reduced PDE are as follows:

(45)[2][2], map(rhs, [op((45)[2][2])]);

x=r&comma;u=vr&comma;s,r&comma;vr&comma;s

(46)

dchange((46)[1], PDE, (46)[2]);

vr,r=0

(47)
  

Solving this transformed PDE (now it is an easier ODE problem) and changing variables back, you obtain the solution to the original PDE. Alternatively, you can perform this entire symmetry cycle automatically, in one step, by directly computing the group invariant solutions of PDE via

[InvariantSolutions(PDE)];

u=c__1&comma;u=c__1x+c__2&comma;u=c__1t&ExponentialE;x24t&comma;u=c__1+erf12tx2c__2&comma;u=c__1x+c__2tttxx&ExponentialE;x24t

(48)
  

This new Maple InvariantSolutions command is unique among the symbolic computation implementations available.  The command automatically reduces, in one step, up to all but one of the independent variables of a PDE system transforming it into an ODE problem.  The command accomplishes this by making simultaneous use of the whole symmetry group, not just one symmetry at a time. It is also possible to provide InvariantSolutions with the dependency expected in the solutions, as well as various other indicators that could make the solution more suitable to your application.

Numeric Solutions and related Plotting

• 

An important number of efficiency improvements have been made for computing numerical solutions of initial value problems for ODEs and differential-algebraic equations (DAEs). These are primarily focused on difficult or large problems.

• 

The DEplot command has been enhanced and given a graphical interface.

Large DAE problems

  

The following example is a rather large and complex system of DAEs for the Stewart platform. This is a platform that is connected to the ground by 6 arms, all of which can expand and contract, and the connections (joints) provide several degrees of freedom.

restart;

stEQs := {z_P(t)-1.*cos(beta_2(t))*s_2(t)-.9659258263*cos(xi_P(t))*sin(eta_P(t
))*cos(zeta_P(t))+.9659258263*sin(xi_P(t))*sin(zeta_P(t))+.2588190451*sin(xi_P
(t))*sin(eta_P(t))*cos(zeta_P(t))+.2588190451*cos(xi_P(t))*sin(zeta_P(t))-.8288000000*cos(beta_2(t)), x_P(t)-1.*sin(beta_2(t))*cos(alpha_2(t))*s_2(t)-1.931851653+.9659258263*cos(xi_P(t))*cos(eta_P(t))-.2588190451*sin(xi_P(t))*cos(


eta_P(t))-.8288000000*sin(beta_2(t))*cos(alpha_2(t)), -.3200000000e-6*(-1.+cos
(beta_1(t))^2)*(2135569.+3125000.*s_1(t)^2+2590000.*s_1(t))*diff(diff(alpha_1(
t),t),t)+(sin(beta_1(t))*sin(alpha_1(t))*s_1(t)+.8288000000*sin(beta_1(t))*sin
(alpha_1(t)))*lambda1(t)+(-1.*sin(beta_1(t))*cos(alpha_1(t))*s_1(t)-.8288000000*sin(beta_1(t))*cos(alpha_1(t)))*lambda2(t)+.6400000000e-6*sin(beta_1(t))*(

1295000.*sin(beta_1(t))*diff(s_1(t),t)+3125000.*s_1(t)^2*cos(beta_1(t))*diff(
beta_1(t),t)+3125000.*sin(beta_1(t))*s_1(t)*diff(s_1(t),t)+2590000.*s_1(t)*cos
(beta_1(t))*diff(beta_1(t),t)+2135569.*cos(beta_1(t))*diff(beta_1(t),t))*diff(
alpha_1(t),t), -.3200000000e-6*(-1.+cos(beta_5(t))^2)*(2590000.*s_5(t)+2135569.+3125000.*s_5(t)^2)*diff(diff(alpha_5(t),t),t)+(sin(beta_5(t))*sin(alpha_5(t)

)*s_5(t)+.8288000000*sin(beta_5(t))*sin(alpha_5(t)))*lambda13(t)+(-1.*sin(
beta_5(t))*cos(alpha_5(t))*s_5(t)-.8288000000*sin(beta_5(t))*cos(alpha_5(t)))*
lambda14(t)+.6400000000e-6*sin(beta_5(t))*(2135569.*cos(beta_5(t))*diff(beta_5
(t),t)+3125000.*s_5(t)^2*cos(beta_5(t))*diff(beta_5(t),t)+3125000.*sin(beta_5(
t))*s_5(t)*diff(s_5(t),t)+2590000.*s_5(t)*cos(beta_5(t))*diff(beta_5(t),t)+
1295000.*sin(beta_5(t))*diff(s_5(t),t))*diff(alpha_5(t),t), -.3200000000e-6*(-1.+cos(beta_4(t))^2)*(2135569.+2590000.*s_4(t)+3125000.*s_4(t)^2)*diff(diff(

alpha_4(t),t),t)+(sin(beta_4(t))*sin(alpha_4(t))*s_4(t)+.8288000000*sin(beta_4
(t))*sin(alpha_4(t)))*lambda10(t)+(-1.*sin(beta_4(t))*cos(alpha_4(t))*s_4(t)-.8288000000*sin(beta_4(t))*cos(alpha_4(t)))*lambda11(t)+.6400000000e-6*sin(

beta_4(t))*(1295000.*sin(beta_4(t))*diff(s_4(t),t)+3125000.*s_4(t)^2*cos(
beta_4(t))*diff(beta_4(t),t)+3125000.*sin(beta_4(t))*s_4(t)*diff(s_4(t),t)+
2590000.*s_4(t)*cos(beta_4(t))*diff(beta_4(t),t)+2135569.*cos(beta_4(t))*diff(
beta_4(t),t))*diff(alpha_4(t),t), -.3200000000e-6*(-1.+cos(beta_3(t))^2)*(
2135569.+2590000.*s_3(t)+3125000.*s_3(t)^2)*diff(diff(alpha_3(t),t),t)+(sin(
beta_3(t))*sin(alpha_3(t))*s_3(t)+.8288000000*sin(beta_3(t))*sin(alpha_3(t)))*
lambda7(t)+(-1.*sin(beta_3(t))*cos(alpha_3(t))*s_3(t)-.8288000000*sin(beta_3(t
))*cos(alpha_3(t)))*lambda8(t)+.6400000000e-6*sin(beta_3(t))*(1295000.*sin(
beta_3(t))*diff(s_3(t),t)+2135569.*cos(beta_3(t))*diff(beta_3(t),t)+3125000.*
s_3(t)^2*cos(beta_3(t))*diff(beta_3(t),t)+2590000.*s_3(t)*cos(beta_3(t))*diff(
beta_3(t),t)+3125000.*sin(beta_3(t))*s_3(t)*diff(s_3(t),t))*diff(alpha_3(t),t)
, -.3200000000e-6*(-1.+cos(beta_2(t))^2)*(2590000.*s_2(t)+2135569.+3125000.*
s_2(t)^2)*diff(diff(alpha_2(t),t),t)+(sin(beta_2(t))*sin(alpha_2(t))*s_2(t)+.8288000000*sin(beta_2(t))*sin(alpha_2(t)))*lambda4(t)+(-1.*sin(beta_2(t))*cos(

alpha_2(t))*s_2(t)-.8288000000*sin(beta_2(t))*cos(alpha_2(t)))*lambda5(t)+.6400000000e-6*sin(beta_2(t))*(2135569.*cos(beta_2(t))*diff(beta_2(t),t)+

3125000.*sin(beta_2(t))*s_2(t)*diff(s_2(t),t)+3125000.*s_2(t)^2*cos(beta_2(t))
*diff(beta_2(t),t)+2590000.*s_2(t)*cos(beta_2(t))*diff(beta_2(t),t)+1295000.*
sin(beta_2(t))*diff(s_2(t),t))*diff(alpha_2(t),t), -.3200000000e-6*(-1.+cos(
beta_6(t))^2)*(2590000.*s_6(t)+2135569.+3125000.*s_6(t)^2)*diff(diff(alpha_6(t
),t),t)+(sin(beta_6(t))*sin(alpha_6(t))*s_6(t)+.8288000000*sin(beta_6(t))*sin(
alpha_6(t)))*lambda16(t)+(-1.*sin(beta_6(t))*cos(alpha_6(t))*s_6(t)-.8288000000*sin(beta_6(t))*cos(alpha_6(t)))*lambda17(t)+.6400000000e-6*sin(beta_6(t))*(

1295000.*sin(beta_6(t))*diff(s_6(t),t)+3125000.*s_6(t)^2*cos(beta_6(t))*diff(
beta_6(t),t)+3125000.*sin(beta_6(t))*s_6(t)*diff(s_6(t),t)+2590000.*s_6(t)*cos
(beta_6(t))*diff(beta_6(t),t)+2135569.*cos(beta_6(t))*diff(beta_6(t),t))*diff(
alpha_6(t),t), x_P(t)-1.*sin(beta_5(t))*cos(alpha_5(t))*s_5(t)+1.931851653-.8288000000*sin(beta_5(t))*cos(alpha_5(t))-.9659258263*cos(xi_P(t))*cos(eta_P(t))

-.2588190451*sin(xi_P(t))*cos(eta_P(t)), (-5*cos(eta_P(t))^2+10)*diff(diff(
zeta_P(t),t),t)+10*diff(diff(xi_P(t),t),t)*sin(eta_P(t))+(.7071067812*cos(xi_P
(t))*sin(eta_P(t))*cos(zeta_P(t))-.7071067812*sin(xi_P(t))*sin(zeta_P(t))-.7071067812*sin(xi_P(t))*sin(eta_P(t))*cos(zeta_P(t))-.7071067812*cos(xi_P(t))*sin

(zeta_P(t)))*lambda2(t)+(.7071067812*cos(xi_P(t))*sin(eta_P(t))*sin(zeta_P(t))
+.7071067812*sin(xi_P(t))*cos(zeta_P(t))-.7071067812*sin(xi_P(t))*sin(eta_P(t)
)*sin(zeta_P(t))+.7071067812*cos(xi_P(t))*cos(zeta_P(t)))*lambda3(t)+(.9659258263*cos(xi_P(t))*sin(eta_P(t))*cos(zeta_P(t))-.9659258263*sin(xi_P(t))*sin(

zeta_P(t))-.2588190451*sin(xi_P(t))*sin(eta_P(t))*cos(zeta_P(t))-.2588190451*
cos(xi_P(t))*sin(zeta_P(t)))*lambda5(t)+(.9659258263*cos(xi_P(t))*sin(eta_P(t)
)*sin(zeta_P(t))+.9659258263*sin(xi_P(t))*cos(zeta_P(t))-.2588190451*sin(xi_P(
t))*sin(eta_P(t))*sin(zeta_P(t))+.2588190451*cos(xi_P(t))*cos(zeta_P(t)))*
lambda6(t)+(.2588190451*cos(xi_P(t))*sin(eta_P(t))*cos(zeta_P(t))-.2588190451*
sin(xi_P(t))*sin(zeta_P(t))+.9659258263*sin(xi_P(t))*sin(eta_P(t))*cos(zeta_P(
t))+.9659258263*cos(xi_P(t))*sin(zeta_P(t)))*lambda8(t)+(.2588190451*cos(xi_P(
t))*sin(eta_P(t))*sin(zeta_P(t))+.2588190451*sin(xi_P(t))*cos(zeta_P(t))+.9659258263*sin(xi_P(t))*sin(eta_P(t))*sin(zeta_P(t))-.9659258263*cos(xi_P(t))*cos(

zeta_P(t)))*lambda9(t)+(-.2588190451*cos(xi_P(t))*sin(eta_P(t))*cos(zeta_P(t))
+.2588190451*sin(xi_P(t))*sin(zeta_P(t))+.9659258263*sin(xi_P(t))*sin(eta_P(t)
)*cos(zeta_P(t))+.9659258263*cos(xi_P(t))*sin(zeta_P(t)))*lambda11(t)+(-.2588190451*cos(xi_P(t))*sin(eta_P(t))*sin(zeta_P(t))-.2588190451*sin(xi_P(t))*cos(

zeta_P(t))+.9659258263*sin(xi_P(t))*sin(eta_P(t))*sin(zeta_P(t))-.9659258263*
cos(xi_P(t))*cos(zeta_P(t)))*lambda12(t)+(-.9659258263*cos(xi_P(t))*sin(eta_P(
t))*cos(zeta_P(t))+.9659258263*sin(xi_P(t))*sin(zeta_P(t))-.2588190451*sin(
xi_P(t))*sin(eta_P(t))*cos(zeta_P(t))-.2588190451*cos(xi_P(t))*sin(zeta_P(t)))
*lambda14(t)+(-.9659258263*cos(xi_P(t))*sin(eta_P(t))*sin(zeta_P(t))-.9659258263*sin(xi_P(t))*cos(zeta_P(t))-.2588190451*sin(xi_P(t))*sin(eta_P(t))*sin(

zeta_P(t))+.2588190451*cos(xi_P(t))*cos(zeta_P(t)))*lambda15(t)+(-.7071067812*
cos(xi_P(t))*sin(eta_P(t))*cos(zeta_P(t))+.7071067812*sin(xi_P(t))*sin(zeta_P(
t))-.7071067812*sin(xi_P(t))*sin(eta_P(t))*cos(zeta_P(t))-.7071067812*cos(xi_P
(t))*sin(zeta_P(t)))*lambda17(t)+(-.7071067812*cos(xi_P(t))*sin(eta_P(t))*sin(
zeta_P(t))-.7071067812*sin(xi_P(t))*cos(zeta_P(t))-.7071067812*sin(xi_P(t))*
sin(eta_P(t))*sin(zeta_P(t))+.7071067812*cos(xi_P(t))*cos(zeta_P(t)))*lambda18
(t)+10*cos(eta_P(t))*(sin(eta_P(t))*diff(zeta_P(t),t)+diff(xi_P(t),t))*diff(
eta_P(t),t), z_P(t)-1.*cos(beta_5(t))*s_5(t)-.8288000000*cos(beta_5(t))+.9659258263*cos(xi_P(t))*sin(eta_P(t))*cos(zeta_P(t))-.9659258263*sin(xi_P(t))*sin(

zeta_P(t))+.2588190451*sin(xi_P(t))*sin(eta_P(t))*cos(zeta_P(t))+.2588190451*
cos(xi_P(t))*sin(zeta_P(t)), 5*diff(diff(eta_P(t),t),t)+(-.7071067812*cos(xi_P
(t))*sin(eta_P(t))+.7071067812*sin(xi_P(t))*sin(eta_P(t)))*lambda1(t)+(.7071067812*sin(zeta_P(t))*cos(xi_P(t))*cos(eta_P(t))-.7071067812*sin(zeta_P(t))*sin(

xi_P(t))*cos(eta_P(t)))*lambda2(t)+(-.7071067812*cos(zeta_P(t))*cos(xi_P(t))*
cos(eta_P(t))+.7071067812*cos(zeta_P(t))*sin(xi_P(t))*cos(eta_P(t)))*lambda3(t
)+(-.9659258263*cos(xi_P(t))*sin(eta_P(t))+.2588190451*sin(xi_P(t))*sin(eta_P(
t)))*lambda4(t)+(.9659258263*sin(zeta_P(t))*cos(xi_P(t))*cos(eta_P(t))-.2588190451*sin(zeta_P(t))*sin(xi_P(t))*cos(eta_P(t)))*lambda5(t)+(-.9659258263*cos(

zeta_P(t))*cos(xi_P(t))*cos(eta_P(t))+.2588190451*cos(zeta_P(t))*sin(xi_P(t))*
cos(eta_P(t)))*lambda6(t)+(-.2588190451*cos(xi_P(t))*sin(eta_P(t))-.9659258263
*sin(xi_P(t))*sin(eta_P(t)))*lambda7(t)+(.2588190451*sin(zeta_P(t))*cos(xi_P(t
))*cos(eta_P(t))+.9659258263*sin(zeta_P(t))*sin(xi_P(t))*cos(eta_P(t)))*
lambda8(t)+(-.2588190451*cos(zeta_P(t))*cos(xi_P(t))*cos(eta_P(t))-.9659258263
*cos(zeta_P(t))*sin(xi_P(t))*cos(eta_P(t)))*lambda9(t)+(.2588190451*cos(xi_P(t
))*sin(eta_P(t))-.9659258263*sin(xi_P(t))*sin(eta_P(t)))*lambda10(t)+(-.2588190451*sin(zeta_P(t))*cos(xi_P(t))*cos(eta_P(t))+.9659258263*sin(zeta_P(t))*sin(

xi_P(t))*cos(eta_P(t)))*lambda11(t)+(.2588190451*cos(zeta_P(t))*cos(xi_P(t))*
cos(eta_P(t))-.9659258263*cos(zeta_P(t))*sin(xi_P(t))*cos(eta_P(t)))*lambda12(
t)+(.9659258263*cos(xi_P(t))*sin(eta_P(t))+.2588190451*sin(xi_P(t))*sin(eta_P(
t)))*lambda13(t)+(-.9659258263*sin(zeta_P(t))*cos(xi_P(t))*cos(eta_P(t))-.2588190451*sin(zeta_P(t))*sin(xi_P(t))*cos(eta_P(t)))*lambda14(t)+(.9659258263*cos

(zeta_P(t))*cos(xi_P(t))*cos(eta_P(t))+.2588190451*cos(zeta_P(t))*sin(xi_P(t))
*cos(eta_P(t)))*lambda15(t)+(.7071067812*cos(xi_P(t))*sin(eta_P(t))+.7071067812*sin(xi_P(t))*sin(eta_P(t)))*lambda16(t)+(-.7071067812*sin(zeta_P(t))*cos(

xi_P(t))*cos(eta_P(t))-.7071067812*sin(zeta_P(t))*sin(xi_P(t))*cos(eta_P(t)))*
lambda17(t)+(.7071067812*cos(zeta_P(t))*cos(xi_P(t))*cos(eta_P(t))+.7071067812
*cos(zeta_P(t))*sin(xi_P(t))*cos(eta_P(t)))*lambda18(t)-5*cos(eta_P(t))*diff(
zeta_P(t),t)*(2*diff(xi_P(t),t)+sin(eta_P(t))*diff(zeta_P(t),t)), x_P(t)-1.*
sin(beta_6(t))*cos(alpha_6(t))*s_6(t)-.8288000000*sin(beta_6(t))*cos(alpha_6(t
))+.5176380902-.7071067812*cos(xi_P(t))*cos(eta_P(t))-.7071067812*sin(xi_P(t))
*cos(eta_P(t)), y_P(t)-1.*sin(beta_3(t))*sin(alpha_3(t))*s_3(t)+1.414213562+.2588190451*cos(xi_P(t))*sin(eta_P(t))*sin(zeta_P(t))+.2588190451*sin(xi_P(t))*

cos(zeta_P(t))+.9659258263*sin(xi_P(t))*sin(eta_P(t))*sin(zeta_P(t))-.9659258263*cos(xi_P(t))*cos(zeta_P(t))-.8288000000*sin(beta_3(t))*sin(alpha_3(t)), y_P

(t)-1.*sin(beta_6(t))*sin(alpha_6(t))*s_6(t)-.8288000000*sin(beta_6(t))*sin(
alpha_6(t))-1.931851653-.7071067812*cos(xi_P(t))*sin(eta_P(t))*sin(zeta_P(t))-.7071067812*sin(xi_P(t))*cos(zeta_P(t))-.7071067812*sin(xi_P(t))*sin(eta_P(t))

*sin(zeta_P(t))+.7071067812*cos(xi_P(t))*cos(zeta_P(t)), (.6833820800+s_1(t)^2
+.8288000000*s_1(t))*diff(diff(beta_1(t),t),t)+(-1.*cos(alpha_1(t))*cos(beta_1
(t))*s_1(t)-.8288000000*cos(beta_1(t))*cos(alpha_1(t)))*lambda1(t)+(-1.*sin(
alpha_1(t))*cos(beta_1(t))*s_1(t)-.8288000000*cos(beta_1(t))*sin(alpha_1(t)))*
lambda2(t)+(sin(beta_1(t))*s_1(t)+.8288000000*sin(beta_1(t)))*lambda3(t)+2.*
s_1(t)*diff(beta_1(t),t)*diff(s_1(t),t)-12.19579200*sin(beta_1(t))-.6833820800
*cos(beta_1(t))*sin(beta_1(t))*diff(alpha_1(t),t)^2-9.810000000*sin(beta_1(t))
*s_1(t)+.8288000000*diff(beta_1(t),t)*diff(s_1(t),t)-.8288000000*cos(beta_1(t)
)*s_1(t)*diff(alpha_1(t),t)^2*sin(beta_1(t))-cos(beta_1(t))*s_1(t)^2*diff(
alpha_1(t),t)^2*sin(beta_1(t)), 10*sin(eta_P(t))*diff(diff(zeta_P(t),t),t)+10*
diff(diff(xi_P(t),t),t)+(-.7071067812*sin(xi_P(t))*cos(eta_P(t))-.7071067812*
cos(xi_P(t))*cos(eta_P(t)))*lambda1(t)+(-.7071067812*cos(xi_P(t))*sin(eta_P(t)
)*sin(zeta_P(t))-.7071067812*sin(xi_P(t))*cos(zeta_P(t))-.7071067812*sin(xi_P(
t))*sin(eta_P(t))*sin(zeta_P(t))+.7071067812*cos(xi_P(t))*cos(zeta_P(t)))*
lambda2(t)+(.7071067812*sin(xi_P(t))*sin(eta_P(t))*cos(zeta_P(t))+.7071067812*
cos(xi_P(t))*sin(zeta_P(t))+.7071067812*cos(xi_P(t))*sin(eta_P(t))*cos(zeta_P(
t))-.7071067812*sin(xi_P(t))*sin(zeta_P(t)))*lambda3(t)+(-.9659258263*sin(xi_P
(t))*cos(eta_P(t))-.2588190451*cos(xi_P(t))*cos(eta_P(t)))*lambda4(t)+(.9659258263*cos(xi_P(t))*cos(zeta_P(t))-.9659258263*sin(xi_P(t))*sin(eta_P(t))*sin(

zeta_P(t))-.2588190451*sin(xi_P(t))*cos(zeta_P(t))-.2588190451*cos(xi_P(t))*
sin(eta_P(t))*sin(zeta_P(t)))*lambda5(t)+(.2588190451*cos(xi_P(t))*sin(eta_P(t
))*cos(zeta_P(t))-.2588190451*sin(xi_P(t))*sin(zeta_P(t))+.9659258263*sin(xi_P
(t))*sin(eta_P(t))*cos(zeta_P(t))+.9659258263*cos(xi_P(t))*sin(zeta_P(t)))*
lambda6(t)+(-.2588190451*sin(xi_P(t))*cos(eta_P(t))+.9659258263*cos(xi_P(t))*
cos(eta_P(t)))*lambda7(t)+(.9659258263*cos(xi_P(t))*sin(eta_P(t))*sin(zeta_P(t
))+.9659258263*sin(xi_P(t))*cos(zeta_P(t))-.2588190451*sin(xi_P(t))*sin(eta_P(
t))*sin(zeta_P(t))+.2588190451*cos(xi_P(t))*cos(zeta_P(t)))*lambda8(t)+(.2588190451*sin(xi_P(t))*sin(eta_P(t))*cos(zeta_P(t))+.2588190451*cos(xi_P(t))*sin(

zeta_P(t))-.9659258263*cos(xi_P(t))*sin(eta_P(t))*cos(zeta_P(t))+.9659258263*
sin(xi_P(t))*sin(zeta_P(t)))*lambda9(t)+(.2588190451*sin(xi_P(t))*cos(eta_P(t)
)+.9659258263*cos(xi_P(t))*cos(eta_P(t)))*lambda10(t)+(-.2588190451*cos(xi_P(t
))*cos(zeta_P(t))+.2588190451*sin(xi_P(t))*sin(eta_P(t))*sin(zeta_P(t))+.9659258263*sin(xi_P(t))*cos(zeta_P(t))+.9659258263*cos(xi_P(t))*sin(eta_P(t))*sin(

zeta_P(t)))*lambda11(t)+(-.9659258263*cos(xi_P(t))*sin(eta_P(t))*cos(zeta_P(t)
)+.9659258263*sin(xi_P(t))*sin(zeta_P(t))-.2588190451*sin(xi_P(t))*sin(eta_P(t
))*cos(zeta_P(t))-.2588190451*cos(xi_P(t))*sin(zeta_P(t)))*lambda12(t)+(.9659258263*sin(xi_P(t))*cos(eta_P(t))-.2588190451*cos(xi_P(t))*cos(eta_P(t)))*

lambda13(t)+(-.2588190451*cos(xi_P(t))*sin(eta_P(t))*sin(zeta_P(t))-.2588190451*sin(xi_P(t))*cos(zeta_P(t))+.9659258263*sin(xi_P(t))*sin(eta_P(t))*sin(

zeta_P(t))-.9659258263*cos(xi_P(t))*cos(zeta_P(t)))*lambda14(t)+(-.9659258263*
sin(xi_P(t))*sin(eta_P(t))*cos(zeta_P(t))-.9659258263*cos(xi_P(t))*sin(zeta_P(
t))+.2588190451*cos(xi_P(t))*sin(eta_P(t))*cos(zeta_P(t))-.2588190451*sin(xi_P
(t))*sin(zeta_P(t)))*lambda15(t)+(.7071067812*sin(xi_P(t))*cos(eta_P(t))-.7071067812*cos(xi_P(t))*cos(eta_P(t)))*lambda16(t)+(-.7071067812*cos(xi_P(t))*cos(

zeta_P(t))+.7071067812*sin(xi_P(t))*sin(eta_P(t))*sin(zeta_P(t))-.7071067812*
sin(xi_P(t))*cos(zeta_P(t))-.7071067812*cos(xi_P(t))*sin(eta_P(t))*sin(zeta_P(
t)))*lambda17(t)+(.7071067812*cos(xi_P(t))*sin(eta_P(t))*cos(zeta_P(t))-.7071067812*sin(xi_P(t))*sin(zeta_P(t))-.7071067812*sin(xi_P(t))*sin(eta_P(t))*cos(

zeta_P(t))-.7071067812*cos(xi_P(t))*sin(zeta_P(t)))*lambda18(t)+10*cos(eta_P(t
))*diff(zeta_P(t),t)*diff(eta_P(t),t), 10*diff(diff(x_P(t),t),t)+lambda1(t)+
lambda4(t)+lambda7(t)+lambda10(t)+lambda13(t)+lambda16(t), 10*diff(diff(y_P(t)
,t),t)+lambda2(t)+lambda5(t)+lambda8(t)+lambda11(t)+lambda14(t)+lambda17(t),
x_P(t)-1.*sin(beta_1(t))*cos(alpha_1(t))*s_1(t)-.5176380902+.7071067812*cos(
xi_P(t))*cos(eta_P(t))-.7071067812*sin(xi_P(t))*cos(eta_P(t))-.8288000000*sin(
beta_1(t))*cos(alpha_1(t)), (.8288000000*s_2(t)+.6833820800+s_2(t)^2)*diff(
diff(beta_2(t),t),t)+(-1.*cos(alpha_2(t))*cos(beta_2(t))*s_2(t)-.8288000000*
cos(beta_2(t))*cos(alpha_2(t)))*lambda4(t)+(-1.*sin(alpha_2(t))*cos(beta_2(t))
*s_2(t)-.8288000000*cos(beta_2(t))*sin(alpha_2(t)))*lambda5(t)+(sin(beta_2(t))
*s_2(t)+.8288000000*sin(beta_2(t)))*lambda6(t)-.6833820800*cos(beta_2(t))*sin(
beta_2(t))*diff(alpha_2(t),t)^2-9.810000000*sin(beta_2(t))*s_2(t)-cos(beta_2(t
))*s_2(t)^2*diff(alpha_2(t),t)^2*sin(beta_2(t))-.8288000000*cos(beta_2(t))*s_2
(t)*diff(alpha_2(t),t)^2*sin(beta_2(t))+2.*s_2(t)*diff(beta_2(t),t)*diff(s_2(t
),t)+.8288000000*diff(beta_2(t),t)*diff(s_2(t),t)-12.19579200*sin(beta_2(t)),
z_P(t)-1.*cos(beta_1(t))*s_1(t)-.7071067812*cos(xi_P(t))*sin(eta_P(t))*cos(
zeta_P(t))+.7071067812*sin(xi_P(t))*sin(zeta_P(t))+.7071067812*sin(xi_P(t))*
sin(eta_P(t))*cos(zeta_P(t))+.7071067812*cos(xi_P(t))*sin(zeta_P(t))-.8288000000*cos(beta_1(t)), y_P(t)-1.*sin(beta_1(t))*sin(alpha_1(t))*s_1(t)-1.931851653

+.7071067812*cos(xi_P(t))*sin(eta_P(t))*sin(zeta_P(t))+.7071067812*sin(xi_P(t)
)*cos(zeta_P(t))-.7071067812*sin(xi_P(t))*sin(eta_P(t))*sin(zeta_P(t))+.7071067812*cos(xi_P(t))*cos(zeta_P(t))-.8288000000*sin(beta_1(t))*sin(alpha_1(t)), (

.6833820800+.8288000000*s_3(t)+s_3(t)^2)*diff(diff(beta_3(t),t),t)+(-1.*cos(
alpha_3(t))*cos(beta_3(t))*s_3(t)-.8288000000*cos(beta_3(t))*cos(alpha_3(t)))*
lambda7(t)+(-1.*sin(alpha_3(t))*cos(beta_3(t))*s_3(t)-.8288000000*cos(beta_3(t
))*sin(alpha_3(t)))*lambda8(t)+(sin(beta_3(t))*s_3(t)+.8288000000*sin(beta_3(t
)))*lambda9(t)+.8288000000*diff(beta_3(t),t)*diff(s_3(t),t)+2.*s_3(t)*diff(
beta_3(t),t)*diff(s_3(t),t)-9.810000000*sin(beta_3(t))*s_3(t)-.8288000000*cos(
beta_3(t))*s_3(t)*diff(alpha_3(t),t)^2*sin(beta_3(t))-cos(beta_3(t))*s_3(t)^2*
diff(alpha_3(t),t)^2*sin(beta_3(t))-.6833820800*cos(beta_3(t))*sin(beta_3(t))*
diff(alpha_3(t),t)^2-12.19579200*sin(beta_3(t)), (.6833820800+.8288000000*s_4(
t)+s_4(t)^2)*diff(diff(beta_4(t),t),t)+(-1.*cos(alpha_4(t))*cos(beta_4(t))*s_4
(t)-.8288000000*cos(beta_4(t))*cos(alpha_4(t)))*lambda10(t)+(-1.*sin(alpha_4(t
))*cos(beta_4(t))*s_4(t)-.8288000000*cos(beta_4(t))*sin(alpha_4(t)))*lambda11(
t)+(sin(beta_4(t))*s_4(t)+.8288000000*sin(beta_4(t)))*lambda12(t)+.8288000000*
diff(beta_4(t),t)*diff(s_4(t),t)-9.810000000*sin(beta_4(t))*s_4(t)-12.19579200
*sin(beta_4(t))+2.*s_4(t)*diff(beta_4(t),t)*diff(s_4(t),t)-cos(beta_4(t))*s_4(
t)^2*diff(alpha_4(t),t)^2*sin(beta_4(t))-.8288000000*cos(beta_4(t))*s_4(t)*
diff(alpha_4(t),t)^2*sin(beta_4(t))-.6833820800*cos(beta_4(t))*sin(beta_4(t))*
diff(alpha_4(t),t)^2, (.8288000000*s_5(t)+.6833820800+s_5(t)^2)*diff(diff(
beta_5(t),t),t)+(-1.*cos(alpha_5(t))*cos(beta_5(t))*s_5(t)-.8288000000*cos(
beta_5(t))*cos(alpha_5(t)))*lambda13(t)+(-1.*sin(alpha_5(t))*cos(beta_5(t))*
s_5(t)-.8288000000*cos(beta_5(t))*sin(alpha_5(t)))*lambda14(t)+(sin(beta_5(t))
*s_5(t)+.8288000000*sin(beta_5(t)))*lambda15(t)-9.810000000*sin(beta_5(t))*s_5
(t)-12.19579200*sin(beta_5(t))+.8288000000*diff(beta_5(t),t)*diff(s_5(t),t)-.6833820800*cos(beta_5(t))*sin(beta_5(t))*diff(alpha_5(t),t)^2-cos(beta_5(t))*

s_5(t)^2*diff(alpha_5(t),t)^2*sin(beta_5(t))-.8288000000*cos(beta_5(t))*s_5(t)
*diff(alpha_5(t),t)^2*sin(beta_5(t))+2.*s_5(t)*diff(beta_5(t),t)*diff(s_5(t),t
), z_P(t)-1.*cos(beta_6(t))*s_6(t)-.8288000000*cos(beta_6(t))+.7071067812*cos(
xi_P(t))*sin(eta_P(t))*cos(zeta_P(t))-.7071067812*sin(xi_P(t))*sin(zeta_P(t))+
.7071067812*sin(xi_P(t))*sin(eta_P(t))*cos(zeta_P(t))+.7071067812*cos(xi_P(t))
*sin(zeta_P(t)), y_P(t)-1.*sin(beta_2(t))*sin(alpha_2(t))*s_2(t)+.5176380902+.9659258263*cos(xi_P(t))*sin(eta_P(t))*sin(zeta_P(t))+.9659258263*sin(xi_P(t))*

cos(zeta_P(t))-.2588190451*sin(xi_P(t))*sin(eta_P(t))*sin(zeta_P(t))+.2588190451*cos(xi_P(t))*cos(zeta_P(t))-.8288000000*sin(beta_2(t))*sin(alpha_2(t)), y_P

(t)-1.*sin(beta_5(t))*sin(alpha_5(t))*s_5(t)+.5176380902-.8288000000*sin(
beta_5(t))*sin(alpha_5(t))-.9659258263*cos(xi_P(t))*sin(eta_P(t))*sin(zeta_P(t
))-.9659258263*sin(xi_P(t))*cos(zeta_P(t))-.2588190451*sin(xi_P(t))*sin(eta_P(
t))*sin(zeta_P(t))+.2588190451*cos(xi_P(t))*cos(zeta_P(t)), (.8288000000*s_6(t
)+.6833820800+s_6(t)^2)*diff(diff(beta_6(t),t),t)+(-1.*cos(alpha_6(t))*cos(
beta_6(t))*s_6(t)-.8288000000*cos(beta_6(t))*cos(alpha_6(t)))*lambda16(t)+(-1.
*sin(alpha_6(t))*cos(beta_6(t))*s_6(t)-.8288000000*cos(beta_6(t))*sin(alpha_6(
t)))*lambda17(t)+(sin(beta_6(t))*s_6(t)+.8288000000*sin(beta_6(t)))*lambda18(t
)-12.19579200*sin(beta_6(t))+.8288000000*diff(beta_6(t),t)*diff(s_6(t),t)-9.810000000*sin(beta_6(t))*s_6(t)-.6833820800*cos(beta_6(t))*sin(beta_6(t))*diff(

alpha_6(t),t)^2-.8288000000*cos(beta_6(t))*s_6(t)*diff(alpha_6(t),t)^2*sin(
beta_6(t))-cos(beta_6(t))*s_6(t)^2*diff(alpha_6(t),t)^2*sin(beta_6(t))+2.*s_6(
t)*diff(beta_6(t),t)*diff(s_6(t),t), x_P(t)-1.*sin(beta_3(t))*cos(alpha_3(t))*
s_3(t)-1.414213562+.2588190451*cos(xi_P(t))*cos(eta_P(t))+.9659258263*sin(xi_P
(t))*cos(eta_P(t))-.8288000000*sin(beta_3(t))*cos(alpha_3(t)), x_P(t)-1.*sin(
beta_4(t))*cos(alpha_4(t))*s_4(t)-.8288000000*sin(beta_4(t))*cos(alpha_4(t))+1.414213562-.2588190451*cos(xi_P(t))*cos(eta_P(t))+.9659258263*sin(xi_P(t))*cos

(eta_P(t)), y_P(t)-1.*sin(beta_4(t))*sin(alpha_4(t))*s_4(t)-.8288000000*sin(
beta_4(t))*sin(alpha_4(t))+1.414213562-.2588190451*cos(xi_P(t))*sin(eta_P(t))*
sin(zeta_P(t))-.2588190451*sin(xi_P(t))*cos(zeta_P(t))+.9659258263*sin(xi_P(t)
)*sin(eta_P(t))*sin(zeta_P(t))-.9659258263*cos(xi_P(t))*cos(zeta_P(t)), z_P(t)
-1.*cos(beta_4(t))*s_4(t)-.8288000000*cos(beta_4(t))+.2588190451*cos(xi_P(t))*
sin(eta_P(t))*cos(zeta_P(t))-.2588190451*sin(xi_P(t))*sin(zeta_P(t))-.9659258263*sin(xi_P(t))*sin(eta_P(t))*cos(zeta_P(t))-.9659258263*cos(xi_P(t))*sin(

zeta_P(t)), z_P(t)-1.*cos(beta_3(t))*s_3(t)-.2588190451*cos(xi_P(t))*sin(eta_P
(t))*cos(zeta_P(t))+.2588190451*sin(xi_P(t))*sin(zeta_P(t))-.9659258263*sin(
xi_P(t))*sin(eta_P(t))*cos(zeta_P(t))-.9659258263*cos(xi_P(t))*sin(zeta_P(t))-.8288000000*cos(beta_3(t)), diff(diff(s_6(t),t),t)-1.*sin(beta_6(t))*cos(

alpha_6(t))*lambda16(t)-1.*sin(beta_6(t))*sin(alpha_6(t))*lambda17(t)-1.*cos(
beta_6(t))*lambda18(t)-28-diff(beta_6(t),t)^2*s_6(t)+1.*diff(alpha_6(t),t)^2*
s_6(t)*cos(beta_6(t))^2+.4144000000*diff(alpha_6(t),t)^2*cos(beta_6(t))^2-diff
(alpha_6(t),t)^2*s_6(t)-.4144000000*diff(beta_6(t),t)^2-.4144000000*diff(
alpha_6(t),t)^2+9.810000000*cos(beta_6(t)), diff(diff(s_5(t),t),t)-1.*sin(
beta_5(t))*cos(alpha_5(t))*lambda13(t)-1.*sin(beta_5(t))*sin(alpha_5(t))*
lambda14(t)-1.*cos(beta_5(t))*lambda15(t)+1.*diff(alpha_5(t),t)^2*s_5(t)*cos(
beta_5(t))^2+9.810000000*cos(beta_5(t))-28-.4144000000*diff(beta_5(t),t)^2-
diff(beta_5(t),t)^2*s_5(t)-diff(alpha_5(t),t)^2*s_5(t)+.4144000000*diff(
alpha_5(t),t)^2*cos(beta_5(t))^2-.4144000000*diff(alpha_5(t),t)^2, diff(diff(
s_4(t),t),t)-1.*sin(beta_4(t))*cos(alpha_4(t))*lambda10(t)-1.*sin(beta_4(t))*
sin(alpha_4(t))*lambda11(t)-1.*cos(beta_4(t))*lambda12(t)+.4144000000*diff(
alpha_4(t),t)^2*cos(beta_4(t))^2-.4144000000*diff(alpha_4(t),t)^2+9.810000000*
cos(beta_4(t))-28-3.0*sin(4*Pi*t)-diff(alpha_4(t),t)^2*s_4(t)+1.*diff(alpha_4(
t),t)^2*s_4(t)*cos(beta_4(t))^2-diff(beta_4(t),t)^2*s_4(t)-.4144000000*diff(
beta_4(t),t)^2, diff(diff(s_3(t),t),t)-1.*sin(beta_3(t))*cos(alpha_3(t))*
lambda7(t)-1.*sin(beta_3(t))*sin(alpha_3(t))*lambda8(t)-1.*cos(beta_3(t))*
lambda9(t)+9.810000000*cos(beta_3(t))-.4144000000*diff(alpha_3(t),t)^2-28-3.0*
sin(4*Pi*t)-.4144000000*diff(beta_3(t),t)^2-diff(beta_3(t),t)^2*s_3(t)-diff(
alpha_3(t),t)^2*s_3(t)+1.*diff(alpha_3(t),t)^2*s_3(t)*cos(beta_3(t))^2+.4144000000*diff(alpha_3(t),t)^2*cos(beta_3(t))^2, diff(diff(s_2(t),t),t)-1.*sin(

beta_2(t))*cos(alpha_2(t))*lambda4(t)-1.*sin(beta_2(t))*sin(alpha_2(t))*
lambda5(t)-1.*cos(beta_2(t))*lambda6(t)-28-5.0*cos(2*Pi*t)+1.*diff(alpha_2(t),
t)^2*s_2(t)*cos(beta_2(t))^2-diff(beta_2(t),t)^2*s_2(t)-diff(alpha_2(t),t)^2*
s_2(t)+9.810000000*cos(beta_2(t))-.4144000000*diff(alpha_2(t),t)^2-.4144000000
*diff(beta_2(t),t)^2+.4144000000*diff(alpha_2(t),t)^2*cos(beta_2(t))^2, diff(
diff(s_1(t),t),t)-1.*sin(beta_1(t))*cos(alpha_1(t))*lambda1(t)-1.*sin(beta_1(t
))*sin(alpha_1(t))*lambda2(t)-1.*cos(beta_1(t))*lambda3(t)+1.*diff(alpha_1(t),
t)^2*s_1(t)*cos(beta_1(t))^2+.4144000000*diff(alpha_1(t),t)^2*cos(beta_1(t))^2
-diff(alpha_1(t),t)^2*s_1(t)-.4144000000*diff(alpha_1(t),t)^2-diff(beta_1(t),t
)^2*s_1(t)-.4144000000*diff(beta_1(t),t)^2+9.810000000*cos(beta_1(t))-28-5.0*
cos(2*Pi*t), 10*diff(diff(z_P(t),t),t)+lambda3(t)+lambda6(t)+lambda9(t)+
lambda12(t)+lambda15(t)+lambda18(t)+98.10000000}:

stICs := {D(alpha_5)(0) = 0, D(beta_5)(0) = 0, D(alpha_6)(0) = 0, D(beta_2)(0)
= 0, D(s_5)(0) = 0, D(beta_6)(0) = 0, D(s_6)(0) = 0, D(beta_3)(0) = 0, D(
zeta_P)(0) = 0, D(alpha_2)(0) = 0, D(eta_P)(0) = 0, D(x_P)(0) = 0, D(y_P)(0) =
0, D(alpha_3)(0) = 0, D(xi_P)(0) = 0, D(z_P)(0) = 0, x_P(0) = 0, y_P(0) = 0,
z_P(0) = 3, s_1(0) = 2.417104864, s_2(0) = 2.417104864, s_3(0) = 2.417104864,
s_4(0) = 2.417104864, s_5(0) = 2.417104864, s_6(0) = 2.417104864, D(alpha_4)(0
) = 0, D(s_1)(0) = 0, D(alpha_1)(0) = 0, D(s_2)(0) = 0, zeta_P(0) = 0, eta_P(0
) = 0, xi_P(0) = 0, alpha_1(0) = -1.417312476, beta_1(0) = .3917521166,
alpha_2(0) = 2.464510027, beta_2(0) = .3917521166, alpha_3(0) = 2.771477729,
beta_3(0) = .3917521166, alpha_4(0) = .3701149244, beta_4(0) = .3917521166,
alpha_5(0) = .6770826264, beta_5(0) = .3917521166, alpha_6(0) = -1.724280178,
beta_6(0) = .3917521166, D(s_3)(0) = 0, D(beta_4)(0) = 0, D(beta_1)(0) = 0, D(
s_4)(0) = 0}:

  

In Maple 10, this system could not be brought into the required DAE form and numerically solved in anything under an hour. Now however:

tt := time():

dsn := dsolve(stEQs union stICs, numeric, implicit=true):

dsn(1);

t=1.&comma;alpha_1t=−1.41404214314427&comma;&DifferentialD;&DifferentialD;talpha_1t=0.00849473764599224&comma;alpha_2t=2.46031987052574&comma;&DifferentialD;&DifferentialD;talpha_2t=−0.00427806160361022&comma;alpha_3t=2.76667851916170&comma;&DifferentialD;&DifferentialD;talpha_3t=−0.00224383300346078&comma;alpha_4t=0.373189183301449&comma;&DifferentialD;&DifferentialD;talpha_4t=−0.00345824250260502&comma;alpha_5t=0.676103669084111&comma;&DifferentialD;&DifferentialD;talpha_5t=−0.0126138578651418&comma;alpha_6t=−1.72061386826404&comma;&DifferentialD;&DifferentialD;talpha_6t=0.0145971225866554&comma;beta_1t=0.399270368285695&comma;&DifferentialD;&DifferentialD;tbeta_1t=0.0168856342949760&comma;beta_2t=0.395754725106839&comma;&DifferentialD;&DifferentialD;tbeta_2t=0.00854132788090072&comma;beta_3t=0.390747038446706&comma;&DifferentialD;&DifferentialD;tbeta_3t=0.0115082928413281&comma;beta_4t=0.394386856258629&comma;&DifferentialD;&DifferentialD;tbeta_4t=0.0218166276873265&comma;beta_5t=0.407005428357208&comma;&DifferentialD;&DifferentialD;tbeta_5t=0.0389816800352700&comma;beta_6t=0.406830807784426&comma;&DifferentialD;&DifferentialD;tbeta_6t=0.0366850260352228&comma;eta_Pt=−0.0422894233981797&comma;&DifferentialD;&DifferentialD;teta_Pt=−0.105787671732061&comma;&lambda;1t=−2.10919513521311&comma;&lambda;10t=−10.4984524640529&comma;&lambda;11t=−4.15923393631383&comma;&lambda;12t=−16.1304320809642&comma;&lambda;13t=−9.02987945391216&comma;&lambda;14t=−7.31967620701184&comma;&lambda;15t=−16.2842801493366&comma;&lambda;16t=1.74150222732886&comma;&lambda;17t=11.3835236220591&comma;&lambda;18t=−15.8371858927843&comma;&lambda;2t=13.0088590989426&comma;&lambda;3t=−19.0184790709148&comma;&lambda;4t=10.0748474937966&comma;&lambda;5t=−8.24678875909166&comma;&lambda;6t=−18.9498152109804&comma;&lambda;7t=10.2437486523454&comma;&lambda;8t=−4.11354762594244&comma;&lambda;9t=−15.6064276928410&comma;s_1t=2.35208666252732&comma;&DifferentialD;&DifferentialD;ts_1t=−0.118897311472506&comma;s_2t=2.38757724404409&comma;&DifferentialD;&DifferentialD;ts_2t=−0.0809563086256958&comma;s_3t=2.42579134679721&comma;&DifferentialD;&DifferentialD;ts_3t=−0.105119915874658&comma;s_4t=2.40699590616213&comma;&DifferentialD;&DifferentialD;ts_4t=−0.150615237136879&comma;s_5t=2.31403328393115&comma;&DifferentialD;&DifferentialD;ts_5t=−0.262431777006591&comma;s_6t=2.29725446223393&comma;&DifferentialD;&DifferentialD;ts_6t=−0.255377349604363&comma;x_Pt=0.00332650598171118&comma;&DifferentialD;&DifferentialD;tx_Pt=0.0111590920256401&comma;xi_Pt=−0.00124602843239396&comma;&DifferentialD;&DifferentialD;txi_Pt=−0.00409946001158069&comma;y_Pt=0.00368685194005796&comma;&DifferentialD;&DifferentialD;ty_Pt=−0.00260844663960202&comma;z_Pt=2.94202195913911&comma;&DifferentialD;&DifferentialD;tz_Pt=−0.177126550474447&comma;zeta_Pt=−0.0583839149780300&comma;&DifferentialD;&DifferentialD;tzeta_Pt=−0.0400305840410983

(49)

time()-tt;

3.827

(50)

DAE Efficiency

  

This example demonstrates that the DAE engine has also been made more efficient. This stiff DAE system is for a 10 weight multiple pendulum problem, where xi,yi represent the position of the ith weight as a function of time.

  

The DAE system is:

dsys := {diff(diff(x[1](t),t),t)+2*lambda[1](t)*x[1](t)+2*lambda[2](t)*(-x[2](t)+x[1](
t)), diff(diff(y[1](t),t),t)+49/5+2*lambda[1](t)*y[1](t)+2*lambda[2](t)*(-y[2]
(t)+y[1](t)), diff(diff(x[2](t),t),t)-2*lambda[2](t)*(-x[2](t)+x[1](t))+2*
lambda[3](t)*(-x[3](t)+x[2](t)), diff(diff(y[2](t),t),t)+49/5-2*lambda[2](t)*(
-y[2](t)+y[1](t))+2*lambda[3](t)*(-y[3](t)+y[2](t)), diff(diff(x[3](t),t),t)-2
*lambda[3](t)*(-x[3](t)+x[2](t))+2*lambda[4](t)*(-x[4](t)+x[3](t)), diff(diff(
y[3](t),t),t)+49/5-2*lambda[3](t)*(-y[3](t)+y[2](t))+2*lambda[4](t)*(-y[4](t)+
y[3](t)), diff(diff(x[4](t),t),t)-2*lambda[4](t)*(-x[4](t)+x[3](t))+2*lambda[5
](t)*(-x[5](t)+x[4](t)), diff(diff(y[4](t),t),t)+49/5-2*lambda[4](t)*(-y[4](t)
+y[3](t))+2*lambda[5](t)*(-y[5](t)+y[4](t)), diff(diff(x[5](t),t),t)-2*lambda[
5](t)*(-x[5](t)+x[4](t))+2*lambda[6](t)*(-x[6](t)+x[5](t)), diff(diff(y[5](t),
t),t)+49/5-2*lambda[5](t)*(-y[5](t)+y[4](t))+2*lambda[6](t)*(-y[6](t)+y[5](t))
, diff(diff(x[6](t),t),t)-2*lambda[6](t)*(-x[6](t)+x[5](t))+2*lambda[7](t)*(-x
[7](t)+x[6](t)), diff(diff(y[6](t),t),t)+49/5-2*lambda[6](t)*(-y[6](t)+y[5](t)
)+2*lambda[7](t)*(-y[7](t)+y[6](t)), diff(diff(x[7](t),t),t)+2*lambda[7](t)*(x
[7](t)-x[6](t))+2*lambda[8](t)*(-x[8](t)+x[7](t)), diff(diff(y[7](t),t),t)+49/
5+2*lambda[7](t)*(y[7](t)-y[6](t))+2*lambda[8](t)*(-y[8](t)+y[7](t)), diff(
diff(x[8](t),t),t)-2*lambda[8](t)*(-x[8](t)+x[7](t))+2*lambda[9](t)*(-x[9](t)+
x[8](t)), diff(diff(y[8](t),t),t)+49/5-2*lambda[8](t)*(-y[8](t)+y[7](t))+2*
lambda[9](t)*(-y[9](t)+y[8](t)), diff(diff(x[9](t),t),t)-2*lambda[9](t)*(-x[9]
(t)+x[8](t))+2*lambda[10](t)*(-x[10](t)+x[9](t)), diff(diff(y[9](t),t),t)+49/5
-2*lambda[9](t)*(-y[9](t)+y[8](t))+2*lambda[10](t)*(-y[10](t)+y[9](t)), diff(
diff(x[10](t),t),t)-2*lambda[10](t)*(-x[10](t)+x[9](t)), diff(diff(y[10](t),t)
,t)+49/5-2*lambda[10](t)*(-y[10](t)+y[9](t)), x[1](t)^2+y[1](t)^2-1, (x[2](t)-
x[1](t))^2+(y[2](t)-y[1](t))^2-1, (x[3](t)-x[2](t))^2+(y[3](t)-y[2](t))^2-1, (
x[4](t)-x[3](t))^2+(y[4](t)-y[3](t))^2-1, (x[5](t)-x[4](t))^2+(y[5](t)-y[4](t)
)^2-1, (x[6](t)-x[5](t))^2+(y[6](t)-y[5](t))^2-1, (x[7](t)-x[6](t))^2+(y[7](t)
-y[6](t))^2-1, (x[8](t)-x[7](t))^2+(y[8](t)-y[7](t))^2-1, (x[9](t)-x[8](t))^2+
(y[9](t)-y[8](t))^2-1, (x[10](t)-x[9](t))^2+(y[10](t)-y[9](t))^2-1,
x[1](0) = 1, x[2](0) = 2, x[3](0) = 3, x[4](0) = 4, x[5](0) = 5, x[6](0) = 6,
x[7](0) = 7, x[8](0) = 8, x[9](0) = 9, x[10](0) = 10, y[1](0) = 0, y[2](0) = 0
, y[3](0) = 0, y[4](0) = 0, y[5](0) = 0, y[6](0) = 0, y[7](0) = 0, y[8](0) = 0
, y[9](0) = 0, y[10](0) = 0, D(x[1])(0) = 0, D(x[2])(0) = 0, D(x[3])(0) = 0, D
(x[4])(0) = 0, D(x[5])(0) = 0, D(x[6])(0) = 0, D(x[7])(0) = 0, D(x[8])(0) = 0,
D(x[9])(0) = 0, D(x[10])(0) = 0, D(y[1])(0) = 0, D(y[2])(0) = 0, D(y[3])(0) =
0, D(y[4])(0) = 0, D(y[5])(0) = 0, D(y[6])(0) = 0, D(y[7])(0) = 0, D(y[8])(0)
= 0, D(y[9])(0) = 0, D(y[10])(0) = 0}:

  

And the solution:

dsn := dsolve(dsys,numeric,stiff=true,implicit=true);

dsn:=procx_rosenbrock_dae...end proc

(51)

tt := time():

dsn(1):

time()-tt;

1.461

(52)
  

which is less than half the time required in Maple 10.

See Also

DEtools

dsolve

Index of New Maple 11 Features

PDEtools

pdsolve