The simplification is not always valid for complex v. So you need an extra assumption:

simplify(sqrt(r^2*(R+r*cos(v))^2)) assuming 0 < r < R, v::real;
                        r (R + r cos(v))

If you use exact quantities, the examples work:

with(geom3d):

point(A, 7, 29, 7): point(B, 10, 27, 13/2): point(C, 30, 47, 13/2): point(D1, 32, 49, 13/2):
line(f, [D1, C]): line(g, [A, B]):
intersection(P, f, g):

coordinates(P); # the lines intersect at B
                                [10, 27, 13/2]

But floating-point computations fail, which is indeed a bug, as there is an unresolved variable in the output:

point(A, 7, 29, 7): point(B, 10, 27, 13/2.): point(C, 30, 47, 13/2.): point(D1, 32, 49, 13/2.):
line(f, [D1, C]): line(g, [A, B]):
intersection(P, f, g):

coordinates(P);
           [32 - 0.7071067811 t1, 49 - 0.7071067811 t1, 6.500000000]

The fact that this doesn't find an intersection is possibly a bug too:

Digits := 20:
point(A, 7, 29, 7): point(B, 10, 27, 13/2.): point(C, 30, 47, 13/2.): point(D1, 32, 49, 13/2.):
line(f, [D1, C]): line(g, [A, B]):
intersection(P, f, g): # careful, P is now NULL
intersection: the given objects do not intersect

patmatch is structural, that is, the op()s in the pattern should match the op()s in the expression.

op([0 .. -1], 1/y^n::anything);
                          `^`, y^n::anything, -1

op([0 .. -1], 1/y^2);
                                 ^, y, -2

1/y^n::anything is y^n::anything to the power of -1. In 1/y^2 there is no minus first power of anything, so there won't be a match.

y^(-n::anything) won't match 1/y^2 for the same reason:

op([0 .. -1], y^(-n::anything));
                           ^, y, -n::anything

op([0 .. -1], -n::anything);
                           *, -1, n::anything

The exponent is -1 times n::anything. In 1/y^2 there is no -1 times anything.

(An exception to the general rule is that the pattern matcher can try to rewrite x as x+0 or 1*x or x^1.)

This will work:

patmatch(x^2/y^2, conditional(x^m::integer*y^n::integer, m+n = 0));
                                true

(although there is a strange disclaimer in ?patmatch that "it is not possible to use `=` or `<>`").

@AndreaAlp

When you're eliminating a variable from the equations of two surfaces, you're not getting the equation of their intersection. What you can do is rotate the intersection plane so that it becomes horizontal, and then you have a 2D problem:

ccent:=proc(quadeq, lineq, vars)
  local nv, rotm, rotsub, zsub, c;
  nv := coeff~(lineq, vars);
  rotm := Student:-LinearAlgebra:-RotationMatrix(
    -VectorAngle(nv, [0, 0, 1]),
    LinearAlgebra:-CrossProduct(nv, [0, 0, 1]));
  rotsub := vars =~ convert(rotm.Vector(vars), list);
  zsub := solve(subs(rotsub, vars[1..2] =~ 0, lineq), {vars[3]});
  conic(c, subs(rotsub, zsub, quadeq), vars[1..2]);
  rotm.Vector([(op@coordinates@center)(c), subs(zsub, vars[3])])
end proc:

pp := ccent(y-model_M_femoral_condyle, z-model_M_tibial_plateau, [x, y, z]);
            Vector[column](3, [50.3331780498402, 493.551723510568, 82.5916599287289])

@AndreaAlp

In LinearFit, the dependent variable is the last column. So if you're looking for an equation of the form y=f(x,z), you need to transpose the columns. It has nothing to do with how you name the variables.

with(Statistics):

clip := proc(arr, xr, yr, zr)
  local crits, ans;
  crits := map(proc(ir)
      if ir[2, 1] = -infinity then
        m -> min(m[.., ir[1]]) < ir[2, 2]
      elif ir[2, 2] = infinity then
        m -> max(m[.., ir[1]]) > ir[2, 1]
      else m -> ir[2, 1] < Mean(m[.., ir[1]]) < ir[2, 2]
      end if end proc,
    [[1, xr], [2, yr], [3, zr]]);
  ans := select(
    m -> crits[1](m) and crits[2](m) and crits[3](m),
    [seq](arr[i], i = ArrayDims(arr)[1]));
  Array(1 .. nops(ans), 1 .. 3, 1 .. 3, (i, j, k) -> ans[i][j, k])
