sand15

797 Reputation

11 Badges

9 years, 323 days

MaplePrimes Activity


These are answers submitted by sand15

As @tomleslie said, one BC is missing.

An example with a periodic BC for D[1](u) ( = diff(u(x, t), x) )

restart:
PDE  := diff(u(x, t), t) = -u(x, t)*diff(u(x, t), x) + 0.1 * diff(u(x, t), x$2):
IC   := u(x, 0) = sin(x):
BC   := u(0, t)=u(2*Pi, t), D[1](u)(0, t)=D[1](u)(2*Pi, t):
IBC  := {IC, BC}:

pds := pdsolve(PDE, IBC, numeric, time=t, range=0..2*Pi):  # just do as explained in the help page
plots:-display(
  seq(
    pds:-plot(
      t=tau
      , numpoints=50
      , color=ColorTools:-Color([rand()/10^12, rand()/10^12, rand()/10^12])
      , legend=('t'=nprintf("%1.2e", tau))
    )
    , tau in [seq](0..2, 0.25)
  )
)

Explanation

P := plot3d( [ b, ( b^2 ) / 4, -b / 2 ], b = -4 .. 4, c = -4 .. 4, thickness = 5, color = red  ); 

P is made of two elements

op(P) = MESH(...., COLOR(....)), THICKNESS(...)

When a MESH structure is plotted the "construction lines" of the mesh are too (default choice).
To avoid plotting them (see PLOT3D help page) you have to use the option STYLE(HIDDEN).

Then (no need to the THICKNESS option)

PLOT3D( MESH( op(op(P)), STYLE(HIDDEN) ) )

does the job.

Not that simple but this explains why the curve you got is black: the black boundaries of each element of the MESH "overwrite" the red on their interiors.


First thing to do is to load (ExcelTools:-Import ; look the corresponding help page) your data 

Data := ExcelTools:-Import(....):

Determine the number of lines:

N := numelems(M[.., 1])

To split Data into a train base and a test base :

  • define the fraction tau of data in the 
  • set NT := round(tau*N)  (or floor or ceil instead)
  • realise a random splitting:
    combinat:-randperm([$1..N]):
    TrainSet := %[1..NT]:
    TestSet  := %%[NT+1..N]
  • then:
    TrainData := Data[TrainSet];
    TestData  := Data[TestSet];

Apply Statistics with the fit method (linear or nonlinear) to TrainData..
Evaluate the quality of the fit on TestData.

Provide an Excel data file if you want more details.


The reason why you get a wong picture for your second test case is detailed in the mw file I'm nor able to load right now (it will be the case in a few hours).


Concerning this second test case (Maple 2020)

restart:
with(plots):
with(ComputationalGeometry):

# Feasibility domain

Q := inequal(
       Or(u^2+4*v^2 <= 4, And(u^2+v^2 < 4, (u+2)^2+2*(v-1)^2 <= 2))
       , u=-2..2, v=-2..2
       , nolines
     ):

# Extract the polygons that <inequal> has built

POL := select(has, [op(Q)], POLYGONS):
numelems(POL): 
                         2

# Mapping & plotting

F := plottools:-transform((u, v)  -> [u^3-v^2, u^2-v^3]):
display(F(Q)):

# Searching noon at 2pm
# The analysis of the poorness of the above plot shows it comes from the closure of the drawing
# of some of the elemmentary polygons POL contains (this will be clear in the mw file).
# All these concerned polygons having more than 3 sides, the basic idea is to triangulate 
# the feasibility domain using DalaunayTriangulation.

# first component of POL

pol  := op(1..-4, op(1, POL)):
pts  :=  `<,>`(pol)
tris := DelaunayTriangulation(pts):
display(map(x -> plottools:-polygon(pts[x], style='polygon'), tris), axes=none):
d1 := display(F(%)):

# second component of POL

pol  := op(1..-4, op(2, POL)):
pts  :=  `<,>`(pol)
tris := DelaunayTriangulation(pts):
display(map(x -> plottools:-polygon(pts[x], style='polygon'), tris), axes=none):
d2 := display(F(%)):

# drawings

display(`<,>`(d1, d2, display(d1, d2))):


The result is already very good but could probably be improved if the construction of Q was refined.
Use optionsimplicit=[grid=[50, 50]] in inequal(...) ... but beware of the computation time for finer grids !!!


File to come (either from this logname or from my non professional logname "mmcdara")
Here it is


DrawingIssue.mw
Maple

Mathematica


Unless I'm mistaken (it's early in the morning and I haven't get my coffee) the answer is 260

Reasoning.

  1. Choose a color for the right triangle : 5 choices.
  2. Choose a color for top and bottom triangles : 4 choices for each.
  3. For the left triangle :
    1. if colors of top and bottom triangles are the same, 4 choices are left.
    2. If they are not the same 3 choices are left.

