Product Suggestions

Post your suggestions on new features and products.

eulermac(1/(n*ln(n)^2),n=2..N,1);  #Error
Error, (in SumTools:-DefiniteSum:-ClosedForm) summand is singular in the interval of summation


eulermac(1/(n*ln(n)^2+1),n=2..N,1);  #nonsense

 

 

The method of solving underdetermined systems of equations, and universal method for calculating link mechanisms. It is based on the Draghilev’s method for solving systems of nonlinear equations. 
When calculating link mechanisms we can use geometrical relationships to produce their mathematical models without specifying the “input link”. The new method allows us to specify the “input link”, any link of mechanism.

Example.
Three-bar mechanism.  The system of equations linkages in this mechanism is as follows:

f1 := x1^2+(x2+1)^2+(x3-.5)^2-R^2;
f2 := x1-.5*x2+.5*x3;
f3 := (x1-x4)^2+(x2-x5)^2+(x3-x6)^2-19;
f4 := sin(x4)-x5;
f5 := sin(2*x4)-x6;

Coordinates green point x'i', i = 1..3, the coordinates of red point x'i', i = 4..6.
Set of x0'i', i = 1..6 searched arbitrarily, is the solution of the system of equations and is the initial point for the solution of the ODE system. The solution of ODE system is the solution of system of equations linkages for concrete assembly linkage.
Two texts of the program for one mechanism. In one case, the “input link” is the red-green, other case the “input link” is the green-blue.
After the calculation trajectories of points, we can always find the values of other variables, for example, the angles.
Animation displays the kinematics of the mechanism.
MECAN_3_GR_P_bar.mw 
MECAN_3_Red_P_bar.mw

(if to use another color instead of color = "Niagara Dark Orchid", the version of Maple <17)

Method_Mechan_PDF.pdf






The method of solving underdetermined systems of equations, and universal method for calculating link mechanisms. It is based on the Draghilev’s method for solving systems of nonlinear equations. 
When calculating link mechanisms we can use geometrical relationships to produce their mathematical models without specifying the “input link”. The new method allows us to specify the “input link”, any link of mechanism.

Example.
Three-bar mechanism.  The system of equations linkages in this mechanism is as follows:

f1 := x1^2+(x2+1)^2+(x3-.5)^2-R^2;
f2 := x1-.5*x2+.5*x3;
f3 := (x1-x4)^2+(x2-x5)^2+(x3-x6)^2-19;
f4 := sin(x4)-x5;
f5 := sin(2*x4)-x6;

Coordinates green point x'i', i = 1..3, the coordinates of red point x'i', i = 4..6.
Set of x0'i', i = 1..6 searched arbitrarily, is the solution of the system of equations and is the initial point for the solution of the ODE system. The solution of ODE system is the solution of system of equations linkages for concrete assembly linkage.
Two texts of the program for one mechanism. In one case, the “input link” is the red-green, other case the “input link” is the green-blue.
After the calculation trajectories of points, we can always find the values of other variables, for example, the angles.
Animation displays the kinematics of the mechanism.
MECAN_3_GR_P_bar.mw 
MECAN_3_Red_P_bar.mw

(if to use another color instead of color = "Niagara Dark Orchid", the version of Maple <17)

Method_Mechan_PDF.pdf






As the year draws to a close, we start looking forward to a new year and a new release of Maple. With every new release comes many new features and updates to explore.

We are looking for several new beta testers with a good working knowledge of Maple; We need your input, your ideas, and your experience with our products to help us improve the software and get it ready for general release.

There are many benefits to becoming a beta tester:

  • You’ll get to use the new software before anyone else does.
  • You’ll help us make our software better in ways that work for you.
  • Your suggestions could determine the future direction of the software.
  • You’ll get feedback right from the development team.

If you are interested in becoming a beta tester for the next version of Maple, please email: beta (at) maplesoft.com for more information.

As the year draws to a close, we start looking forward to a new year and a new release of Maple. With every new release comes many new features and updates to explore.

We are looking for several new beta testers with a good working knowledge of Maple; We need your input, your ideas, and your experience with our products to help us improve the software and get it ready for general release.

There are many benefits to becoming a beta tester:

  • You’ll get to use the new software before anyone else does.
  • You’ll help us make our software better in ways that work for you.
  • Your suggestions could determine the future direction of the software.
  • You’ll get feedback right from the development team.

If you are interested in becoming a beta tester for the next version of Maple, please email: beta (at) maplesoft.com for more information.

Mathematica 10.3.0 was announced yesterday. This is the 6th release of Mathematica 10 during 16 months. I wonder its  MathematicaFunctionData and   FindFormula . At first sight, the former is an analog of FunctionAdvisor of Maple, but the latter isn't any analog. Also compare the outputs of

Residue[Binomial[n,k],{n,-j}]

(-1)^j/(j!*k!*(-j-k)!)

 and

>`assuming`([residue(binomial(n, k), n = -j)], [integer, j > 0]);

                residue(binomial(n, k), n = -j)
Let us wait for Maple 2016.

 

Solving DAEs in Maple

