Updates to Differential Equation (DE) Solvers in Maple 11
Summary
Exact Solutions
Numeric Solutions and related Plotting
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.
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);
y⁡x⁢will now be displayed as⁢y
derivatives with respect to⁢x⁢of functions of one variable will now be displayed with '
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''=4⁢y'x+6⁢y2x2−6⁢yx2−x2⁢a2
dsolve((2));
y=WeierstrassP⁡x+c__1,a,c__2⁢x2
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);
dsolve((4));
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.
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;
−2⁢yx3+y'x22−1−b2⁢y4x8−−b2−1⁢y2x4=0
dsolve((6), singsol = none);
y=JacobiSN⁡x+c__1,b⁢x2
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;
−2⁢yx3+y'x22+y⁢ax2−4⁢y3x6+b=0
dsolve((8), implicit); # implicit avoids computing the roots of the cube
y3−b⁢x64−y⁢a⁢x44=0,y=WeierstrassP⁡x+c__1,a,b⁢x2
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+a6⁢y+y'⁢x3+a32⁢b3−27⁢x4⁢y42−27⁢b⁢x2⁢y2c−27⁢b22⁢c2−3⁢a⁢x+a4⁢y+y'⁢x22⁢b=0
[dsolve((10), implicit)];
y4−−54⁢b4⁢x2⁢y2⁢c+a3⁢c2−27⁢b527⁢x4⁢b3⁢c2=0,y=WeierstrassPPrime⁡1x+a+c__1,ab,bcx
map(odetest, (11), (10)); # verify these results for correctness
0,0
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.
Two examples involving non-rational coefficients
diff(y(x),x) = y(x)/x/(ln(x)*y(x)+(-ln(x)+b)*(-ln(x)+c));
y'=yx⁢ln⁡x⁢y+−ln⁡x+b⁢−ln⁡x+c
dsolve((13));
c__1+ln⁡x⁢KummerM⁡−b,1+c−b,y−b⁢KummerM⁡−b+1,1+c−b,yln⁡x⁢KummerU⁡−b,1+c−b,y+b⁢c⁢KummerU⁡−b+1,1+c−b,y=0
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'=y⁢y−1⁢ⅇxy⁢a+ⅇ2⁢x−ⅇx⁢c
sol := dsolve((15));
sol≔c__1+HeunCPrime⁡0,−c,−1,0,a+12,yy−1⁢y−ⅇ2⁢x⁢y−1⁢HeunC⁡0,−c,−1,0,a+12,yy−1ⅇx⁢y−1cy−1⁢−ⅇ2⁢xⅇx+c⁢HeunC⁡0,c,−1,0,a+12,yy−1+y⁢HeunCPrime⁡0,c,−1,0,a+12,yy−1⁢yc=0
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'=3⁢x4−6⁢x3+3⁢x2+6⁢x+3⁢y2+2⁢x8−15⁢x7−26⁢x6+3⁢x5+30⁢x4−3⁢x3⁢y−x7⁢x5+5⁢x4+29⁢x3+21⁢x2−27⁢x−31x4+x2+1⁢y+x3⁢x5−x4+x3−3⁢x2+x−1⁢x
dsolve((17));
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).
Kamke's example number 1.83
diff(y(x),x) = tan(x*y(x));
y'=tan⁡x⁢y
dsolve((19), implicit);
−c__1+erf⁡2⁢x+I⁢y2−erf⁡2⁢x−I⁢y2=0
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);
dsolve((21));
y=c__1⁢hypergeom⁡38,58,78,98,34,54,32,−3x4x32+c__2⁢hypergeom⁡18,38,58,78,12,34,54,−3x4x+c__3⁢x⁢hypergeom⁡−18,18,38,58,14,12,34,−3x4+c__4⁢hypergeom⁡58,78,98,118,54,32,74,−3x4x52
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,ChangeSymmetry,CharacteristicQ,CharacteristicQInvariants,ConservedCurrentTest,ConservedCurrents,ConsistencyTest,D_Dx,DeterminingPDE,Eta_k,Euler,FirstIntegralSolver,FromJet,FunctionFieldSolutions,InfinitesimalGenerator,Infinitesimals,IntegratingFactorTest,IntegratingFactors,InvariantEquation,InvariantSolutions,InvariantTransformation,Invariants,Laplace,Library,PDEplot,PolynomialSolutions,ReducedForm,SimilaritySolutions,SimilarityTransformation,Solve,SymmetryCommutator,SymmetryGauge,SymmetrySolutions,SymmetryTest,SymmetryTransformation,TWSolutions,ToJet,ToMissingDependentVariable,build,casesplit,charstrip,dchange,dcoeffs,declare,diff_table,difforder,dpolyform,dsubs,mapde,separability,splitstrip,splitsys,undeclare
declare((u,v)(x,y,z));
u⁡x,y,z⁢will now be displayed as⁢u
v⁡x,y,z⁢will now be displayed as⁢v
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];
sys≔ux,y,z⁢v+ux,y⁢vz+ux,z⁢vy+ux⁢vy,z+uy,z⁢vx+uy⁢vx,z+uz⁢vx,y+u⁢vx,y,z=0,ux,y,z+vx,y,z=0
As in the case of a single PDE, for PDE systems, unless indicated otherwise, the test is regarding separability by sum.
separability(sys);
0
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, `*`)];
remain≔ux⁢vx⁢uz2⁢v2+vz2⁢u2⁢uy2⁢v2+vy2⁢u2⁢ux⁢uy⁢uz⁢v2+vx⁢vy⁢vz⁢u2,uy⁢vy⁢uz2⁢v2+vz2⁢u2⁢ux2⁢v2+vx2⁢u2⁢ux⁢uy⁢uz⁢v2+vx⁢vy⁢vz⁢u2,uz⁢vz⁢uy2⁢v2+vy2⁢u2⁢ux2⁢v2+vx2⁢u2⁢ux⁢uy⁢uz⁢v2+vx⁢vy⁢vz⁢u2
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)];
restriction≔u=u⁡x,z,v=v⁡y,z
Evaluate remain at this restriction
eval(remain, restriction);
0,0,0
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__1⁡x+f__2⁡y+f__3⁡z,v=f__4⁡x+f__5⁡y+f__6⁡z
pdetest((30), sys); # verify results - see ?pdetest
pdsolve(sys, HINT = `*`);
u=f__1⁡x⁢f__2⁡y⁢c__1,v=f__4⁡x⁢f__5⁡y⁢c__2
pdetest((32), sys);
Symmetry analysis and group invariant solutions for PDE systems
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));
u⁡x,t⁢will now be displayed as⁢u
PDE := diff(u(x,t),x,x) - diff(u(x,t), t) = 0;
PDE≔ux,x−ut=0
The determining PDE system whose solutions are the infinitesimals
DetSys := DeterminingPDE(PDE);
DetSys≔_ξtt,t,t=0,_ηut,u=−_ξtt,t4,_ηuu,u=0,_ηux,x=_ηut,_ξtu=0,_ξtx=0,_ξxt=−2⁢_ηuu,x,_ξxu=0,_ξxx=_ξtt2
You can now compute the infinitesimals of the generators of point symmetry transformations of PDE by solving this system with pdsolve.
pdsolve(DetSys);
_ηu⁡x,t,u=−x2−2⁢t⁢c__18+c__4⁢x+c__5⁢u+c__8⁢ⅇ_c1⁢t⁢c__6⁢ⅇ_c1⁢x+c__8⁢ⅇ_c1⁢t⁢c__7ⅇ_c1⁢x,_ξt⁡x,t,u=12⁢c__1⁢t2+c__2⁢t+c__3,_ξx⁡x,t,u=c__1⁢x−4⁢c__4⁢t2+c__2⁢x2+c__9
Alternatively, you can call the new Infinitesimals command with PDE as the argument.
G := Infinitesimals(PDE);
G≔_ξx⁡x,t,u=0,_ξt⁡x,t,u=1,_ηu⁡x,t,u=0,_ξx⁡x,t,u=1,_ξt⁡x,t,u=0,_ηu⁡x,t,u=0,_ξx⁡x,t,u=0,_ξt⁡x,t,u=0,_ηu⁡x,t,u=u,_ξx⁡x,t,u=x2,_ξt⁡x,t,u=t,_ηu⁡x,t,u=0,_ξx⁡x,t,u=−2⁢t,_ξt⁡x,t,u=0,_ηu⁡x,t,u=u⁢x,_ξx⁡x,t,u=0,_ξt⁡x,t,u=0,_ηu⁡x,t,u=ⅇ_c1⁢x⁢ⅇ_c1⁢t,_ξx⁡x,t,u=0,_ξt⁡x,t,u=0,_ηu⁡x,t,u=ⅇ_c1⁢tⅇ_c1⁢x,_ξx⁡x,t,u=t⁢x2,_ξt⁡x,t,u=t22,_ηu⁡x,t,u=−x2+2⁢t⁢u8
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,s=t+_ε,v⁡r,s=u
Solve this transformation for {x, t, u}
TR := solve((39), {x,t, u(x,t)});
TR≔t=s−_ε,x=r,u=v⁡r,s
Applying now this transformation to PDE (see dchange) you obtain again PDE, that is, the meaning of "transformation leaving invariant the PDE."
PDE;
ux,x−ut=0
dchange(TR, (41), [r,s, v(r,s)]);
vr,r−vs=0
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,s=t,v⁡r=u,t=s,x=r,u=v⁡r
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
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,v⁡r,s=u,t=r,u=v⁡r,s,r=x,v⁡r,s=u,x=r,u=v⁡r,s,r=tx2,s=t,v⁡r,s=u,t=s,x=sr,u=v⁡r,s,r=t,s=x,v⁡r,s=uⅇ−x24⁢t,t=r,x=s,u=v⁡r,s⁢ⅇ−s24⁢r,r=tx,s=t,v⁡r,s=u⁢xⅇ−x24⁢t,t=s,x=sr,u=v⁡r,s⁢ⅇ−s4⁢r2sr
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,u=v⁡r,s,r,v⁡r,s
dchange((46)[1], PDE, (46)[2]);
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,u=c__1⁢x+c__2,u=c__1t⁢ⅇx24⁢t,u=c__1+erf⁡12⁢tx2⁢c__2,u=c__1⁢x+c__2⁢tt⁢tx⁢x⁢ⅇx24⁢t
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.
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.,alpha_1⁡t=−1.41404214314427,ⅆⅆtalpha_1⁡t=0.00849473764599224,alpha_2⁡t=2.46031987052574,ⅆⅆtalpha_2⁡t=−0.00427806160361022,alpha_3⁡t=2.76667851916170,ⅆⅆtalpha_3⁡t=−0.00224383300346078,alpha_4⁡t=0.373189183301449,ⅆⅆtalpha_4⁡t=−0.00345824250260502,alpha_5⁡t=0.676103669084111,ⅆⅆtalpha_5⁡t=−0.0126138578651418,alpha_6⁡t=−1.72061386826404,ⅆⅆtalpha_6⁡t=0.0145971225866554,beta_1⁡t=0.399270368285695,ⅆⅆtbeta_1⁡t=0.0168856342949760,beta_2⁡t=0.395754725106839,ⅆⅆtbeta_2⁡t=0.00854132788090072,beta_3⁡t=0.390747038446706,ⅆⅆtbeta_3⁡t=0.0115082928413281,beta_4⁡t=0.394386856258629,ⅆⅆtbeta_4⁡t=0.0218166276873265,beta_5⁡t=0.407005428357208,ⅆⅆtbeta_5⁡t=0.0389816800352700,beta_6⁡t=0.406830807784426,ⅆⅆtbeta_6⁡t=0.0366850260352228,eta_P⁡t=−0.0422894233981797,ⅆⅆteta_P⁡t=−0.105787671732061,λ1⁡t=−2.10919513521311,λ10⁡t=−10.4984524640529,λ11⁡t=−4.15923393631383,λ12⁡t=−16.1304320809642,λ13⁡t=−9.02987945391216,λ14⁡t=−7.31967620701184,λ15⁡t=−16.2842801493366,λ16⁡t=1.74150222732886,λ17⁡t=11.3835236220591,λ18⁡t=−15.8371858927843,λ2⁡t=13.0088590989426,λ3⁡t=−19.0184790709148,λ4⁡t=10.0748474937966,λ5⁡t=−8.24678875909166,λ6⁡t=−18.9498152109804,λ7⁡t=10.2437486523454,λ8⁡t=−4.11354762594244,λ9⁡t=−15.6064276928410,s_1⁡t=2.35208666252732,ⅆⅆts_1⁡t=−0.118897311472506,s_2⁡t=2.38757724404409,ⅆⅆts_2⁡t=−0.0809563086256958,s_3⁡t=2.42579134679721,ⅆⅆts_3⁡t=−0.105119915874658,s_4⁡t=2.40699590616213,ⅆⅆts_4⁡t=−0.150615237136879,s_5⁡t=2.31403328393115,ⅆⅆts_5⁡t=−0.262431777006591,s_6⁡t=2.29725446223393,ⅆⅆts_6⁡t=−0.255377349604363,x_P⁡t=0.00332650598171118,ⅆⅆtx_P⁡t=0.0111590920256401,xi_P⁡t=−0.00124602843239396,ⅆⅆtxi_P⁡t=−0.00409946001158069,y_P⁡t=0.00368685194005796,ⅆⅆty_P⁡t=−0.00260844663960202,z_P⁡t=2.94202195913911,ⅆⅆtz_P⁡t=−0.177126550474447,zeta_P⁡t=−0.0583839149780300,ⅆⅆtzeta_P⁡t=−0.0400305840410983
time()-tt;
3.827
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
dsn(1):
1.461
which is less than half the time required in Maple 10.
See Also
DEtools
dsolve
Index of New Maple 11 Features
PDEtools
pdsolve
Download Help Document