end proc:

polyplot := arr -> plots:-polygonplot3d([seq](arr[i], i = ArrayDims(arr)[1]), style = surface):

femur := Import("Femur.stl", output = triangles);

femur := Array(ArrayDims(femur), 1 .. 3, 1 .. 3, (i, j, k) -> femur[i][j, k]);

XYZ_COM_femur := [seq]((Mean@@2)(femur[.., .., i]), i = 1 .. 3);

min_condyle := min(femur(.., .., 2));

L_femoral_condyle := clip(femur,
  [XYZ_COM_femur[1]-40, XYZ_COM_femur[1]-10],
  [-infinity, min_condyle+15],
  [-infinity, infinity]);

n_points_L_femoral_condyle :=
  upperbound(L_femoral_condyle, 1)-lowerbound(L_femoral_condyle, 1)+1;

model_L_femoral_condyle := LinearFit(
  [1, x, z, x^2, z^2],
  ArrayTools:-Reshape(L_femoral_condyle, n_points_L_femoral_condyle*3, 3)[
    .., [1, 3, 2]], # xzy transposition happens here
  [x, z]);

plots:-display(
  polyplot(L_femoral_condyle),
  plot3d([x, model_L_femoral_condyle, z],
    x = min(L_femoral_condyle[.., .., 1]) .. max(L_femoral_condyle[.., .., 1]),
    z = min(L_femoral_condyle[..,  .., 3]) .. max(L_femoral_condyle[ .., .., 3]),
   color = red));

Groebner:-NormalForm computes the remainder from polynomial division:

sb := Groebner:-NormalForm(Phi[1], [tobesbstd] -~ [sbeqmaple18], (tdeg@op@indets)(Phi[1])): # in 2017.3
simplify(sb-`sbΦ`[1]);
                                   0

But in general, if the list of the polynomials is not a Groebner basis, the remainder will not be unique: if you rearrange the list, the answer may be different.

Regarding algsubs, the help does mention that it needs the expanded form, but I'm not sure what it's supposed to do exactly:

algsubs(x^2 = x*y, x^3);
                                 y*x^2

Was it giving the same result (as opposed to y^2*x) in Maple 18? If yes, I don't understand what is meant by "this is repeated until the leading monomial in f is less than the leading monomial in a" in the documentation (and how that would work with something like x^2=x^3*y).

You're plotting my_plane=0, but then you're solving my_plane=z. And in my_center_point you're doing [xcenter, zcenter, ...] instead of [xcenter, solve(subs(x=xcenter, z=zcenter, my_plane)), zcenter].

Regarding implicitplot3d, my_quadric is bounded in x and z, so you can estimate the plot range easily after doing CompleteSquare. Or just by doing Optimization:-Maximize(x, {my_quadric=0}), etc. If the plot range is too wide, all sampling points may happen to lie outside the cylinder, and the plot will be empty. numpoints or grid can help, e.g., grid=[100, 5, 100]. 'view' won't help because that just clips unwanted points, and without dense enough sampling you won't have any points at all.

You can construct the continuous function by integrating the derivative:

f := u -> arctan((u-sqrt(u^2+81)*cot(sqrt(u^2+81))*tan(u))/(u*tan(u)+sqrt(u^2+81)*cot(sqrt(u^2+81))));

f2 := unapply(Int(simplify(D(f)(xi)), xi = 0 .. u), u);

plot(f2(u), u = 0 .. 50);

Or you can write the same as a dsolve.

Alternatively, you can compute the positions of the discontinuities:

rts := ([fsolve])(u*tan(u)+sqrt(u^2+81)*cot(sqrt(u^2+81)), u = 0 .. 50, maxsols = infinity);
            [1.519814744, 5.673379054, 24.50237178]

rts := RootOf~(u*tan(u)+sqrt(u^2+81)*cot(sqrt(u^2+81)), u, rts);

f3 := unapply(f(u)-Pi*add(Heaviside~(u-~rts)), u);

plot(f3(u), u = 0 .. 50);

Could be written more naturally as a piecewise, but RootOf has issues with less/greater comparisons.

Just Distribution(ProbabilityFunction= fY) works, since fY is already a function.

You have a few typos in your code: fX(X) should be fY(X); opt=(x->x) needs parentheses, otherwise it's interpreted as (opt=x)->x; x->unapply(y,x) becomes x->x->y when you need just x->y:

simp := e -> (simplify(expand(e)) assuming n::integer);
BB := Distribution(ProbabilityFunction = unapply(convert(fY(k), GAMMA), k));

simp(Mean(BB));
       n*a/(a+b)

sum((k-simp(Mean(BB)))^2*ProbabilityFunction(BB, k), k = 0 .. infinity);
       n^2*a^2*GAMMA(n+b)*GAMMA(a+b)*hypergeom([a, -n, (-a*n+a+b)/(a+b), (-a*n+a+b)/(a+b)],
       [-n+1-b, -n*a/(a+b), -n*a/(a+b)], 1)/((a+b)^2*GAMMA(b)*GAMMA(a+b+n))

This shows two possible methods. For small positive c, there are two positive solutions, asymptotic to 1+5/9*c and to (3/5*1/c)^(1/4)-1/12*sqrt(15*c). For small negative c, there is only the first solution. The leading term is immediately obvious from the dominant balance.

The value of the overshoot is exactly the same for other functions. That is, if f and f' are continuous except for a finite number of finite jumps, then if the value of a jump of f is A, the difference between the maximum and the minimum of the nth Fourier sum near the jump will tend to 2/Pi*Si(Pi)*A when n goes to infinity. I'm not sure what you mean about the rectangle and 1-x being different.

The infinite Fourier series either converges to the midpoint or diverges at the jump. In the latter case, it's still summable to the midpoint. So there is no Gibbs phenomenon for the infinite series.

If you want to compare the convergence, consider:

simplify(1/2*int((1-abs(x))*exp(-I*Pi*k*x), x = -1 .. 1));
                             (1-cos(k*Pi))/(Pi^2*k^2)

simplify(1/2*int((1-x)*exp(-I*Pi*k*x), x = 0 .. 1));
                   (1/2)*(-I*k*Pi-exp(-I*k*Pi)+1)/(Pi^2*k^2)

The first one is O(1/k^2), the second is only O(1/k). No surprise that the first one is going to converge faster. You can also compare the L2 norms as a way to measure the differences. Plus, for discontinuous functions there is such thing as the Gibbs phenomenon.

I suppose the explanation is that the bug hasn't been fixed yet. The Fourier integral exists in the sense of ordinary functions, you can verify numerically that Pi*GAMMA(0, abs(w)) is correct. Most likely happens because of this:

int(arctan(x)*cos(w*x)/x, x = 0 .. infinity, method = MeijerG);
                                   (1/4)*Pi*exp(-sqrt(w^2))

Probably isn't too difficult to fix, as this is a standard MeijerG convolution case, giving sqrt(Pi)/4*MeijerG([[], [1]], [[1/2, 0, 0], []], w^2/4). Also:

int(arctan(x)*exp(-I*w*x)/x, x = 0 .. infinity); # incorrect for w<0
         (1/8*I)*w*MeijerG([[0], [1/2]], [[0, 0, -1/2, -1/2], []], -(1/4)*w^2)/sqrt(Pi)

This one is more tricky.

The inverse transform can be computed like this:

int(GAMMA(0, w)*cos(x*w), w = 0 .. infinity);
                       -(1/2)*(-Pi*sqrt(x^2)+2*x*arctan(1/x))/x^2

This is explained in ?algcurves,Weierstrassform. By default the output is the curve y^2+(polynomial in x)=0. It doesn't say anything about the sign of the leading term in x. If you want Weierstrass normal form, you need to set the option:

algcurves:-Weierstrassform(x^3+y^2-7*x+88, x, y, x0, y0);
                                          [x0^3+y0^2-7*x0+88, x, y, x0, y0]

algcurves:-Weierstrassform(x^3+y^2-7*x+88, x, y, x0, y0, Weierstrass);
                                          [-4*x0^3+y0^2+28*x0+352, -x, 2*y, -x0, (1/2)*y0]

Indeed, there are zeros of the denominator with Re(s)>0:

A11f := factor(A11):

a0 := max(Re~([fsolve](denom(A11f), complex)));
               0.00353314997991911033113733406816

So that'll be the dominant exponential term:

A11pf := convert(A11f, fullparfrac, s):

B11 := inttrans[invlaplace](A11pf, s, t):

plot(B11*exp(-a0*t), t = 0 .. 20000); # bounded oscillations

The root can be verified by setting Digits higher or by working with convert(A11, rational, exact), so factor() can be trusted here.