acer

32490 Reputation

29 Badges

20 years, 8 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are answers submitted by acer

Here is something you might experiment with, simply using indexing,

M := Matrix([[1,2],[-2,3],[3,18],[-5,6]]);

                                     [ 1   2]
                                     [      ]
                                     [-2   3]
                                M := [      ]
                                     [ 3  18]
                                     [      ]
                                     [-5   6]

m,n := op(1,M);

                                m, n := 4, 2

rowlist:=[seq(`if`(not M[i,1]<0, i, NULL), i=1..m)];

                              rowlist := [1, 3]

M[rowlist,1..n];

                                   [1   2]
                                   [     ]
                                   [3  18]

acer

This below is a simple minded approach. It's likely slower than necessary for very some long lists that in principle could be failed early, since it builds all the conditions before checking any of them.

palin:=proc(T::list)
  evalb(`and`(seq(T[i]=T[nops(T)-i+1],i=1..nops(T)/2)));
end proc:

palin( [a,b,c,d,c,b,a] ); # odd number of entries

                              true

palin( [a,b,c,c,b,a] ); # even number of entries

                              true

palin( [a,b,c,d,b,a] );

                             false

palin( [17,[a,b],0.3,0.3,[a,b],17] );

                              true

palin( [P] );

                              true

palin( [] ); # add special case for this, if undesirable

                              true

acer

Here's something that you could try out. The basic idea is that the assignment is not distributed through the translated object. The hope is that the Maple `piecewise` could be used in usual ways (here, assignment).

Sorry if there are any syntax errors to fix, I am not in front of a C compiler this minute.

I reversed the conditions in your example, just to make it perhaps a little more interesting.

restart:

with(CodeGeneration):
with(LanguageDefinition):
    
    LanguageDefinition:-Define("NewC", extend="C",
        AddFunction("piecewise", anything::numeric,
            proc()
                local i;
                Printer:-Print("( (",_passed[1],") ? ",_passed[2]);
                for i from 3 to _npassed-2 by 2 do
                    Printer:-Print(" : (",_passed[i],") ? ",_passed[i+1]);
                end do;
                Printer:-Print(" : ",_passed[_npassed]," ) ");
            end proc,
        numeric=double)
    );

myp:=proc(x::numeric) local result::numeric;
        result := piecewise(x>3,3*x,x>2,2*x,x>1,1*x,0);
    end proc:

Translate(myp, language="NewC");

double myp (double x)
{
  double result;
  result = ( (0.3e1 < x) ? 0.3e1 * x : (0.2e1 < x) ? 0.2e1 * x : (0.1e1 < x) ? x : 0 ) ;
  return(result);
}

acer

Try the Maple 15 User Manual which is available from here.

acer

You could change settings so that by default your get the 1D instead of the 2D Math input mode while working in the Standard GUI. See here. That would make "whole equation stays inline," like in the Classic GUI.

