Carl Love

Carl Love

26658 Reputation

25 Badges

11 years, 228 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are answers submitted by Carl Love

The problem is that EmpiricalDistribution is a discrete distribution, yet the sample is drawn from a continuous distribution. Something almost analogous to EmpiricalDistribution that works for continuous distributions is KernelDensity. Specifically, KernelDensity returns a probability density function (PDF) deduced from a sample drawn from a continuous distribution. Once you have the PDF, you can create a custom Distribution, and then a RandomVariable from that.

restart;

with(Statistics):

N:= RandomVariable(Normal(10,10)):

P:= Sample(N, 2^9):

A much larger sample would work better, but I'm sticking close to your 500.

FrequencyPlot(P);

F:= KernelDensity(P, method= exact):

plot(F, -20..40);

EmpDist:= Distribution(PDF= F):

X:= RandomVariable(EmpDist):

The DensityPlot is just the same PDF.

DensityPlot(X, range= -20..40);

Compare with the underlying distribution:

DensityPlot(N);

 

 

Download KernelDensity.mw

Answering your first question, I strongly suspect that in the code before the code that you show that you assigned a plot to the variable a. Doing a restart should fix that.

I believe that the following will work in Maple 12. Let me know if it does not.

restart;

a:= 0:  b:= 0.8:  step:= 0.2:

IVP:= {diff(y(t),t)=y(t)-2, y(a)=1}, {y(t)}:

soln1:= dsolve(
   IVP, numeric, method= classical[foreuler],
   stepsize= step,
   output= Array([seq](T, T= a..b, step))
)[2][1];

soln1 := Matrix(5, 2, {(1, 1) = 0., (1, 2) = 1., (2, 1) = .2, (2, 2) = .800000000000000, (3, 1) = .4, (3, 2) = .560000000000000, (4, 1) = .6, (4, 2) = .272000000000000, (5, 1) = .8, (5, 2) = -0.735999999999998e-1})

soln2:= dsolve(IVP, numeric):  #high-accuracy numerical solution

soln3:= dsolve(IVP);  #exact solution

y(t) = 2-exp(t)

P1:= plot(
   [soln1 $ 2],
   style= [point, line],
   symbol= cross, symbolsize= 16,
   thickness= 0,
   color= red
):

P2:= plots:-odeplot(soln2, t= a..b, color= blue):

P3:= plot(eval(y(t), soln3), t= a..b, color= green):

plots:-display([P1,P2,P3]);

 

 

Download Euler.mw

Here (attached) is your code with my (minor) modifications to make it run in hardware floating point. I ran it up to size 512 in less than a minute, and got reasonable accuracy (norm of the residuals 1e-16 seems reasonable to me). Let me know if this fulfills your needs.

Markov_a2.mw

In your program, what is the purpose of the following lines?

if type(y, rational) = false then
   y:= evalf(y);
   y:= convert(y, rational)
end if;

In other words, why are you using exact rational arithmetic, even after floating point approximation?

I ran your code for m=50, and the answers have more than 10000 digits in the numerator and denominator. It could be over a million digits for m=300.

Maple's sign is the sign of the leading coefficient of a polynomial, so sign(x) is always 1. Maple's signum is fine if you are dealing with strictly real expressions; for complex values, it may be counterintuitive (I'm not saying that it's wrong---it just may not be what you expect). An alternative is csgn. See ?signum, ?csgn, and ?sign.

Your IVP can be solved symbolically. The solution is x(t) = (1+exp(-t))/2. Thus .5 < x(t) < 1 for all t > 0. Thus both piecewises and the abs make no difference at all and the whole ODE reduces to

diff(x(t), t) = .5 - x(t).

When you make assumptions on a variable, a copy of that variable is made. If you then assign to that same variable, then you'll have the original and the copy being used in different ways, and things get confusing. The takeaway advice for the inexperienced user is Don't assign to assumed variables.

Your last two lines were

  • a:=solve(z1,a);
  • z2;

You should change those to

  • soln:= solve(z1, {a});
  • eval(z2, soln);

When used in this way, eval has most of the desired effect of an assignment statement without the side-effect of changing the global state. You can of course assign the result of eval to another variable.

