5635 Reputation

17 Badges

7 years, 189 days

MaplePrimes Activity

These are replies submitted by mmcdara


I was kidding, don't take it at face value.

What I wanted to say is that as soon as you have more equations than unknowns the existence of a solution is hazardous, unless those equations are linked in some way such that their number reduces after simple algebraic transformations.

In the three attach files you will find:

  • a simple proof that T1 has solutions and that your original system can be reduced to a system of 2 equations in 2 unknowns (file , file proves there exist 4 solutions one can easily obtain by hand)
  • a simple proof that T2 doesn't have any solution even if you add d as an unknown.  and

I recognize in your question a divided difference approximation of a diffusion equation with diffusion coefficient 1+a*gamma*alpha(x).
Am I right?

If it is so, then your discrete approximation of the diffusion scheme is not consistent

This can easily seen by putting alpha(x) = Constant.
The proof is given at the end of the attached file whose first part is dedicated to the construction of a consistent approximation of the diffusion term

diff(alpha(x)*diff(u(x, t), x), x)



Consistent divided difference approximation of a diffusion term

diffusion := Diff(alpha(x)*Diff(u(x, t), x), x)

Diff(alpha(x)*(Diff(u(x, t), x)), x)


# Introducing the flux phi(x, t) = alpha(x)*Diff(u(x, t), x),
# rewrite diffusion this way

flux_div := Diff(phi(x, t), x)

Diff(phi(x, t), x)


# Assuming the spatial mesh is uniform with (constant) step
# h = x[i+1]-x[i] whatever the value of i (let's put apart the
# problem of boundary conditions).
# At point x=x[i] a centered divided difference approximation
# of flux_div is

cdd_flux_diff := i -> (phi(x[i]+h/2, t) - phi(x[i]-h/2, t))/h

proc (i) options operator, arrow; (phi(x[i]+(1/2)*h, t)-phi(x[i]-(1/2)*h, t))/h end proc


# A centered divided difference approximation of phi(x[i]+(1/2)*h, t) at
# point x=x[i]+h/2 is

phi(x[i]+(1/2)*h, t) = alpha(x[i]+(1/2)*h)*(u(x[i]+h, t) - u(x[i], t))/h

phi(x[i]+(1/2)*h, t) = alpha(x[i]+(1/2)*h)*(u(x[i]+h, t)-u(x[i], t))/h


# A centered divided difference approximation of phi(x[i]-(1/2)*h, t) at
# point x=x[i]-h/2 is

phi(x[i]-(1/2)*h, t) = alpha(x[i]-(1/2)*h)*(u(x[i], t) - u(x[i]-h, t))/h

phi(x[i]-(1/2)*h, t) = alpha(x[i]-(1/2)*h)*(u(x[i], t)-u(x[i]-h, t))/h


# Combining relations 5, 6 and 7 gives

cdd_diffusion := eval(cdd_flux_diff(i), [(4), (5)])

(alpha(x[i]+(1/2)*h)*(u(x[i]+h, t)-u(x[i], t))/h-alpha(x[i]-(1/2)*h)*(u(x[i], t)-u(x[i]-h, t))/h)/h


# Write cdd_diffusion in a more convenient form

select(has, indets(cdd_diffusion, function), t);

cdd_diffusion := collect(cdd_diffusion, %);

{u(x[i], t), u(x[i]-h, t), u(x[i]+h, t)}


(-alpha(x[i]+(1/2)*h)/h-alpha(x[i]-(1/2)*h)/h)*u(x[i], t)/h+alpha(x[i]-(1/2)*h)*u(x[i]-h, t)/h^2+alpha(x[i]+(1/2)*h)*u(x[i]+h, t)/h^2


# Assuming alpha(x) is continuous within each cell x[i]..x[i+1] whatever the
# value of i, simple approximations of alpha at "mid" points x[i]-(1/2)*h and
# then are