The probability that the top and bottom triangles share the same color is 1/4.
Thus 
N = 5 * 4 * 4 * (1/4 * 4 + 3/4 * 3) = 5 * 4 * (4 + 9) = 20 * 13 = 260

Given the 2 elementary events R = "a Red ball is drawn" and W = "a White ball is drawn", this will give you all the compound events of an experiment in which 3 balls are drawn from the box.

restart:
alias(C=combinat):
Box := [R$3, W$3]:
(op @ C:-permute)~(C:-choose(Box, 3));


BTW, you code can be simply written this way

combinat:-powerset({$1..3})

 

As isolve seems (IMO) more powerful and versatile than msolve I have build a simple procedure _msolve  which uses syntax of msolve and the power of the isolve .
 

restart:
_msolve := proc(eqs, m)
  local v := convert(indets(eqs, name), list);
  local e := lhs~(eqs) =~ rhs~(eqs) +~ m *~ parse~(cat~(_K, [$1..numelems(v)]));
  local s := isolve(convert(e, set));
  map(u -> v =~ eval(v, u), [s]);
end proc:

_msolve(x^2=-1, 5);  # "your" problem
  [ [ x = 3 - 5 _Z1 ], [ x = 2 - 5 _Z1 ] ]

_msolve(2^i=3, 19);  # from msolve help page
  [ [ i = 13 - 18 _Z1 ] ]

_msolve({3*x-4*y=1, 7*x+y=2}, 19); # from msolve help page
  [ [ x = 15 + 57 _Z1 + 76 _Z2,  y = 11 + 38 _Z1 + 57 _Z2 ] ]

# A generalization of the previous system where 2 modulii are used:
  eqs := {3*x-4*y=1, 7*x+y=2}:
  m   := [7, 19]
_msolve(eqs, m);
  [ [ x = 15 + 37 _Z1 + 76 _Z2, y = 11 + 26 _Z1 + 57 _Z2 ] ]

# a specific instance
eval(%[], [_Z1=1, _Z2=2]);
   [ x = 204, y = 151 ]

# check
eval(irem~(lhs~(eqs), m), %);
   [ 1, 2 ]

 


Let P the list of (2D) points:

ch  := simplex:-convexhull(P, output=true):
pch := plottools:-getdata(ch)[1][1];    
n:=2:
while pch[n, 1] > pch[n_1, 1] do
  n := n+1:
end do:
NewtonPolygon := pch[1..n-1]


Explanation
pch contains the set of points (a matrix) ordered with respect to an increasing curvilinear abscissa, starting from the leftmost point and in counter clockwise sense.
The Newton polygon is the submatrix of pch wich contains all the points with increasing "x" (1st column of pch).
The loop is aimed to determine the nth point whose abscissa "x" is lower than the previous one

A minimalist answer you can decorate as you want.
I hope there is no typos for I do not have Maple right now.