As I had mentioned in many posts, Maple cannot solve nonlinear DAEs. As of today (Maple 2015), given a system of index 1 DAE dy/dt = f(y,z); 0 = g(y,z), Maple “extends” g(y,z) to get dz/dt = g1(y,z). So, any given index 1 DAE is converted to a system of ODEs dy/dt = f, dz/dt = g1 with the constraint g = 0, before it solves. This is true for all the solvers in Maple despite the wrong claims in the help files. It is also true for MEBDFI, the FORTRAN implementation of which actually solves index 2 DAEs directly. In addition, the initial condition for the algebraic variable has to be consistent. The problem with using fsolve is that one cannot specify tolerance. Often times one has to solve DAEs at lower tolerances (1e-3 or 1e-4), not 1e-15. In addition, one cannot use compile =true for high efficiency.

The approach in Maple fails for many DAEs of practical interest. In addition, this can give wrong answers. For example, dy/dt = z, y^2+z^2-1=0, with y(0)=0 and z(0)=1 has a meaningful solution only from 0 to Pi/2, Maple’s dsolve/numeric will convert this to dy/dt = z and dz/dt = -y and integrate in time for ever. This is wrong as at t = Pi/2, the system becomes index 2 DAE and there is more than 1 acceptable solution.

We just recently got our paper accepted that helps Maple's dsolve and many ODE solvers in other languages handle DAEs directly. The approach is rather simple, the index 1 DAE is perturbed as dy/dt = f.S ; -mu.diff(g,t) = g. The switch function is a tanh function which is zero till t = tinit (initialization time). Mu is chosen to be a small number. In addition, lhs of the perturbed equation is simplified using approximate initial guesses as Maple cannot handle non-constant Mass matrix. The paper linked below gives below more details.

http://depts.washington.edu/maple/pubs/U_Apprach_FULL_DRAFT.pdf  

Next, I discuss different examples (code attached).

Example 1: Simple DAE (one ODE and one AE), the proposed approach helps solve this system efficiently without knowing the exact initial condition.

Hint: The code is given with a semicolon at the end of the most of the statements for educational purposes for new users. Using semicolon after dsolve numeric in classic worksheet crashes the code (Maplesoft folks couldn’t fix this despite my request).

Example 2:

This is a nickel battery model (one ODE and one AE). This fails with many solvers unless exact IC is given. The proposed approach works well. In particular, stiff=true, implicit=true option is found to be very efficient. The code given in example 1 can be used to solve example 2 or any other DAEs by just entering ODEs, AEs, ICs in proper places.

Example 3:

