mmcdara

5259 Reputation

17 Badges

7 years, 111 days

MaplePrimes Activity


These are replies submitted by mmcdara

@Carl Love 

Thanks Carl.

I agree with you that the question is very vague.
However, I thought that an example, let's say a generic one, might encourage the OP to provide more details.
The problem with this type of question is that they can sometimes be like Russian dolls: as the exchanges progress, the initial question fades, new ones appear and the final problem can be very different from what we (or I) thought it was.
So let's wait and see.

Thanks again for your comment.


@AylLivedo
I don't know if you've done anything in Maple about your problem?
If so, you can download your mw file by clicking on the big green arrow in the menu bar.

What is this list you are talking about? What is its structure?

What is dataplot intended to do? Why do you mention it?

You are completely mistaken: Statistics:-Histogram is not at all aimed to plot binned data but to plot data after having itself binned them.
If you have data there is no need to bind them and seek next for some specific Maple procedure to plot these binned data: just apply Statistics:-Histogram to the original data. Here is an example
Binning.mw

"Doing the binning is trivial" ??? Absolutely not: if you want a binning which correctly approximates the underlying distribution the sample is drawn from, this is all but trivial (choice of the anchor, binwidth, method of splitting, ...)

Please provide an example of your data if what my file contains doesn't suit you.

@Traruh Synred 

I do not understand what is complicated and why you say "Might be easier just to count it "
It's pure nonsense as soon as N is large (see below).

restart

with(Statistics):

N := 10^6:
S := Sample(Normal(0, 1), N):

# How many points in the bin [-3.1, -2.87] ?

TallyInto(S, [-3.1 .. -2.87], bins=1);
rhs~(%);

[HFloat(-3.1) .. HFloat(-2.87) = 1074]

 

[1074]

(1)

# How many points over 3.5 ?

TallyInto(S, [3.5 .. max(S)], bins=1);
rhs~(%);

[HFloat(3.5) .. HFloat(5.3723343754921675) = 227]

 

[227]

(2)

# Another possibility (not to use for large values of N)

Select(s -> is(s > 2), S[1..1000])^+;
numelems(%);

