1)

I think that 1) is showing a bug in the _d01ajc integrator. How to workaround that best could depend on what you are really trying to do. I'll try and describe some workarounds.

You've set up your procedure `f` as recursive, just to get periodic behaviour. But if you know the desired range of integration in advance, then you could build just a simpler non-recursive integrand up front. (I can't tell whether your posted example is what you really want to work with, or just a illustrative toy example.)

Many (but not all) quadrature schemes will have accuracy that depends upon the integrand's being smooth, and you integrand is only piecewise smooth. Now, evaf/Int can try and split the integration across "discontinuties" but it might not do that if either the method is specified or if the integrand is supplied as a procedure (or maybe even an unevaluated function call). That makes your recursive `f` problematic, because the ways to normally prevent Maple from going into infinite recursion while trying to eveluate the integrand as procedure argument are also things that make evalf/Int not look for discontinuties and split the region.

You might want to look these over, and take your pick.

a) Force _Dexp method, which produces a rough result and is not very fast,

restart:
f := x->piecewise(x>=0 and 4>x,50,x>=4 and 5>x,0,x>=5,f(x-5)):
plot(Int(f, 0 .. t, method=_Dexp), t=0 .. 20);

b) Force use of _d01akc which is NAG integrator for periodic integrands. This seems to do well, being a smooth result and pretty fast. (Is that by accident, or because integrand actually is periodic if not differentiable everywhere?)

restart:
f := x->piecewise(x>=0 and 4>x,50,x>=4 and 5>x,0,x>=5,f(x-5)):
plot(Int(f, 0 .. t, method=_d01akc), t=0 .. 20);

c) Rewrite the integrand as nonrecursive piecewise (which requires knowing the extend of the domain in advance, if you wanted to program it). Note use of lowercase `int`, not `Int`. This works because Maple can do this integration symbolically, not numerically. The result is smooth, and quite fast. But you may not be able to do this for other examples.

restart:
F:=t->piecewise(t>19,0,t>15,1,t>14,0,t>10,1,t>9,0,t>5,1,t>4,0,1):
plot(int(F,0..t),t=0..20);

d) Rewrite the integrand to not use piecewise. Below the integration is done out to point 17.9, with informational messages temporarily enabled to show the domain splitting.

restart:
f:=50*min(1,irem(ceil(x),5)):

infolevel[`evalf/int`]:=1:
evalf(Int(f, x=0..17.9));
evalf(%);
infolevel[`evalf/int`]:=0:

plot(Int(f, x=0..t), t=0..20);

2)

I don't understand part 2). Are both the operators supposed to be identical, and if so then could you explain what result you hope for?

I used Maple 13 for all the above, as I recall you specified that in a previous Question.

acer

Are you trying to get  as a name with an underscore, in 2D Math? Or do you really want it to appear as an indexed name when you enter it as 1D input?

If you want to use just  as the name with trailing underscore in 2D Math, and don't need to have that change (as a reference) when `i` changes, then you can use an atomic identifier. Perhaps the simplest way to input that is to use the keystrokes,

o m e g a Ctrl-Shift-minus i

where Ctrl-Shift-_ means pressing the Control, Shift, and minus keys together (presuming that your keyboard has underscore as the Shift-minus). This should produce something that typesets just like the indexed name but is actually a distinct name.

If, on the other hand, you just want it to look exactly like omega[i] as 1D input, then you can turn that into a name distinct from omega by wrapping it in single-left (name) quotes,

          `omega[i]`

You cannot assign to both the name omega[i] and also (separately) to the name omega, and have the first of those change as `i` changes value. The methods described above relate to using omega[i] as a distinct name (well, distinct from the name omega), which entails that the `i` in it be fixed as `i` and not be variable.

acer

There are a variety of ways to get this.

with(DocumentTools):

# variant 1
Do(%MathContainer0=cos(evalf(Do(%Slider0(value))*Pi/180))):

# variant 2
Do(%MathContainer0=Units:-Standard:-cos(evalf(Do(%Slider0(value))*Unit('deg')))):

# variant 3
with(Units:-Standard,cos):
Do(%MathContainer0=cos(evalf(Do(%Slider0(value))*Unit('deg')))):