This is a nonlinear implicit ODE (posted in Mapleprimes earlier by joha, (http://www.mapleprimes.com/questions/203096-Solving-Nonlinear-ODE#answer211682 ). This can be converted to index 1 DAE and solved using the proposed approach.

Example 4:

This example was posted earlier by me in Mapleprimes (http://www.mapleprimes.com/posts/149877-ODEs-PDEs-And-Special-Functions) . Don’t try to solve this in Maple using its dsolve numeric solver for N greater than 5 directly. The proposed approach can handle this well. Tuning the perturbation parameters and using compile =true will help in solving this system more efficiently.

Finally example 1 is solved for different perturbation parameters to show how one can arrive at convergence. Perhaps more sophisticated users and Maplesoft folks can enable this as automatically tuned parameters in dsolve/numeric.

Note:

This post should not be viewed as just being critical on Maple. I have wasted a lot of time assuming that Maple does what it claims in its help file. This post is to bring awareness about the difficulty in dealing with DAEs. It is perfectly fine to not have a DAE solver, but when literature or standard methods are cited/claimed, they should be implemented in proper form. I will forever remain a loyal Maple user as it has enabled me to learn different topics efficiently. Note that Maplesim can solve DAEs, but it is a pain to create a Maplesim model/project just for solving a DAE. Also, events is a pain with Maplesim. The reference to Mapleprimes links are missing in the article, but will be added before the final printing stage. The ability of Maple to find analytical Jacobian helps in developing many robust ODE and DAE solvers and I hope to post my own approaches that will solve more complicated systems.

I strongly encourage testing of the proposed approach and implementation for various educational/research purposes by various users. The proposed approach enables solving lithium-ion and other battery/electrochemical storage models accurately in a robust manner. A disclosure has been filed with the University of Washington to apply for a provisional patent for battery models and Battery Management System for transportation, storage and other applications because of the current commercial interest in this topic (for batteries). In particular, use of this single step avoids intialization issues/(no need to initialize separately) for parameter estimation, state estimation or optimal control of battery models.

 

Appendix A

Maple code for Examples 1-4 from "Extending Explicit and Linealry Implicit ODE Solvers for Index-1 DAEs", M. T. Lawder,

V. Ramadesigan, B. Suthar and V. R. Subramanian, Computers and Chemical Engineering, in press (2015).

Use y1, y2, etc. for all differential variables and z1, z2, etc. for all algebraic variables

 

Example 1

restart;

with(plots):

Enter all ODEs in eqode

eqode:=[diff(y1(t),t)=-y1(t)^2+z1(t)];

eqode := [diff(y1(t), t) = -y1(t)^2+z1(t)]

(1)

Enter all AEs in eqae

eqae:=[cos(y1(t))-z1(t)^0.5=0];

eqae := [cos(y1(t))-z1(t)^.5 = 0]

(2)

Enter all initial conditions for differential variables in icodes

icodes:=[y1(0)=0.25];

icodes := [y1(0) = .25]

(3)

Enter all intial conditions for algebraic variables in icaes

icaes:=[z1(0)=0.8];

icaes := [z1(0) = .8]

(4)

Enter parameters for perturbation value (mu), switch function (q and tint), and runtime (tf)

pars:=[mu=0.1,q=1000,tint=1,tf=5];

pars := [mu = .1, q = 1000, tint = 1, tf = 5]

(5)

Choose solving method (1 for explicit, 2 for implicit)

Xexplicit:=2;

Xexplicit := 2

(6)

Standard solver requires IC z(0)=0.938791 or else it will fail

solx:=dsolve({eqode[1],eqae[1],icodes[1],icaes[1]},numeric):

Error, (in dsolve/numeric/DAE/checkconstraints) the initial conditions do not satisfy the algebraic constraints
  error = .745e-1, tolerance = .559e-6, constraint = cos(y1(t))-z1(t)^.5000000000000000000000

 

ff:=subs(pars,1/2+1/2*tanh(q*(t-tint)));

ff := 1/2+(1/2)*tanh(1000*t-1000)

(7)

NODE:=nops(eqode);NAE:=nops(eqae);

NODE := 1

NAE := 1

(8)

for XX from 1 to NODE do
EQODE||XX:=lhs(eqode[XX])=rhs(eqode[XX])*ff;
end do;

EQODE1 := diff(y1(t), t) = (-y1(t)^2+z1(t))*(1/2+(1/2)*tanh(1000*t-1000))

(9)

for XX from 1 to NAE do
EQAE||XX:=subs(pars,-mu*(diff(rhs(eqae[XX])-lhs(eqae[XX]),t))=rhs(eqae[XX])-lhs(eqae[XX]));
end do;

EQAE1 := -.1*sin(y1(t))*(diff(y1(t), t))-0.5e-1*(diff(z1(t), t))/z1(t)^.5 = -cos(y1(t))+z1(t)^.5

(10)

 

Dvars1:={seq(diff(z||x(t),t)=D||x,x=1..NAE)};

Dvars1 := {diff(z1(t), t) = D1}

(11)

Dvars2:={seq(rhs(Dvars1[x])=lhs(Dvars1[x]),x=1..NAE)};

Dvars2 := {D1 = diff(z1(t), t)}

(12)

icsn:=seq(subs(y||x(0)=y||x(t),icodes[x]),x=1..NODE),seq(subs(z||x(0)=z||x(t),icaes[x]),x=1..NAE);

icsn := y1(t) = .25, z1(t) = .8

(13)

for j from 1 to NAE do

EQAEX||j:=subs(Dvars1,eqode,icsn,Dvars2,lhs(EQAE||j))=rhs(EQAE||j);

end do:

Sys:={seq(EQODE||x,x=1..NODE),seq(EQAEX||x,x=1..NAE),seq(icodes[x],x=1..NODE),seq(icaes[x],x=1..NAE)};

Sys := {-0.1824604200e-1-0.5590169945e-1*(diff(z1(t), t)) = -cos(y1(t))+z1(t)^.5, y1(0) = .25, z1(0) = .8, diff(y1(t), t) = (-y1(t)^2+z1(t))*(1/2+(1/2)*tanh(1000*t-1000))}

(14)

if Xexplicit=1 then

sol:=dsolve(Sys,numeric):

else

sol:=dsolve(Sys,numeric,stiff=true,implicit=true):
end if:

 

for XX from 1 to NODE do
a||XX:=odeplot(sol,[t,y||XX(t)],0..subs(pars,tf),color=red);
end do:

for XX from NODE+1 to NODE+NAE do
a||XX:=odeplot(sol,[t,z||(XX-NODE)(t)],0..subs(pars,tf),color=blue);
end do:

display(seq(a||x,x=1..NODE+NAE),axes=boxed);

 

End Example 1

 

Example 2

restart;

with(plots):

eq1:=diff(y1(t),t)=j1*W/F/rho/V;

eq1 := diff(y1(t), t) = j1*W/(F*rho*V)

(15)

eq2:=j1+j2=iapp;

eq2 := j1+j2 = iapp

(16)

j1:=io1*(2*(1-y1(t))*exp((0.5*F/R/T)*(z1(t)-phi1))-2*y1(t)*exp((-0.5*F/R/T)*(z1(t)-phi1)));

j1 := io1*(2*(1-y1(t))*exp(.5*F*(z1(t)-phi1)/(R*T))-2*y1(t)*exp(-.5*F*(z1(t)-phi1)/(R*T)))

(17)

j2:=io2*(exp((F/R/T)*(z1(t)-phi2))-exp((-F/R/T)*(z1(t)-phi2)));

j2 := io2*(exp(F*(z1(t)-phi2)/(R*T))-exp(-F*(z1(t)-phi2)/(R*T)))

(18)

params:={F=96487,R=8.314,T=298.15,phi1=0.420,phi2=0.303,W=92.7,V=1e-5,io1=1e-4,io2=1e-10,iapp=1e-5,rho=3.4};

params := {F = 96487, R = 8.314, T = 298.15, V = 0.1e-4, W = 92.7, io1 = 0.1e-3, io2 = 0.1e-9, rho = 3.4, iapp = 0.1e-4, phi1 = .420, phi2 = .303}

(19)

eqode:=[subs(params,eq1)];

eqode := [diff(y1(t), t) = 0.5651477584e-2*(1-y1(t))*exp(19.46229155*z1(t)-8.174162450)-0.5651477584e-2*y1(t)*exp(-19.46229155*z1(t)+8.174162450)]

(20)

eqae:=[subs(params,eq2)];

eqae := [0.2e-3*(1-y1(t))*exp(19.46229155*z1(t)-8.174162450)-0.2e-3*y1(t)*exp(-19.46229155*z1(t)+8.174162450)+0.1e-9*exp(38.92458310*z1(t)-11.79414868)-0.1e-9*exp(-38.92458310*z1(t)+11.79414868) = 0.1e-4]

(21)

icodes:=[y1(0)=0.05];

icodes := [y1(0) = 0.5e-1]

(22)

icaes:=[z1(0)=0.7];

icaes := [z1(0) = .7]

(23)

solx:=dsolve({eqode[1],eqae[1],icodes[1],icaes[1]},type=numeric):

Error, (in dsolve/numeric/DAE/checkconstraints) the initial conditions do not satisfy the algebraic constraints
  error = .447e9, tolerance = .880e4, constraint = -2000000*(-1+y1(t))*exp(19.46229155000000000000*z1(t)-8.174162450000000000000)-2000000*y1(t)*exp(-19.46229155000000000000*z1(t)+8.174162450000000000000)+exp(38.92458310000000000000*z1(t)-11.79414868000000000000)-exp(-38.92458310000000000000*z1(t)+11.79414868000000000000)-100000

 

pars:=[mu=0.00001,q=1000,tint=1,tf=5001];

pars := [mu = 0.1e-4, q = 1000, tint = 1, tf = 5001]

(24)

Xexplicit:=2;

Xexplicit := 2

(25)

ff:=subs(pars,1/2+1/2*tanh(q*(t-tint)));

ff := 1/2+(1/2)*tanh(1000*t-1000)

(26)

NODE:=nops(eqode):NAE:=nops(eqae);

NAE := 1

(27)

for XX from 1 to NODE do
EQODE||XX:=lhs(eqode[XX])=rhs(eqode[XX])*ff:
end do;

EQODE1 := diff(y1(t), t) = (0.5651477584e-2*(1-y1(t))*exp(19.46229155*z1(t)-8.174162450)-0.5651477584e-2*y1(t)*exp(-19.46229155*z1(t)+8.174162450))*(1/2+(1/2)*tanh(1000*t-1000))

(28)

for XX from 1 to NAE do
EQAE||XX:=subs(pars,-mu*(diff(rhs(eqae[XX])-lhs(eqae[XX]),t))=rhs(eqae[XX])-lhs(eqae[XX]));
end do;

EQAE1 := -0.2e-8*(diff(y1(t), t))*exp(19.46229155*z1(t)-8.174162450)+0.3892458310e-7*(1-y1(t))*(diff(z1(t), t))*exp(19.46229155*z1(t)-8.174162450)-0.2e-8*(diff(y1(t), t))*exp(-19.46229155*z1(t)+8.174162450)+0.3892458310e-7*y1(t)*(diff(z1(t), t))*exp(-19.46229155*z1(t)+8.174162450)+0.3892458310e-13*(diff(z1(t), t))*exp(38.92458310*z1(t)-11.79414868)+0.3892458310e-13*(diff(z1(t), t))*exp(-38.92458310*z1(t)+11.79414868) = 0.1e-4-0.2e-3*(1-y1(t))*exp(19.46229155*z1(t)-8.174162450)+0.2e-3*y1(t)*exp(-19.46229155*z1(t)+8.174162450)-0.1e-9*exp(38.92458310*z1(t)-11.79414868)+0.1e-9*exp(-38.92458310*z1(t)+11.79414868)

(29)

Dvars1:={seq(diff(z||x(t),t)=D||x,x=1..NAE)};

Dvars1 := {diff(z1(t), t) = D1}

(30)

Dvars2:={seq(rhs(Dvars1[x])=lhs(Dvars1[x]),x=1..NAE)};

Dvars2 := {D1 = diff(z1(t), t)}

(31)

icsn:=seq(subs(y||x(0)=y||x(t),icodes[x]),x=1..NODE),seq(subs(z||x(0)=z||x(t),icaes[x]),x=1..NAE);

icsn := y1(t) = 0.5e-1, z1(t) = .7

(32)

for j from 1 to NAE do

EQAEX||j:=subs(Dvars1,eqode,icsn,Dvars2,lhs(EQAE||j))=rhs(EQAE||j);

end do;

EQAEX1 := -0.2e-8*(0.5368903705e-2*exp(5.449441630)-0.2825738792e-3*exp(-5.449441630))*exp(5.449441630)+0.3697835394e-7*(diff(z1(t), t))*exp(5.449441630)-0.2e-8*(0.5368903705e-2*exp(5.449441630)-0.2825738792e-3*exp(-5.449441630))*exp(-5.449441630)+0.1946229155e-8*(diff(z1(t), t))*exp(-5.449441630)+0.3892458310e-13*(diff(z1(t), t))*exp(15.45305949)+0.3892458310e-13*(diff(z1(t), t))*exp(-15.45305949) = 0.1e-4-0.2e-3*(1-y1(t))*exp(19.46229155*z1(t)-8.174162450)+0.2e-3*y1(t)*exp(-19.46229155*z1(t)+8.174162450)-0.1e-9*exp(38.92458310*z1(t)-11.79414868)+0.1e-9*exp(-38.92458310*z1(t)+11.79414868)

(33)

Sys:={seq(EQODE||x,x=1..NODE),seq(EQAEX||x,x=1..NAE),seq(icodes[x],x=1..NODE),seq(icaes[x],x=1..NAE)};

Sys := {-0.5810962488e-6+0.8802389238e-5*(diff(z1(t), t)) = 0.1e-4-0.2e-3*(1-y1(t))*exp(19.46229155*z1(t)-8.174162450)+0.2e-3*y1(t)*exp(-19.46229155*z1(t)+8.174162450)-0.1e-9*exp(38.92458310*z1(t)-11.79414868)+0.1e-9*exp(-38.92458310*z1(t)+11.79414868), y1(0) = 0.5e-1, z1(0) = .7, diff(y1(t), t) = (0.5651477584e-2*(1-y1(t))*exp(19.46229155*z1(t)-8.174162450)-0.5651477584e-2*y1(t)*exp(-19.46229155*z1(t)+8.174162450))*(1/2+(1/2)*tanh(1000*t-1000))}

(34)

if Xexplicit=1 then

sol:=dsolve(Sys,numeric,maxfun=0):

else

sol:=dsolve(Sys,numeric,stiff=true,implicit=true,maxfun=0):

end if:

 

for XX from 1 to NODE do
a||XX:=odeplot(sol,[t,y||XX(t)],0..subs(pars,tf),color=red);
end do:

for XX from NODE+1 to NODE+NAE do
a||XX:=odeplot(sol,[t,z||(XX-NODE)(t)],0..subs(pars,tf),color=blue);
end do:

b1:=odeplot(sol,[t,y1(t)],0..1,color=red):
b2:=odeplot(sol,[t,z1(t)],0..1,color=blue):

display(b1,b2,axes=boxed);

 

display(seq(a||x,x=1..NODE+NAE),axes=boxed);

 

End Example 2

 

Example 3

restart;

with(plots):

eq1:=diff(y1(t),t)^2+diff(y1(t),t)*(y1(t)+1)+y1(t)=cos(diff(y1(t),t));

eq1 := (diff(y1(t), t))^2+(diff(y1(t), t))*(y1(t)+1)+y1(t) = cos(diff(y1(t), t))

(35)

solx:=dsolve({eq1,y1(0)=0},numeric):

Error, (in dsolve/numeric/make_proc) Could not convert to an explicit first order system due to 'RootOf'

 

eqode:=[diff(y1(t),t)=z1(t)];

eqode := [diff(y1(t), t) = z1(t)]

(36)

eqae:=[subs(eqode,eq1)];

eqae := [z1(t)^2+z1(t)*(y1(t)+1)+y1(t) = cos(z1(t))]

(37)

icodes:=[y1(0)=0.0];

icodes := [y1(0) = 0.]

(38)

icaes:=[z1(0)=0.0];

icaes := [z1(0) = 0.]

(39)

pars:=[mu=0.1,q=1000,tint=1,tf=4];

pars := [mu = .1, q = 1000, tint = 1, tf = 4]

(40)

Xexplicit:=2;

Xexplicit := 2

(41)

ff:=subs(pars,1/2+1/2*tanh(q*(t-tint)));

ff := 1/2+(1/2)*tanh(1000*t-1000)

(42)

NODE:=nops(eqode);NAE:=nops(eqae);

NODE := 1

NAE := 1

(43)

for XX from 1 to NODE do
EQODE||XX:=lhs(eqode[XX])=rhs(eqode[XX])*ff:
end do;

EQODE1 := diff(y1(t), t) = z1(t)*(1/2+(1/2)*tanh(1000*t-1000))

(44)

for XX from 1 to NAE do
EQAE||XX:=subs(pars,-mu*(diff(rhs(eqae[XX])-lhs(eqae[XX]),t))=rhs(eqae[XX])-lhs(eqae[XX]));
end do;

EQAE1 := .1*sin(z1(t))*(diff(z1(t), t))+.2*z1(t)*(diff(z1(t), t))+.1*(diff(z1(t), t))*(y1(t)+1)+.1*z1(t)*(diff(y1(t), t))+.1*(diff(y1(t), t)) = cos(z1(t))-z1(t)^2-z1(t)*(y1(t)+1)-y1(t)

(45)

 

Dvars1:={seq(diff(z||x(t),t)=D||x,x=1..NAE)};

Dvars1 := {diff(z1(t), t) = D1}

(46)

Dvars2:={seq(rhs(Dvars1[x])=lhs(Dvars1[x]),x=1..NAE)};

Dvars2 := {D1 = diff(z1(t), t)}

(47)

icsn:=seq(subs(y||x(0)=y||x(t),icodes[x]),x=1..NODE),seq(subs(z||x(0)=z||x(t),icaes[x]),x=1..NAE);

icsn := y1(t) = 0., z1(t) = 0.

(48)

for j from 1 to NAE do

EQAEX||j:=subs(Dvars1,eqode,icsn,Dvars2,lhs(EQAE||j))=rhs(EQAE||j);

end do;

EQAEX1 := .1*sin(0.)*(diff(z1(t), t))+.1*(diff(z1(t), t)) = cos(z1(t))-z1(t)^2-z1(t)*(y1(t)+1)-y1(t)

(49)

Sys:={seq(EQODE||x,x=1..NODE),seq(EQAEX||x,x=1..NAE),seq(icodes[x],x=1..NODE),seq(icaes[x],x=1..NAE)};

Sys := {.1*(diff(z1(t), t)) = cos(z1(t))-z1(t)^2-z1(t)*(y1(t)+1)-y1(t), y1(0) = 0., z1(0) = 0., diff(y1(t), t) = z1(t)*(1/2+(1/2)*tanh(1000*t-1000))}

(50)

if Xexplicit=1 then

sol:=dsolve(Sys,numeric):

else

sol:=dsolve(Sys,numeric,stiff=true,implicit=true):

end if:

 

for XX from 1 to NODE do
a||XX:=odeplot(sol,[t,y||XX(t)],0..subs(pars,tf),color=red);
end do:

for XX from NODE+1 to NODE+NAE do
a||XX:=odeplot(sol,[t,z||(XX-NODE)(t)],0..subs(pars,tf),color=blue);
end do:

display(seq(a||x,x=1..NODE+NAE),axes=boxed);

 

End Example 3

 

Example 4

restart;

with(plots):

N:=11:h:=1/(N+1):

for i from 2 to N+1 do eq1[i]:=diff(y||i(t),t)=(y||(i+1)(t)-2*y||i(t)+y||(i-1)(t))/h^2-y||i(t)*(1+z||i(t));od:

for i from 2 to N+1 do eq2[i]:=0=(z||(i+1)(t)-2*z||i(t)+z||(i-1)(t))/h^2-(1-y||i(t)^2)*(exp(-z||i(t)));od:

eq1[1]:=(3*y1(t)-4*y2(t)+y3(t))/(2*h)=0: eq1[N+2]:=y||(N+2)(t)-1=0:

eq2[1]:=(3*z1(t)-4*z2(t)+1*z3(t))/(2*h)=0: eq2[N+2]:=z||(N+2)(t)=0:

eq1[1]:=subs(y1(t)=z||(N+3)(t),eq1[1]):

eq1[N+2]:=subs(y||(N+2)(t)=z||(N+4)(t),eq1[N+2]):

eqode:=[seq(subs(y1(t)=z||(N+3)(t),y||(N+2)(t)=z||(N+4)(t),eq1[i]),i=2..N+1)]:

eqae:=[eq1[1],eq1[N+2],seq(eq2[i],i=1..N+2)]:

icodes:=[seq(y||j(0)=1,j=2..N+1)]:

icaes:=[seq(z||j(0)=0,j=1..N+2),z||(N+3)(0)=1,z||(N+4)(0)=1]:

pars:=[mu=0.00001,q=1000,tint=1,tf=2]:

Xexplicit:=2:

ff:=subs(pars,1/2+1/2*tanh(q*(t-tint))):

NODE:=nops(eqode):NAE:=nops(eqae):

for XX from 1 to NODE do

EQODE||XX:=lhs(eqode[XX])=rhs(eqode[XX])*ff: end do:

for XX from 1 to NAE do

EQAE||XX:=subs(pars,-mu*(diff(rhs(eqae[XX])-lhs(eqae[XX]),t))=rhs(eqae[XX])-lhs(eqae[XX])); end do:

Dvars1:={seq(diff(z||x(t),t)=D||x,x=1..NAE)}:

Dvars2:={seq(rhs(Dvars1[x])=lhs(Dvars1[x]),x=1..NAE)}:

icsn:=seq(subs(y||x(0)=y||x(t),icodes[x]),x=1..NODE),seq(subs(z||x(0)=z||x(t),icaes[x]),x=1..NAE):

for j from 1 to NAE do

EQAEX||j:=subs(Dvars1,eqode,icsn,Dvars2,lhs(EQAE||j))=rhs(EQAE||j):

end do:

Sys:={seq(EQODE||x,x=1..NODE),seq(EQAEX||x,x=1..NAE),seq(icodes[x],x=1..NODE),seq(icaes[x],x=1..NAE)}:

if Xexplicit=1 then

sol:=dsolve(Sys,numeric,maxfun=0):

else

sol:=dsolve(Sys,numeric,stiff=true,implicit=true,maxfun=0):

end if:

 

for XX from 1 to NODE do

a||XX:=odeplot(sol,[t,y||(XX+1)(t)],1..subs(pars,tf),color=red): end do:

for XX from NODE+1 to NODE+NAE do

a||XX:=odeplot(sol,[t,z||(XX-NODE)(t)],1..subs(pars,tf),color=blue): end do:

display(seq(a||x,x=1..NODE),a||(NODE+NAE-1),a||(NODE+NAE),axes=boxed);

 

End of Example 4

 

Sometimes the parameters of the switch function and perturbation need to be tuned to obtain propoer convergence. Below is Example 1 shown for several cases using the 'parameters' option in Maple's dsolve to compare how tuning parameters affects the solution

restart:

with(plots):

eqode:=[diff(y1(t),t)=-y1(t)^2+z1(t)]: eqae:=[cos(y1(t))-z1(t)^0.5=0]:

icodes:=[y1(0)=0.25]: icaes:=[z1(0)=0.8]:

pars:=[tf=5]:

Xexplicit:=2;

Xexplicit := 2

(51)

ff:=subs(pars,1/2+1/2*tanh(q*(t-tint))):

NODE:=nops(eqode):NAE:=nops(eqae):

for XX from 1 to NODE do
EQODE||XX:=lhs(eqode[XX])=rhs(eqode[XX])*ff:
end do:

for XX from 1 to NAE do
EQAE||XX:=subs(pars,-mu*(diff(rhs(eqae[XX])-lhs(eqae[XX]),t))=rhs(eqae[XX])-lhs(eqae[XX]));
end do:

 

Dvars1:={seq(diff(z||x(t),t)=D||x,x=1..NAE)}:

Dvars2:={seq(rhs(Dvars1[x])=lhs(Dvars1[x]),x=1..NAE)}:

icsn:=seq(subs(y||x(0)=y||x(t),icodes[x]),x=1..NODE),seq(subs(z||x(0)=z||x(t),icaes[x]),x=1..NAE):

for j from 1 to NAE do

EQAEX||j:=subs(Dvars1,eqode,icsn,Dvars2,lhs(EQAE||j))=rhs(EQAE||j):

end do:

Sys:={seq(EQODE||x,x=1..NODE),seq(EQAEX||x,x=1..NAE),seq(icodes[x],x=1..NODE),seq(icaes[x],x=1..NAE)}:

if Xexplicit=1 then

sol:=dsolve(Sys,numeric,'parameters'=[mu,q,tint],maxfun=0):

else

sol:=dsolve(Sys,numeric,'parameters'=[mu,q,tint],stiff=true,implicit=true):

end if:

 

sol('parameters'=[0.1,1000,1]):

plot1:=odeplot(sol,[t-1,y1(t)],0..4,color=red):
plot2:=odeplot(sol,[t-1,z1(t)],0..4,color=blue):

sol('parameters'=[0.001,10,1]):

plot3:=odeplot(sol,[t-1,y1(t)],0..4,color=red,linestyle=dot):
plot4:=odeplot(sol,[t-1,z1(t)],0..4,color=blue,linestyle=dot):

display(plot1,plot2,plot3,plot4,axes=boxed);

 

In general, one has to decrease mu, and increase q and tint until convergence (example at t=3)

sol('parameters'=[0.001,10,1]):sol(3+1);

[t = 4., y1(t) = .738587929442734, z1(t) = .546472878850096]

(52)

sol('parameters'=[0.0001,100,10]):sol(3+10);

[t = 13., y1(t) = .738684397167344, z1(t) = .546618936273638]

(53)

sol('parameters'=[0.00001,1000,20]):sol(3+20);

[t = 23., y1(t) = .738694113087217, z1(t) = .546633473784526]

(54)

 

The results have converged to 4 digits after the decimal. Of course, absolute and relative tolerances of the solvers can be modified if needed

 

Download SolvingDAEsinMaple.mws

This is a very simple suggestion that weights heavily on my enjoyment of Maple.

I'm not sure if it is only me and my students but it is really tricky to resize a code edit region. Trying to get a hold of the contour in an exercise in patience! Can anyone fix this?

 

Thanks!

Stéphane

I do not think the current API for dsolve when asking for series solution is done right. If one wants to obtain a series solution for an ODE, but wants different order than the default 6, now one must set this value using a global setting before making the call, like this:

eq:=diff(y(x),x$2)+y(x)=0;
Order:=10;
dsolve({eq,y(0)=1,D(y)(0)=0},y(x),type='series');

It would be better if options to calls are passed along with the other parameters in the call itself. Something like

dsolve({eq,y(0)=1,D(y)(0)=0},y(x),type='series',order=10);

This should also apply to any command that takes Order, such as series(cos(x),x=0,order=10);

Passing options and values to functions using global and environment variables is not safe and not a good way to go about it as it can cause programming errors.

@ecterrab 

I figured I'd start a new thread for odd things I come across whilst using the new physics package. 

I have found this, and am not sure if it is expected. 

 


restart

with(Physics):

Setup(mathematicalnotation = true):

``

Setup(Commutator(Psigma[i], Psigma[j]) = Physics:-`*`(Physics:-`*`(I, ep_[i, j, k]), Psigma[k]), AntiCommutator(Psigma[i], Psigma[j]) = Physics:-`*`(2, kd_[i, j]));

[algebrarules = {%AntiCommutator(Physics:-Psigma[i], Physics:-Psigma[j]) = 2*Physics:-KroneckerDelta[i, j], %Commutator(Physics:-Psigma[i], Physics:-Psigma[j]) = I*Physics:-LeviCivita[i, j, k]*Physics:-Psigma[k]}]

(1)

NULL

Psigma[1].Psigma[1]

Physics:-Psigma[1]^2

(2)

Simplify(%)

Physics:-Psigma[1]^2

(3)

Simplify(Physics:-Psigma[1]^2)

1

(4)

``


Download Simplify2.mw

Did Maplesoft make a statistics about questions on MaplePrimes?

There must be some commands frequently meet problems.

I listed 10 commands or packages need to be enhanced much more.

1.fsolve.

eg, roots of transcendental equation on complex plane & system of nonlinear equations.

2.int.

Both numerical integration & symbolic integration need further developed.

3.Eigenvalues&LinearSolve.

large scale computing, everyone needs it.

4.dsolve&pdsolve.

dsolve needs to be enhanced for solve more equations in parallel.

As mathematica 10 contains the FEM package, I think Maple needs to do something, develop a new package needs more time than make the existing pdsolve command better.

5.Threads package.

parallel programing is the future. BTW, why don't make the task model much easier to use, like the mathematica, just a command?

6.plot&plot3d

there is much space to go further about the two commands.

7.Optimization package or Statistics package.

In fact, both packages need to be further enhanced, especially for Optimization package.

Although some modern algorithms about optimization&statistics is much easier to see on the internet, the

Maplesoft should not stop, but go further, add them in Maple.

Above is just my opinion.

Thanks for your attention.

The problem is when initializing a Matrix with a list of strings. The worksheet excerpt below shows the normal behavior using a list of integers to initialize a square matrix: the successive list elements fill the matrix by rows.

Then trying the exact same thing with a list of strings instead of integers gives an error message!
This is not right. While it is an odd and likely rare problem, it would be better fixed.

x := [i $ i=1..25];
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25]
Matrix(5,5,x);
Matrix(5,5,[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25]) 

y := map(convert,x,string);
["1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16", "17", "18", "19", "20", "21", "22", "23", "24", "25"]
Matrix(5,5,y);
Error, (in Matrix) initializer defines more columns (25) than column dimension parameter specifies (5)

Greetings to all.

I would like to share a brief observation concerning my experiences with the Euler-Maclaurin summation routine in Maple 17 (X86 64 LINUX). The following Math StackExchange Link shows how to compute a certain Euler-MacLaurin type asymptotic expansion using highly unorthodox divergent series summation techniques. The result that was obtained matches the output from eulermac which is definitely good to know. What follows is the output from said routine.

> eulermac(1/(1+k/n),k=0..n,18);
     1       929569        3202291        691                O(1)
O(- ---) - ----------- + ----------- - --------- + 1/1048576 ----
     19             15            17          11              19
    n      2097152 n     1048576 n     32768 n               n

                                           n
                                          /
        174611      5461        31       |      1           17        1
     - -------- + --------- + ------- +  |   ------- dk - ------- + ------
             19          13         9    |   1 + k/n            7        5
       6600 n     65536 n     4096 n    /                 4096 n    256 n
                                          0

         1       1
     - ------ + ---- + 3/4
            3   16 n
       128 n

While I realize that this is good enough for most purposes I have two minor issues.

  • One could certainly evaluate the integral without leaving it to the user to force evaluation with the AllSolutions option. One can and should make use of what is known about n and k. In particular one can check whether there are singularities on the integration path because we know the range of k/n.
  • Why are there two order terms for the order of the remainder term? There should be at most one and a coefficient times an O(1) term makes little sense as the coefficient would be absorbed.

You might want to fix these so that the output looks a bit more professional which does enter into play when potential future users decide on what CAS to commit to. Other than that it is a very useful routine even for certain harmonic sum computations where one can use Euler-Maclaurin to verify results.

Best regards,

Marko Riedel

Maple 18, document mode. Please try this:

>restart;
>some_long_name:=3:
>f(som)

Now, as the cursor is right after the "m" above, hiting the ESC key to autocomplete, nothing happens. Only a beep, meaning Maple does not know about the variable name I wanted to expand to, which is "some_long_name". What I do now is manually add a space by  pushing the ")" away, like this:

  f(som )

And now put the mouse back in front of the "m" and now hit ESC. Maple now can see the name I wanted.

Really? This is very bad design. You might ask, why did I close the () first? Well, I like to start by writing () then fill it in, or it can be that I am changing things, and wanted to write a new name, and the () was allready there. I use name autocomplete alot, and it is very annoying to not be able to use it if happens that there is a ")" or "]" right next to the name.

Having to keep pushing the ")" away, so that Maple can see the name makes no sense since ")" can't part of a name.

In Mathematica it just works. I can do

someLongName = 3;

Then in a new cell,

f[som] and hit CTRL-K  when mouse at the letter "m" just typed it in, and it will autocomplete.  I did not have to push "]" away first.

Please Maple, fix your user interface so it is less awkward and annoying to use. This is version 18 of the software, and not version 1.

 

A little late in the game, with no Maple 18 wishlist yet I'd like to get it started

granted all wishes for previous versions not yet applied should still, hopefully, be under consideration.  

 

1 - The abitlity to angle x-axis labels.

 

2 - Textplot option for text to rotate with graph as it is rotated.  Currently all textplot texts are held to the horizontal

 

 

 

2 3 4 5 6 7 8 Last Page 4 of 21