Vector(4, {(1) = ` 1 .. 21 `*Vector[column], (2) = `Data Type: `*float[8], (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})

 

21

(3)

 

Download HowToCatchTheBins_2.mw

But if you are patient enough it is indeed a good idea to count yourself the 1074 points out of the million which lie in the bin [-3.1, -2.87]...


More than this your reply has no link with your initial question "How do I find the contents of a bin in a Maple histogram?"

@Rouben Rostamian  

The solution is not continuous. This appears cleary in what the OP writes  "The correct impulse response should be : y(t) = exp(-t/tau)*Heaviside(t)/tau" (a claim you can find in any textbook) and in your solution too.

From the theoritical solution provided by the OP, one can easily check it doesn't verify y(0)=0 because Heaviside(0) is undefined. Taking the limit (from the right for a causal system) at t=0+, doesn't solve the issue

limit(exp(-t/tau)/tau*Heaviside(t), t=0, right)
                               1 
                              ---
                              tau

Thus y(0)=0,  y(t) verifies the ode D(y)(t) + y(t)/tau = Dirac(t)/tau in the open half line (0, +oo), but y(t) is discontinuous at t=0.

Mybe what puzzles you in the solution I got is that absence of the Heaviside function?
I agree it's quite dusturbing, but this is a Maple problem:

laplace(Heaviside(t), t, p);
                               1
                               -
                               p
 invlaplace(%, p, t)
                               1

To find a solution which verifies both the ode and the initial condition on can use a smoothed ode where the Dirac term is replaced by a function f(t) such that int(f(t), t=0..+oo)=1 and f(t) is almost equal to 0 in the open interval (epsilon..+oo) where epsilon is a small strictly positive number.
For instance:

use Statistics in f := PDF('Gamma' (lambda, 2), t) end use:

This function has a peak close to 0, which is higher the smaller lambda is (lambda > 0)

Name it "an engineer trick" or "a physicist's tinkering" if you want, or consider this as a kind of "weak solution"
 

restart

de := D(y)(t) + y(t)/tau = Dirac(t)/tau;

(D(y))(t)+y(t)/tau = Dirac(t)/tau

(1)

smoothed_de := D(y)(t) + y(t)/tau = exp(-lambda*t)/tau;
Y := rhs(dsolve({smoothed_de, y(0)=0}))

(D(y))(t)+y(t)/tau = exp(-lambda*t)/tau

 

(-exp(-t*(lambda*tau-1)/tau)/(lambda*tau-1)+1/(lambda*tau-1))*exp(-t/tau)

(2)

use Statistics in f := PDF('Gamma' (lambda, 2), t) end use:
smoothed_de := D(y)(t) + y(t)/tau = f/tau:
Y := rhs(dsolve({smoothed_de, y(0)=0}))

Y := piecewise(t < 0, 0, 0 <= t, exp(-t/tau)*exp(t*(lambda-tau)/(lambda*tau))*t/(lambda^2-2*lambda*tau+tau^2)-exp(-t/tau)*tau*exp(t*(lambda-tau)/(lambda*tau))/(lambda^2-2*lambda*tau+tau^2)-exp(-t/tau)*tau*exp(t*(lambda-tau)/(lambda*tau))*t/(lambda*(lambda^2-2*lambda*tau+tau^2))+exp(-t/tau)*tau/(lambda^2-2*lambda*tau+tau^2))

(3)

col := () -> ColorTools:-Color([rand()/10^12, rand()/10^12, rand()/10^12]):

plots:-display(
  seq(
    plot(
      eval(Y, [tau=1, lambda=10^k])
      , t=1e-6..1
      , numpoints=10^3
      , color=col()
      , legend=typeset(lambda=10^k)
    )
    , k=-4..-1
  )
  #, axis[1]=[mode=log]
)

 

plots:-display(
  seq(
    plot(
      eval(Y, [tau=1, lambda=10^k, t=10^(z)])
      , z=-6..1
      , numpoints=10^3
      , color=col()
      , legend=typeset(lambda=10^k)
    )
    , k=-5..-1
  )
)

 

plots:-display(
  plot(
    eval(Y, [tau=1, lambda=10^(-6), t=10^(z)])
    , z=-6..1
    , numpoints=10^3
    , color=col()
    , legend=typeset('smoothed'(lambda=10^(-6)))
  ),
  plot(
    eval(exp(-t/tau)/tau, [tau=1, t=10^(z)])
    , z=-6..1
    , numpoints=10^3
    , color=col()
    , legend=typeset(unsmoothed)
  )
)

 

limit(Y, lambda=0, right) assuming t >= 0

1/(tau*exp(t/tau))

(4)

 

Download solution_of_the_smoothed_ode.mw

@dharr 
@rzarouf

Let R a matrix whose elements are iid Uniform(-1, 1) random variables.
The idea is to build the distribution (more precisely an approximation of it) of the eigenvalues of matrix A+epsilon*R.

Let L(R, epsilon)  the eigenvalues of matrix A+epsilon*R. and L0=L(R, 0) those of matrix A.
I use a pointwise representation of the distribution of L(R, epsilon) given a random sample of R, and a 2D-kernel reconstruction of this distribution.

Here is an interesting fact:

The kth eigenvalue  Lk(R, epsilon)  is not necessarily located "around"  L0, k but can be found in the neighbourhood of any other eigenvalue  L0, i.
Please see  (look to the first plot in the attached file and the topic Illustrating the "mixing" of the eigenvalues at the end of it.
So my 2D-kernel based representation, as well as the one produced by your pseudoplot procedure, are probably missing something, or are at least insufficient.

PseudoSpectrum_Another_POV.mw

@gabbbooooo 

From help(dsolve,numeric,BVP):

Thus the default initial value ('continuation'=name item) for lambda is 0 (with a default initial value to 1/100, see 'mincont'=numeric item) and has a final value of 1.
The ode system S must be rewritten S'(lambda) in such a way thar S'(1) = S.
Here I rewrote only the BCS at Y=0: when lambda=1 those BCS are equal to those of your own problem.

Using a continuation method requires an analysis of the reasons of the failures (a thing which often demand a knowledge of the problem to solve) ans is scarcely unique: different reformulations of the same problem do exists and I just use one arbitrarily.

About the shooting method: you will find a lot of Q&A on this site, just use the Search/Advanced Search  features.

@Muhammad Usman 

I just downloaded your last file and changed 1000 by 1e-3 (as I said in my answer): here is the result obtained with Maple 2015:

Download Num_PDE_1000_1e-3.mw

You can easily verify by yourself that Ididn't change the ICs 

@Rouben Rostamian  

"Calculating the tangent plane to a sphere is quite trivial.  Why bother with the geometry package for it?"

I think this thought is directed at me?
Here's my answer: I think you could say the same thing about many of the functions offered by geom and geom3d, and also about other packages. If you take this line of reasoning further, you might even wonder why these packages exist in the first place.
My point of view differs from (what I understand is) yours: as some people did the job to write these packages, or these functions, why shouldn't I use them

Of course you could argue that they may not be computationally efficient (which I do not know) and that ad hoc coding might be preferable. And this would be a good reason not to use geom/geom3d.
Personally I like to use them because they are very close, IMO, to the geometry reasoning (by opposition to the analytic method you proposed).The main criticism I would level against them is that those packages have some weakness in taking assumptions into account, for instance, here, to build the tangent plane at a formal point [u, v, w]:

point(P, [u, v, w]):
TangentPlane(T, P, S) assuming IsOnObject(P, S);
IsOnObject: hint: unable to determine if u^2+v^2+w^2-2*u-6*v-10*w-134 is zero
IsOnObject: hint: unable to determine if u^2+v^2+w^2-2*u-6*v-10*w-134 is zero
Error, (in sprintf) insufficient parameters for algebraic format
TangentPlane(T, P, S) assuming u^2+v^2+w^2-2*u-6*v-10*w-134=0
IsOnObject: hint: unable to determine if u^2+v^2+w^2-2*u-6*v-10*w-134 is zero
Error, (in sprintf) insufficient parameters for algebraic format

The unnamed theorem you refer to is likely due to Lotka
https://www.pnas.org/doi/pdf/10.1073/pnas.18.2.172
(page 4/7

... but "it an be shown" is a quite vague statement and I wasn't capable to get Lotka's complete article)

In case you would be interested in pursuit problems, here is a Maple code:
https://www.hsu.edu/site/assets/files/4579/2006-7afpursuit.pdf

@acer 

forget it, I'll delete my comment anyway

 

Deleted by the autor

You assign xi to 0.65  BEFORE differentiating wrt xi

@vv 

Thank you vv

@Kitonum 

Thanks for your feedback.
Concerning point 2 "But how will you solve a similar problem if the number  n  is large", this was the meaning of my my writting "regardless of its computational efficiency".
I was aware that that "my" solution will present serious problems of computational efficiency when n is large.

@Kitonum 

No offense intended, but I think your program has two shortcomings:

  • firstly, it seems too complex,
  • secondly, it's not a good idea (IMO) to specify a range as the second argument of Partition.

For instance you write 

n := 9:
Partition(n, 0..n, {2,3,6,9});

because you know that there is a one-element partition of 9, but why do you write

Partition(8, 2..8, {2,3,6,9});

Why?
In fact a correct range for partitionning any number n should have been

Partition(n, 0..n, {2,3,6,9});

, which proves the specification of this range is useless.
Unfortunately Partition requires the first operator of the range to be a strictly positive integer

Partition(8, 0..8, L);
Error, (in Partition) invalid input: k_Partition expects its 2nd argument, k, to be of type posint, but received 0


Note that Partition does not always return the correct answer:

Partition(8, 1..8, L);
             [[8], [2, 6], [2, 3, 3], [2, 2, 2, 2]]



To answer the OP's question this simpler code below seems[*]  enough (regardless of its computational efficiency) and doesn't require the specification of any range.

f := proc(n::posint, L::set)
  local P := combinat:-partition(n):
  select((x -> is ({x[]} subset L)), P, L)
end proc:

L := {2, 3, 6, 9}:
f(8, L);
               [[2, 2, 2, 2], [2, 3, 3], [2, 6]]
f(9, L);
             [[2, 2, 2, 3], [3, 3, 3], [3, 6], [9]]
f(3, L);
                             [[3]]
f(1, L)
                               []

[*]  we both implicitely assume then the number to be partionned was a positive integer. What does the OP want if this number is negative and the set of number contains negative numbers, for instance:

n := -1:
L := {-2, 1}:
# then n writes
n := L[1] + L[2]



Your procedure  Partition may be used in a lot of ways and is then more general than the simple code above.
Nevertheless I would suggest you to modify Partition in order to eliminate this unnecessary range argument  range argument and correct the "Partition(8, 1..8, L) error"

I hope I wasn't too blunt in my comment.

1 2 3 4 5 6 7 Last Page 1 of 113