# variant 4
with(Units:-Standard): # this breaks variants 1,2, and 3
SetProperty('MathContainer0', 'value',
            cos(evalf(Do(%Slider0(value))*Unit('deg')))):

I haven't looked yet at the code in DocumentTools:-Do, to figure out why loading all of Units:-Standard breaks for variants 1-3 above.

Never load Units:-Natural, just to get the effects that Units:-Standard also brings and just to avoid calling the Unit() constructor. It robs your session of far too may common single letters from being used easily in their non-unit sense.

acer

Your syntax for calling `getdata` was not correct.

But you'd also get a better result from a dedicated optimization routine. The plot structure will just give you a set of discrete values (the plot driver interpolates to get a smooth surface) so the maximal value from the plot structure won't be a very accurate answer to the question.

> P := plot3d(Qnum, Do=0.05..0.5, N=1..800,
>             axes=boxed, labels = ['Do','N', 'Q'],
>             view = [0.05 .. 0.5,1..800, 0..25000]):

> #plots:-display(P);

> MM := plottools:-getdata(P)[3]:

> max(select(type, MM, numeric)); # sieving out any Float(undefined) nonnumeric values

                              24553.6594800000

> Optimization:-Maximize(Qnum, Do=0.05..0.5, N=1..800);

Warning, undefined value encountered
    [24560.2479846227980, [Do = 0.322249987190719, N = 142.628104640423]]

acer

The 1D Array `colorarray` below has 3*N entries. For HSV coloring schemes this Array will be interpreted as a sequence of triples (ie. Hue, Saturation, and intensity Value). For RGB the triples will be interpreted as the Red, Green, and Blue values. The key is that the 3*N entries represent three entries for each of the N points, laid out in a contiguous sequence.

Download coloringbyz.mw

acer