alpha(x[i]-(1/2)*h) = (alpha(x[i]-h) + alpha(x[i]))/2,
alpha(x[i]+(1/2)*h) = (alpha(x[i]) + alpha(x[i]+h))/2;

alpha(x[i]-(1/2)*h) = (1/2)*alpha(x[i]-h)+(1/2)*alpha(x[i]), alpha(x[i]+(1/2)*h) = (1/2)*alpha(x[i])+(1/2)*alpha(x[i]+h)


# Plug these relations into (7)

cdd_diffusion := map(simplify, eval(cdd_diffusion, [(8)]))

-(1/2)*(2*alpha(x[i])+alpha(x[i]+h)+alpha(x[i]-h))*u(x[i], t)/h^2+(1/2)*(alpha(x[i]-h)+alpha(x[i]))*u(x[i]-h, t)/h^2+(1/2)*(alpha(x[i])+alpha(x[i]+h))*u(x[i]+h, t)/h^2


# You wrote

yours := -(
-(1/2)*(2*alpha(x[i])+alpha(x[i]+h)+alpha(x[i]-h))*u(x[i], t)/h^2+(1/2)*(alpha(x[i]+h)+alpha(x[i]-h))*u(x[i]-h, t)/h^2-(1/2)*(alpha(x[i]+h)-alpha(x[i]-h))*u(x[i]+h, t)/h^2

(1/2)*(2*alpha(x[i])+alpha(x[i]+h)+alpha(x[i]-h))*u(x[i], t)/h^2-(1/2)*(alpha(x[i]+h)+alpha(x[i]-h))*u(x[i]-h, t)/h^2+(1/2)*(alpha(x[i]+h)-alpha(x[i]-h))*u(x[i]+h, t)/h^2


# assuming alpha is a constant C gives

eval(yours, alpha=(z -> C))

2*C*u(x[i], t)/h^2-C*u(x[i]-h, t)/h^2


# as expression (9) gives
# (you should recognize here the correct estimation of C*diff(u(x, t), x$2)
eval((9), alpha=(z -> C))

-2*C*u(x[i], t)/h^2+C*u(x[i]-h, t)/h^2+C*u(x[i]+h, t)/h^2


# Proof that (12) approximates C*diff(u(x, t), x$2) = C*(D[1, 1](u))(x[i], t)
# at point x=x[i] with error O(h^2)

taylor((12), h, 3)

series(C*(D[1, 1](u))(x[i], t)+O(h^1),h,1)


# What does your expression (11) approximates?

taylor((11), h, 3)

Error, does not have a taylor expansion, try series()


series((11), h, 3)

series((C*u(x[i], t))/h^2+(C*(D[1](u))(x[i], t))/h-(1/2)*C*(D[1, 1](u))(x[i], t)+O(h^1),h,1)


# Obviously formula (10) is not a consistent approximation of
# the continuous diffusion operator (1)


Add on



# Mass conservation within cell x=x[i]-h/2..x[i]+h/2.
# For any function psi(x) "regular enough" one gets this weak expression
# of mass conservation

MC := Int(eval(diffusion, [u=((x, t) -> U(x))])*psi(x), x=X[i]-h/2..X[i]+h/2) = 0

Int((Diff(alpha(x)*(Diff(U(x), x)), x))*psi(x), x = X[i]-(1/2)*h .. X[i]+(1/2)*h) = 0


# Integrating by parts gives

MC := Parts(MC, psi(x)) assuming h > 0

alpha(X[i]+(1/2)*h)*(D(U))(X[i]+(1/2)*h)*psi(X[i]+(1/2)*h)-alpha(X[i]-(1/2)*h)*(D(U))(X[i]-(1/2)*h)*psi(X[i]-(1/2)*h)-(Int(alpha(x)*(Diff(U(x), x))*(diff(psi(x), x)), x = X[i]-(1/2)*h .. X[i]+(1/2)*h)) = 0


# A particular choice of psi(x) is psi(x)=constant which, without loss of
# generality, can be chosen equal to 1.
# This gives:

value(eval(MC, psi=(x -> 1)))

alpha(X[i]+(1/2)*h)*(D(U))(X[i]+(1/2)*h)-alpha(X[i]-(1/2)*h)*(D(U))(X[i]-(1/2)*h) = 0


# This is basically the equation cdd_flux_diff(i)=0.
# Now you can use any formula you want to approximate the 4 terms
# relation (17) contains, provided the final result is a consistent
# approximation of the fluxes at walls x=X[i]-h/2 and X[i]+h/2.
# What is used here is a staggered mesh where cells are X[i]-h/2..X[i]+h/2,
# U(x) is expressed at the center of the cells and fluxes at their walls.




The discretization of the term diff(u(x, t), t) is not consistent neither: a term of the form u(x[i], t[n+p]), with p=+1 (explicit scheme), p=-1 (fully implicit scheme), or any other value of p, is missing.

Before thinking to build the matrices and vectors you have in mind, I advice you to verify your discretization scheme.

As I can only use Maple 2015 at the moment and many improvements have since been made to pdsolve, I can just give you a simple advice: browse this site to find if Maple (site edited up to Maple 2021) is capable to handle your problem.

I wasn't capable to find a test case corresponding to your problem but, on the other side, I'm not sure this problem is well posed (I have to get back to my old notes to check that).
I seem to remember that the boundary conditions have to be set at both ends of the bar and, even in this case, you do not necesseraly have a well-posed problem.


I use to use only worksheet mode with 1D inputs.
Here is the file written with Maple 2015

Strangely your reply is not correctly displayed with Safari (see below) as it is with Firefox.
I believe this is the first time I notice this.

The first one is a classical advice to all newbye: download your Maple worksheet by using the green up-arrow in the menu bar.

Given the complexity of your pdes there is no chance at all that you can get a formal solution.
To expect a numerical approximation of the solution you will have to give numerical values to all the sonstants and provide initial and boundary conditions.
Actually, I'm not even sure that Maple can solve numerically this system.

The last advidce: write command lines in separate blocks to understand where the error come from.
If you write the commands in a single block all the errors are displayed stacked at the end of the block.

Take inspiration from this


alias(U[1]=U[1](t, x, y, z)):
alias(U[2]=U[2](t, x, y, z)):
alias(U[3]=U[3](t, x, y, z)):
alias(phi=phi(t, x, y, z)):

PDE1 := c[1, 1]*diff(U[1], y, y) + c[1, 2]*diff(U[2], x, y) + c[1, 3]*diff(U[3], x, z) + e[3, 1]*diff(phi, x, z) + c[6, 6]*diff(U[2], x, y) + c[6, 6]*diff(U[1], y, y) + c[4, 4]*diff(U[3], x, z) + c[4, 4]*diff(U[1], z, z) + e[1, 5]*diff(phi, x, y) = rho*diff(U[1], t, t):


PDE1 := simplify(PDE1, size);

(c[1, 1]+c[6, 6])*(diff(diff(U[1], y), y))+(c[1, 2]+c[6, 6])*(diff(diff(U[2], x), y))+(c[1, 3]+c[4, 4])*(diff(diff(U[3], x), z))+e[1, 5]*(diff(diff(phi, x), y))+e[3, 1]*(diff(diff(phi, x), z))+c[4, 4]*(diff(diff(U[1], z), z)) = rho*(diff(diff(U[1], t), t))


# Proceed this way for the other pdes.



@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.

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

"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).



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);

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




# How many points over 3.5 ?

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

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




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

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

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






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)

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);
 invlaplace(%, p, t)

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"


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

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


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




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))


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

      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]


      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


    eval(Y, [tau=1, lambda=10^(-6), t=10^(z)])
    , z=-6..1
    , numpoints=10^3
    , color=col()
    , legend=typeset('smoothed'(lambda=10^(-6)))
    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






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.


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:


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
(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:


forget it, I'll delete my comment anyway

5 6 7 8 9 10 11 Last Page 7 of 119