## 12705 Reputation

8 years, 270 days

## plot3d...

p1:=plot3d([[x,y,-y^2],[x,y,x^2]], x=0..1,y=0..x):
p2:=plot3d([[1,y,z]], z=-y^2..1,y=0..1,transparency=0.5):
p3:=plot3d([[x,0,z]], x=0..1,z=0..x^2,transparency=0.5):
p4:=plot3d([[x,x,z]], x=0..1,z=-x^2..x^2,transparency=0.5):
plots[display](p1,p2,p3,p4,labels=["x","y","z"]);

# less memory used

## curry...

myPDF := (v,t) -> piecewise(t < 0, 0, t < v, 1/v, 0):
myD:= Distribution(PDF = curry(myPDF, 1)):  #v=1
Mean(myD);
1/2

## Change...

J:=Int( exp(-I*k*x)/cosh(x), x=-infinity..infinity):
value(IntegrationTools:-Change(J,exp(x)=t));

## To be fixed...

1. Due to the wonderful 2D math, a variable named _ appeared. This must be corrected.

2. The equation for solve is too complicated to be solved symbolically.

3. For fsolve you must have numerical values for parameters.

## conversion...

PW1:=proc(e::specfunc(anything,piecewise))
local f, o:=op(e);
f:=proc()
if   nargs=3 then '`if`'(args)
elif nargs=2 then '`if`'(args,0)
else '`if`'(args[1],args[2],f(args[3..-1])) fi
end;
f(o)
end:

PW:=e->subsindets(e, specfunc(anything, piecewise), PW1):
####################

PW(piecewise(x<1, 10, x<2, 20, x<3, 30));

`if`(x < 1, 10, `if`(x < 2, 20, `if`(x < 3, 30, 0)))

f:=piecewise(x < 0, piecewise(y < 0, 0, 1), piecewise(y < 0, 1, 0)):
PW(f);
`if`(x < 0, `if`(y < 0, 0, 1), `if`(y < 0, 1, 0))

## Without packages...

L:=(u,v)->u(v(f_(x,t)))-v(u(f_(x,t))):
alias(``=f_(x,t));
#########################
v:=f->diff(f,x):
w:=f->x*diff(f,t):
L(v,w);

u:=f->x*t*diff(f,x,t):
L(u,w):  simplify(%);

## Compiled version...

I was curious to see how much the compiler can help in such problems.

It seems that in Carl's version (which is very fast) the compiler is not very useful.
I used practically my original version, made compilable.

Here it is:

Tri1:=proc(n::integer,V::Vector(datatype= float[8]),
X::Vector(datatype= float[8]),Y::Vector(datatype= float[8]),
Ax::float[8],Ay::float[8],Bx::float[8],By::float[8],Cx::float[8],Cy::float[8])
local k::integer,a::float[8],b::float[8];
for k to n do
a:=V[k]; b:=V[n+k];
if a+b>1 then a:=1-a;b:=1-b fi;
X[k]:=Ax+(Bx-Ax)*a+(Cx-Ax)*b;
Y[k]:=Ay+(By-Ay)*a+(Cy-Ay)*b
od;
end:

CTri1:=Compiler:-Compile(Tri1):

Tri:=proc(n,Ax,Ay,Bx,By,Cx,Cy)

local
V:=LinearAlgebra:-RandomVector(2*n, generator= 0..1., datatype= float[8]),
X:=Vector(n,datatype= float[8]), Y:=Vector(n,datatype= float[8]);
global Ctri1;
CTri1(n,V,X,Y,Ax,Ay,Bx,By,Cx,Cy);
X,Y
end:
####
Ax:=-1: Ay:=0:  # Triangle T = ABC
Bx:=0:  By:=4:
Cx:=2:  Cy:=1:
n:=6000:           # number of random points
T:=CodeTools:-Usage(Tri(n,Ax,Ay,Bx,By,Cx,Cy),iterations=100):
memory used=195.55KiB, alloc change=18.38MiB, cpu time=940.00us, real time=920.00us, gc time=0ns

Compare with Carl's version (not compiled)

T:= CodeTools:-Usage(SampleTriangle([[-1,0],[0,4],[2,1]], 6000), iterations=100):
memory used=0.86MiB, alloc change=8.21MiB, cpu time=24.18ms, real time=6.52ms, gc time=1.25ms

Note that I have used iterations=100, otherwise the timing was not stable enough.

## Definition...

The definition is simple and clearly presented in the help page.

HeunT is the solution of the differential equation (depending on three parameters):

f(0)=1,  f'(0)=0

It is an entire function. For other properties:

## sum...

Yes, it seems that the numerical algorithm in sum fails here.

But we can verify manually (after a simple estimation of the remainder):

evalf(evalf[600](%-29/2));
-4.142437940*10^(-531)

## triangle and ellipse...

Ax:=-1: Ay:=0:  # Triangle T = ABC
Bx:=0:  By:=4:
Cx:=2:  Cy:=1:
n:=4000:        # number of random points

a:=Vector(n): b:=Vector(n):
r:=rand(0.0 .. 1.0):
for i to n do
x:=r(): y:=r():
if x+y>1 then  x,y := 1-x,1-y  fi;
a[i]:=x: b[i]:=y:
od:
u:=Vector(n,1):

plot(
Ax*u+(Bx-Ax)*a+(Cx-Ax)*b, Ay*u+(By-Ay)*a+(Cy-Ay)*b,
symbol= point,style=point,scaling=constrained);

# The distribution is uniform i.e.
# Prob(X in D) = Area(D)/Area(T).

n:=8000:   #ellipse - nonuniform distribution (naive attempt)
a:=Vector(n,i->r()):
b:=Vector(n,i->r()):
plot(3*a*~cos~(2*Pi*b),2*a*~sin~(2*Pi*b),
symbol= point,style=point,scaling=constrained);

a:=sqrt~(a):   #ellipse - uniform distribution
plot(3*a*~cos~(2*Pi*b),2*a*~sin~(2*Pi*b),
symbol= point,style=point,scaling=constrained);