You can also set another option so that by default you get Worksheets instead of Documents in the Standard GUI. (ie. from the top menubar, Tools->Options->Interface tabm then "Default format for new worksheets". That makes it behave even more like Classic, for entering input.

acer

You can create a procedure that generates the EOM for any value of elenum.

restart:
makeEOM:=proc(elenum)
local i, node, nodes, L, E0, E0M;
   nodes:=elenum*2+1:
   L:=evalf((Pi*2)/(elenum*2)):
   E0:=Matrix([[L/3,2*L/3,L/3],[2*L/3,11*L/15,2*L/3],[L/11,2*L/3,L/110]]):
   E0M:=Matrix(nodes);
   for i from 1 to elenum do
      node:=(2*(i-1))+1;
      E0M[node..node+2,node..node+2]:=E0M[node..node+2,node..node+2]+E0;
   end do:
   E0M:
end proc:

Now you can create then at will,

interface(rtablesize=infinity):
makeEOM(6);

And if you want to you can create them for all the values in that Array, and store them.

V := Array([2,4,6,8,16]):
for en in V do
  eom[en]:=makeEOM(en); # store them for later use
end do:

eom[4];

And you can also access them by position of the value in V,

V[1];

eom[V[1]];

Procedures can improve the flow of your code.

acer

If this is a floating-point problem, and you want a double-precision computation, then you are in luck. Putting shape=band[1,1] and datatype=float[8] on your Matrix and your RHS (Vector, or Matrix) will greatly improve the performance by virtue of using a dedicated solver from CLAPACK for band shape & storage. (There also exist dedicated solvers for tridiagonal Matrices -- which usually means specifically the main, 1st super-, and 1st sub-diagonal in the numerical linear algebra discipline -- but the band solver can do relatively well here too, as illustrated below.)

First, let's look at the case where the datatype=float[8] Matrix A has full rectangular storage.

> restart:
> N:=5000:
> with(LinearAlgebra):                                     

> Arect:=RandomMatrix(N,outputoptions=[datatype=float[8]]):

> kernelopts(bytesalloc); # total memory allocated so far  
                                   201813920

> B:=RandomVector(N,outputoptions=[datatype=float[8]]):

> X := CodeTools:-Usage( LinearSolve(Arect,B) ):
memory used=191.68MiB, alloc change=191.34MiB, cpu time=12.62s, real time=12.67s

> kernelopts(bytesalloc); # total memory allocated so far  
                                   402448408

> Norm(Arect.X-B);    
                                                  -8
                           0.971647295955335721 10

Nor let's compare, but this time Matrix A will also have shape=band[1,1]. Notice that the total memory allocation to store A and to get the computation done is much less, the computation is much faster, and the result has better forward error upon resubstitution.

> restart:                                                 
> N:=5000:                                                 
> with(LinearAlgebra):                                     

> A:=RandomMatrix(N,outputoptions=[shape=band[1,1],datatype=float[8]]):
> kernelopts(bytesalloc); # total memory allocated so far              
                                    2096768

> B:=RandomVector(N,outputoptions=[datatype=float[8]]):                

> X := CodeTools:-Usage( LinearSolve(A,B) ):                           
memory used=1.08MiB, alloc change=0.87MiB, cpu time=11.00ms, real time=32.00ms

> kernelopts(bytesalloc); # total memory allocated so far
                                    3145152

> Norm(A.X-B);                                                         
                                                 -10
                          0.436557456851005554 10

In essence, the computational complexity of the problem has been reduced by a power of N, affecting both memory storage and computation time. That means that the relative benefit of using the specialized shape & storage will increase (ie. get even better, in a relative sense) as the size of the problem increases.

Now, you will be constructing your Matrix A using some process other than calling RandomMatrix, of course. If your code creates A up front, and then populates it in some manner, then make sure that the creating has those options. Eg,

A := Matrix(m,m,datatype=float[8],shape=band[1,1]);

Let us know, if you have problems implementing it is this way.

acer

See the help-page for `is` and related property checking,

is( 5/sqrt(x+4) > 0 ) assuming x>1;

                                       true

acer

restart:

for k to 5 do
   f[k] := unapply(piecewise(x > 0, k, 0),x);
end do:

f[3](1);
f[3](-1);
f[2](1);
f[2](-1);

restart:

f:=proc(x)
      if type(procname,indexed) then
         piecewise(x>0,op(1,procname),0);
      end if:
end proc:

f[3](1);
f[3](-1);
f[2](1);
f[2](-1);

restart:

f:=proc(x,k)
      piecewise(x>0,k,0);
end proc:

f(1, 3);
f(-1, 3);
f(1, 2);
f(-1, 2);

acer

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 omega[i] 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 omega[i] 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.

 

restart:
N:=100:

# You data will have been obtained in some other manner than this,
# which is just randomly generated.
#
x:=LinearAlgebra:-RandomVector(N,generator=0.0..1.0,
                               outputoptions=[datatype=float[8]]):

y:=LinearAlgebra:-RandomVector(N,generator=0.0..1.0,
                               outputoptions=[datatype=float[8]]):

z:=LinearAlgebra:-RandomVector(N,generator=0.05..40.0,
                               outputoptions=[datatype=float[8]]):

# Your data's minimal and maximal heights may be different from this.
#
zmax,zmin:=max(z),min(z);

39.5198810263072318, 1.19553239094103869

# An Array of color values. This particular scheme is similar to shading=zhue
# for 3D plots. But you could also create your own R, G, and B functions
# depending on height z-value.
#
colorarray:=Array(1..3*N,[seq([0.92-abs(z[i]-zmin)/(zmax-zmin),1,1][],i=1..N)],
                  datatype=float[8], order=C_order):

# The 2D pointplot, colored.
#
P:=plots:-pointplot(Matrix([x,y],datatype=float[8]),axes=box,
                    symbol=solidcircle,symbolsize=10):
op(0,P)(op(P),COLOR(HSV,colorarray));

# Now conjoin the data Vectors for a 3D plot.
#
M:=Matrix([x,y,z],datatype=float[8]):

# The shading=none option adds only minimally to the cost of creating this
# initial structure, which of course is not yet visually meaningful.
#
#orien:=[45,75,0]:
orien:=[-90,0,0]:
P:=plots:-pointplot3d(M,axes=box,symbol=solidcircle,symbolsize=8,
                      shading=none,orientation=orien,view=0..zmax):

# Now add out coloring scheme.
#
newP:=subsindets(P,specfunc(anything,SHADING),
           t->COLOR(HSV,colorarray)):
newP;

 

 

Download coloringbyz.mw

acer

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

 

restart:

f:=x->646.0548255+(309.8026777+568.8289152*I)*((.8781989023+.4782956073*I)
*ln(1.+(-.8781989023+I*((-1)*.4782956073))*exp(10000000*((-1)*9.013850411)
*Dif+(61.46924721*I)*x))+(.8781989023+.4782956073*I)
*ln(1.+(-.8781989023+I*((-1)*.4782956073))*exp(10000000*((-1)*9.013850411)
*Dif+I*((-1)*61.46924721)*x))+(-.8781989023+I*((-1)*.4782956072))
*ln((.8781989023+I*((-1)*.4782956073))*(((-1)*1.)*exp(10000000*((-1)*9.013850411)
*Dif+(61.46924721*I)*x)+.8781989023+.4782956073*I))+(-.8781989023+I*((-1)*.4782956072))
*ln((.8781989023+I*((-1)*.4782956073))*(((-1)*1.)*exp(10000000*((-1)*9.013850411)
*Dif+I*((-1)*61.46924721)*x)+.8781989023+.4782956073*I))):

x_values:=[0.4056604928e-2, 0.1242487078e-1, 0.2106033816e-1,
           0.2965936896e-1, 0.3814909006e-1, 0.4673597534e-1]:

y_values:=[3274.140334, 746.1905199, 2.356309641, 0, 0, 0]:

obj:=eval(add(abs(evalc(f(x_values[j]))-y_values[j])^2,j=1..nops(x_values)),Dif=scaledDif*1e-9):

sol:=CodeTools:-Usage( Optimization:-Minimize(obj,scaledDif=0..10) );

memory used=11.00MiB, alloc change=4.50MiB, cpu time=172.00ms, real time=163.00ms

          sol := [98315.0959828365, [scaledDif = 1.55135921397022]]

plot(obj,scaledDif=0..10,view=0..2e6);

isolate(eval(sol[2],scaledDif=Dif*1e9)[],Dif);

                                                  -9
                         Dif = 1.55135921397022 10  

 

 

Download Opt.mw

acer

First 268 269 270 271 272 273 274 Last Page 270 of 337