(If you have any problems I will reply you when I'm home [login = mmcdara])

 

restart:
with(plots):
with(Maplets[Elements]):

f := proc(p1, p2)
local maplet:
maplet := Maplet(
  [
    [
      Plotter['PL1'](p1),
      Plotter['PL2'](p2)
    ]
  ]
):

Maplets[Display](maplet):
end proc:

A := dualaxisplot(sin(x), cos(x), x=0..2*Pi):
B := dualaxisplot(sin(x), cos(x), x=0..2*Pi, color=[green, blue]):

f(A, B);

 

To complete Carl's answer...

If you are as lazzy as I am, here is maybe a simple trick which could interest you:

Let R a convex region (!!!) :
R can be written by the list L= [P1, ..., PN] of its "corner" points, picked in the clockwise or counter clockwise sense.
For instance a list associated to R1 can be L =  [ [p_l, p_B], [epsilon, p_B], [epsilon,1], [p_l,1] ]

The Trick :

with(PolydrehalSets):
L :=  [ [p_l, p_B], [epsilon, p_B], [epsilon,1], [p_l,1] ]:
obj := PolydrehalSet(L):
rel := Relations(obj)  ; # rel contains the desired relations

PS 1: happily PolyhedraSet works well even on flat polyhedral, that is on polygons

PS2 : unhappily PolyhedraSet needs that each coordinate in L be a rational number
          If it's not so, you need to write L := convert~(L, rational) before calling Polyhedral Set.
          This puts some limitations on this trick, but it can remains usefull...

@Otavio 

From my experience the better way to do an action when clicking a button is to set its attribute onclick this way
onclick = Action( Evaluate( 'function' = 'Myproc', Argument('a1'), ..., Argument('an'))
If necessary an attribute 'target'='T' can be added in the definition of Evaluate(...)
Please see Maplets[Elements](Evaluate] for more details.

Here 'a1', ..., 'an' and 'T' are n elements the "master maplet" must contain (for instance n TextField elements referenced 'a1', ..., 'an' and a Plotter referenced 'T'.

When you clic on the desired button the procedure MyProc will do some operations of arbitrary complexity (just as any Maple's procedure) and run a "slave maplet" named M.

Example: The master maplet contains a TextField element named 'TF' whose field contains the value 10.
You want Myproc to open a maplet with N TextFields (named 'TF1', ..., 'TF10' below)

  1. associate to the opening button the action onclick = Action( Evaluate( 'function' = 'Myproc', Argument('TF'))
     
  2. Define MyProc:
    Myproc := proc(N)
    local M:
    M := Maplet(
        ......
       GridRow( seq( GridCell (TextField[sprintf("TF%d", n)](.....)),  n=1..N) )
       .......
    )
    Maplets[Display](M)
    end proc:
     
  3. I find interesting (a personal opinion), that the master maplet be encapsulated itself in a procedure:
    Main := proc():
    local MyProc, MainMaplet:
    ......

    ... write here the definition of Myproc...
    ... write here the definition of MainMaplet...
    Maplets[Display](MainMaplet);
    ......

    end proc

 

Unfortunately I do not have Maple on this machine: all I wrote hereis the result of a succession of back and forth between this machine and the one I have Maple on... I hope I missed nothing.

Let me know if you meet difficulties.

What I understand:

  1. You have a non stationnary random process Y indexed by X such that the conditional expectation of Y given X=x[i] (if it is so, it's not the "mean" as you wrote it) is given by the formula CondExp(Y | X=x[i]) = a*x[i]+b, where a and b are known values.
     
  2. You assume this random process is a Gaussian randon process, then a realization of Y given X=x[i] (or [Y(omega) | X=x[i]] for short) is CondExp(Y | X=x[i]) + R(omega) where R is a gaussian random variable of 0 mean and standard deviation sigma.
     
  3. You want to generate a realization "GP" of this random process


If I'm right, and if X is some row vector of size N, just do

S := Statistics:-Sample(Normal(0, sigma), N);
GP := a *~ X +~ b +~ S

Maube you could try this:


A := [ some list on numerical values ]:
B := sort(A, output=permutation); # returns the position of the elements of A sorted in increasing values

For more info pease look to the 'sort' help page.
You'll see its possible to sort in descending order too

PS : once you have found, from some algorithm, the mosition of the highest element of a list A, the position of the second highest element of A is just the position of the highest element of the list A from which youy have removed the highest element.
 

To complete Preben's answer

Knowing that int( 1/(s*sqrt(2*Pi)) * exp(-1/2*(x/s)^2, x=-infinity..+infinity ) = 1, it's immediate to obtain, "by hand",  the value of this integral.

If You want to integrate over a bow in R3, the result is also obvious.
Just remark that  exp(-u^2/a) = exp(-(1/2) * (u /s))^2)  with s = sqrt(a/2).
Then  int( exp(-u^2/a), u=P..Q)  = s*sqrt(2*Pi) * int( 1/(s*sqrt(2*Pi)) * exp(-(1/2) * (u /s))^2) , u=P..Q) 
                                                   = s*sqrt(2*Pi) * (Probability(U <= Q) -Probability(U <= P) 
where U is a Gaussian random variable of expectation 0 and standard deviation s.

The result involves the error function 

I supppose the data you import are containded in X?
In your case X is column vector (a matrix with its number of column equal to 1).

Doing regression implies at least 2 quantities: the independant variable(s) and the dependant one.
However, in your case, it seems there is only one.
It may be that you have stacked the independant and dependant variables in the structure X?
In this case reshape X as desired, for instance by writting XY := Matrix(2, 1440, convert(X, list))^+ to obtain a 1440x2 matrix with the first column containing the independant variable (for instance) and the second column the dependant one.

Other points:

  • you use x in the expression of sigma and X elsewhere.
  • the quotes are unnecessary: just write gor instance mu := add(X[i], i=1..numelems(X)) / numelems(X)
  • in the computation of sigma, reuse mu: sigma := add((X[i]-mu)^2, i=1..numelems(X)) / numelems(X) *****
  • ***** last: the denominator in the epression of sigma must be numelems(X)-1 and not  numelems(X), unless 2880 represents the entire population
     

To conclude, unlesse you really want to use the explicit expressions for the mean and the standard deviation, it's better to use the correspondind Maple's function (just refer to the package Statistics (for instance))

with(Statistics);
mu:= Mean(X);    # or m := Statistics:-Mean(X) if you do not want to load the whole Statistics package
sigma := StandardDeviation(X)

 

 

Last point, when you Import yout data with ExcelTools:-Import (see the help page) you can specity the value an empty cell must contain.
Suppose this value is some number A which cannot be a "correct"number of your sample (for instance A := I, the imaginary root).
If X is the name of the matrix you inported, you can then extract the submatrix whose the components are real (or different from A in a more general way).

 

1 2 3 4 5 Page 4 of 5