The regular Optimization package can get this too (using Markiyan's objective formulation), and scaling (as Axel suggests),

Download Opt.mw

acer

Something (weird) is wrong with the 2D Math input of your assignment to EnergyDensity.

Correcting that, and and making a few adjustments to try and speed it up, gets a result of approximately 0.7528662772 for the root (taking about 4.5 sec on a fast i7).

fsolve2.mw

acer

It's quite a bad idea to unprotect and reassign to the global name D. For one thing, it would break any of Maple's own Library routines which relied on :-D correctly behaving as the differential operator.

Having said that, one can still unprotect (or rebind) and use D a module export. This allows one to use the name D interactively, while keeping the integrity of the global name :-D as the differential operator.

For example,

restart;
with(LinearAlgebra):

m := module() option package; export D; end module:

with(m):
unprotect(m:-D):

D := DiagonalMatrix([2, 3]);

                                      [2  0]
                                 D := [    ]
                                      [0  3]

D.D;

                                   [4  0]
                                   [    ]
                                   [0  9]

D^(2);

                                   [4  0]
                                   [    ]
                                   [0  9]

f := x->sin(x^2);

                                           / 2\
                              f := x -> sin\x /

:-D(f); # you can still use the differential operator, referenced as a global

                                          / 2\  
                                x -> 2 cos\x / x

dsolve( {diff(h(x),x,x)=x^2, h(0)=0, :-D(h)(0)=0} );

                                       1   4
                                h(x) = -- x 
                                       12   

acer

Is a numeric solution adequate, or did you really need an exact result?

Download Optimum3.mw

acer

Do your data points (well, their 2D independent components) form a regular grid in the x-y, longitude-latitude plane?

If yes, then you can use the plots:-surfdata command directly. You could also try to "improve" the surface by interpolation, like here.

If not, and your points are irregularly laid out in the x-y plane, then you can try to compute a surface using some scheme like triangulation (eg. Delaunay). See the links given by Christopher2222 in his Answer (and look for `Surfplot`).

acer

Replace X as desired.

with(LinearAlgebra):

M:=Matrix([[2,5,7,9,8],[8,5,2,6,8],[5,4,5,5,7],[1,4,9,6,6]]):

J:=Vector([3,5,4,9,6]):

#JJ:=DiagonalMatrix(J):  # not needed to form A

A:=Matrix(4,5, (i,j)->`if`( J[i]-J[j]-2=0, X ,M[i,j]/(J[i]-J[j]-2)) );

acer

restart:

L := proc(n::nonnegint)
     option remember;
        if n=0 then return 1;
        elif n=1 then return x;
        else
           return sort( x^(n)
                        - add(L(i-1)*int(x^(n)*L(i-1),x=-1..1)
                              /int(L(i-1)*L(i-1),x=-1..1),
                              i=1..n) );
        end if
     end proc:

seq(print(L(k)), k=0..9);

                                      1

                                      x

                                    2   1
                                   x  - -
                                        3

                                   3   3  
                                  x  - - x
                                       5  

                                4   6  2   3 
                               x  - - x  + --
                                    7      35

                               5   10  3   5   
                              x  - -- x  + -- x
                                   9       21  

                           6   15  4   5   2    5 
                          x  - -- x  + -- x  - ---
                               11      11      231

                          7   21  5   105  3   35   
                         x  - -- x  + --- x  - --- x
                              13      143      429  

                      8   28  6   14  4   28   2    7  
                     x  - -- x  + -- x  - --- x  + ----
                          15      13      143      1287

                     9   36  7   126  5   84   3    63   
                    x  - -- x  + --- x  - --- x  + ---- x
                         17      85       221      2431  

And, of course, you could also write the key fragment as,

                x^(n)
                - add(L(i)*int(x^(n)*L(i),x=-1..1)
                      /int(L(i)*L(i),x=-1..1),
                      i=0..n-1)

acer

You want the rope to move in the animation, yes?

Presumably you already have formulas for XCoordinate and YCoordinate which involve parameter t.

Pendulum:=proc(x,y)
   plots:-display(
      plots:-pointplot([x,y],color=blue, symbol=solidcircle,symbolsize=25),
      plottools:-line([0,1],[x,y],color=blue,thickness=3)
                  );
end proc:

A:=plots:-animate(Pendulum,[XCoordinate,YCoordinate],t=0..8,frames=10):

plots:-display([A],insequence=true);

acer

The command to load a package is lowercase `with`, not `With`.

Also, use the exponent %T to transpose, not T.

acer

It depends on whether you want just any solution, or all possible solutions.

restart:

eqn:=sin(2*x) = -sqrt(3)/2;

                                  1  (1/2)
                     sin(2 x) = - - 3     
                                  2       

sol1:= solve( eqn );

                               1   
                             - - Pi
                               6   

eval(eqn, x=sol1);

                      1  (1/2)     1  (1/2)
                    - - 3      = - - 3     
                      2            2       

sol2 := solve( eqn, AllSolutions ); # _Bx means 0 or 1, and _Zx means any integer

                     1      5                
                   - - Pi + - Pi _B1 + Pi _Z1
                     6      6                

about(indets(sol2,And(name,Non(constant))));

{_B1, _Z1}:
  is used in the following assumed objects
  [_B1] assumed OrProp(0,1)
  [_Z1] assumed integer

S1 := {seq(eval(sol2,[_B1=0,_Z1=i]),i=-2..2)}; # these are Pi apart

             /  13       7       1     5     11   \ 
            { - -- Pi, - - Pi, - - Pi, - Pi, -- Pi }
             \  6        6       6     6     6    / 

seq( eval(lhs(eqn),x=s), s in S1 );

     1  (1/2)    1  (1/2)    1  (1/2)    1  (1/2)    1  (1/2)
   - - 3     , - - 3     , - - 3     , - - 3     , - - 3     
     2           2           2           2           2       

S2 := {seq(eval(sol2,[_B1=1,_Z1=i]),i=-2..2)}; # these are Pi apart

               /  4       1     2     5     8   \ 
              { - - Pi, - - Pi, - Pi, - Pi, - Pi }
               \  3       3     3     3     3   / 

seq( eval(lhs(eqn),x=s), s in S2 );

     1  (1/2)    1  (1/2)    1  (1/2)    1  (1/2)    1  (1/2)
   - - 3     , - - 3     , - - 3     , - - 3     , - - 3     
     2           2           2           2           2       

acer