Better than what you asked for---the solutions for a static list---we can make a "blackbox" solution that provides the x for any given y. This is what symbolic computation systems excel at.

restart;

eqn:= x = y*(25+2*x)^(2/3)/(25*(.566*(25*x)^(2/3)));

x = 0.282685512367491e-2*y*(25+2*x)^(2/3)*25^(1/3)/x^(2/3)

X:= Y-> fsolve(eval(eqn, y= Y), x);

proc (Y) options operator, arrow; fsolve(eval(eqn, y = Y), x) end proc

X(1);

.205293886670198

X(365);

8.68020335645564

plot(X, 1..365, labels= [y,x]);

 

 

Download blackbox.mw

A solution is to ?overload the procedure `diff/Int`.

I coded the overload of `diff/Int` solution. It was actually much easier than I was anticipating, just 11 lines of code. This solution illustrates a very generally applicable technique, one that can be used to modify the behavior of almost any Maple library or even built-in procedure. This technique is not often discussed on MaplePrimes, so I suggest that most readers save this Answer in their Favourites.

restart;

`diff/Int/orig`:= eval(`diff/Int`):
unprotect(`diff/Int`);
`diff/Int`:= overload([
    proc(f::algebraic, x::{name, name=range}, M::(identical(method)=name))
        option overload;
        subsindets(
            diff(Int(f,x), args[-1]),
            specfunc(anything, 'Int'),
            subs([_DUMMY= args[3..-2]], J-> 'Int'(op(J), _DUMMY))
        )
    end proc,
    
    `diff/Int/orig`
]):

Test

   
 

 

S:= Int(sin(a*x)/x, x= 0..infinity, method= _DEFAULT);

Int(sin(a*x)/x, x = 0 .. infinity, method = _DEFAULT)

diff(S,a);

Int(cos(a*x), x = 0 .. infinity, method = _DEFAULT)

diff(Int(sin(a*x)/x, x), a);

Int(cos(a*x), x)

I have only overloaded the case where the first argument to Int is type algebraic, the second argument is a name or name= range, and the third argument is method= name. All other cases should behave as they did originally. If you need the overload to work for more general cases of Int, please let me know.

 

 

Download diffInt_overload.mw

You wrote:

Notice that the order of the output from the EulerLagrange command when three functions are input does not correspond to the order of the functions listed, and when they are requested separately, the last function is not included:

The two phenomena are unrelated. As to the order: The EulerLagrange command always gives its result as a set, in curly braces. There is no way to control the order of elements in a Maple set: An internal lexicographic ordering is used so that set-membership checks can be done with binary search and set order is preserved across sessions. (Older Maple used an address-based ordering and set order varied by session.)

You wrote:

T3:=(1/2)*m*(diff(r(t),t)^2+r(t)^2*diff(theta(t),t)^2+diff(z(t),t)^2):
U3:=m*g*z(t):
L3:=T3-U3;

As for the second issue: Your equations above use diff(theta(t), t) but they never use theta(t) undifferentiated. This confuses EulerLagrange, and perhaps it should be rejected as invalid input. If you change diff(theta(t),t) to theta_prime(t) and do not include this in the solution variables, then the union of the individual results for z(t) and r(t) will be the same as the result when both of those are requested.

It can be done like this (I just picked a simple curve for the implicitplot).

P:= plots:-implicitplot(y= exp(x), x= -10..10, y= 0..10^4):
data:= op([1,1], P):
plots:-logplot(convert(data, listlist));
 

Change your final display command to

display([P, f(C)], viewpoint= ["circleright", frames= 32]);

See ?plot3d,viewpoint. The path of the "camera" can be arbitrary; I just picked an easy option above.

Add the option style= patchnogrid or style= wireframe to your implicitplot3d command.

It can be done, for example, like this:

# Returns an Array with n dimensions each with extent 1..m
DynArray:= (n,m)-> Array((1..m)$n, args[3..-1]):
A:= DynArray(3,4);
              

Be careful with this! A has m^n elements. 2^33 = 8 Gig of elements.

First 366 367 368 369 370 371 372 Last Page 368 of 384