When the limit approaches from left, the result must be zero. What cause the wrong? Will introducing an extra parameter M affect the result in W1?

Download limit.mwlimit.mw

eqn := B(n) = -sum(binomial(n + 1, k)*B(k), k = 0 .. n - 1)/(n + 1);
init := B(0) = 1, B(1) = -1/2;
sol := rsolve({eqn, init}, B(n));
Why doesn’t this give me any solution ? Thank you.

How to show that the angle QF2P remains constant when M moves on the ellipse ? Perhaps with Explore ?
restart;
with(plots);
with(geometry);
_EnvHorizontalName := 'x';
_EnvVerticalName := 'y';
x0 := 100;
y0 := 40;
a := 7;
b := 5;
c := sqrt(a^2 - b^2);
ellipse(el, x^2/a^2 + y^2/b^2 - 1);
point(F1, -c, 0);
point(F2, c, 0);
eq := simplify((a^2 - x0^2)*(y - y0)^2 + (b^2 - y0^2)*(x - x0)^2 + 2*x0*y0*(x - x0)*(y - y0)) = 0;
sol := solve({eq}, {y});
line(tang1, op(sol[1]));
line(tang2, op(sol[2]));
sol2 := op(solve({op(sol[1]), x^2/a^2 + y^2/b^2 - 1 = 0}, {x, y}));
xM2 := rhs(sol2[1]);
yM2 := rhs(sol2[2]);
point(A, xM2, yM2);
sol3 := op(solve({op(sol[2]), x^2/a^2 + y^2/b^2 - 1 = 0}, {x, y}));
xM3 := rhs(sol3[1]);
yM3 := rhs(sol3[2]);
point(B, xM3, yM3);
line(Pol, [A, B]);
simplify(Equation(Pol));
isolate(%, y);
xM := 4;
sol := solve({subs(x = xM, x^2/a^2 + y^2/b^2 - 1 = 0)}, {y});
yM := rhs(op(sol[1]));
point(M, xM, yM);
line(Tang, x*xM/a^2 + y*yM/b^2 - 1 = 0);
intersection(P, tang1, Tang);
intersection(Q, tang2, Tang);
line(PF2, [P, F2]);
line(QF2, [Q, F2]);
alpha := FindAngle(PF2, QF2);
arctan(alpha);
evalf(%);
display(textplot([[-c, 0, "F1"], [c, 0, "F2"], [coordinates(B)[], "B"], [coordinates(A)[], "A "], [coordinates(M)[], "M "], [coordinates(P)[], "P "], [coordinates(Q)[], "Q "]], align = {"above", 'right'}), draw([el(color = red), A(color = black, symbol = solidcircle, symbolsize = 16), PF2(color = brown), QF2(color = brown), B(color = black, symbol = solidcircle, symbolsize = 16), M(color = black, symbol = solidcircle, symbolsize = 16), P(color = black, symbol = solidcircle, symbolsize = 16), tang1(color = green), tang2(color = green), Tang(color = green), F1(color = blue, symbol = solidcircle, symbolsize = 16), Q(color = blue, symbol = solidcircle, symbolsize = 16), F2(color = red, symbol = solidcircle, symbolsize = 16)]), axes = none, view = [-7 .. 15, -7 .. 12]);

It seems like there exists a bug when taking the following limit in Maple (I tried Maple 2021):

If I run this command:
> evalf(limit(CylinderU(0,CylinderU(0,x)),x=0));
1.2722774800
the result is 1.2722774800, which seems to be incorrect.
;

However, when I run this command:
> evalf(CylinderU(0,limit(CylinderU(0,x),x=0)));
0.5456799403
the result is 0.5456799403, which seems to be correct.

Finally, when I run this command:
> evalf(CylinderU(0,CylinderU(0,0)));
0.5456799403
the result is 0.5456799403 which is also correct.

My expectation is that all three commands must return the same result, thus I consider this a bug.
I also run the following command in WolframAlpha
> limit ParabolicCylinderU(0,ParabolicCylinderU(0,x)) as x->0.0
and obtained the correct result 0.54568, confirming that in Maple this is evaluated incorrectly.

Would appreciate if anybody can confirm that this is a bug.

How such Maple bug should be reported?

restart;
with(combinat);
F := unapply(rsolve({F(1) = 1, F(2) = 1, F(n + 1) = F(n) + F(n - 1)}, F(n)), n);
combine(expand(F(n + 1)*F(n + 2) - F(n)*F(n + 3)));
F := n -> fibonacci(n);
combine(expand(F(n + 1)^2 - F(n)*F(n + 3) + (-1)^n));
is(F(n + 1)^2 = F(n)*F(n + 2) + (-1)^n);#should be true
G := n -> arctan(1/F(n));
is(G(4) = G(5) + G(6));
is(G(4) = G(5) + G(6));
is(G(2*n) = G(2*n + 1) + G(2*n + 2));#should be true
How to establish these formulas ? Thank you.

How to correct this code ? Thank you.

restart;
with(plots);
with(geometry);
_EnvHorizontalName := 'x';
_EnvVerticalName := 'y';
x0 := 10;
y0 := 9;
a := 7;
b := 5;
c := sqrt(a^2 - b^2);
ellipse(el, x^2/a^2 + y^2/b^2 - 1);
point(F1, -c, 0);
point(F2, c, 0);
point(P, x0, y0);
eq := simplify((a^2 - x0^2)*(y - y0)^2 + (b^2 - y0^2)*(x - x0)^2 + 2*x0*y0*(x - x0)*(y - y0)) = 0;
sol := solve({eq}, {y});
line(tang1, op(sol[1]));
slope(tang1);
line(tang2, op(sol[2]));
slope(tang2);
sol2 := op(solve({op(sol[1]), x^2/a^2 + y^2/b^2 - 1 = 0}, {x, y}));
xM2 = rhs(sol2[1]);
yM2 = rhs(sol2[2]);
point(M2, xM2, yM2);
sol3 := op(solve({op(sol[2]), x^2/a^2 + y^2/b^2 - 1 = 0}, {x, y}));
xM3 = rhs(sol3[1]);
yM3 = rhs(sol3[2]);
point(M3, xM3, yM3);
assume(xM2 - xM3 <> 0);
line(Pol, [M2, M3]);
coordinates(M2);
display(textplot([[-c, 0, "F1"], [c, 0, "F2"], [coordinates(P)[], "P"]], align = {"above", 'right'}), draw([el(color = red), P(color = black, symbol = solidcircle, symbolsize = 16), tang1(color = green), tang2(color = green), Pol(color = red, thickness = 3), F1(color = blue, symbol = solidcircle, symbolsize = 16), F2(color = red, symbol = solidcircle, symbolsize = 16)]), axes = none);
Warning, data could not be converted to float Matrix

Why this code does not work. Thank you.
restart;
with(plots);
with(geometry);
_EnvHorizontalName := 'x';
_EnvVerticalName := 'y';
x0 := 10;
y0 := 9;
a := 7;
b := 5;
c := sqrt(a^2 - b^2);
ellipse(el, x^2/a^2 + y^2/b^2 - 1);
point(F1, -c, 0);
point(F2, c, 0);
point(P, x0, y0);
eq := simplify((a^2 - x0^2)*(y - y0)^2 + (b^2 - y0^2)*(x - x0)^2 + 2*x0*y0*(x - x0)*(y - y0)) = 0;
eq := (a^2 - x0^2)*m^2 + 2*x0*y0*m + b^2 - y0^2 = 0;
sol := solve(%, m);
m1 := sol[1];
m2 := sol[2];
(y - y0 - m1) + (-x + x0);
(y - y0 - m2) + (-x + x0);
line(tang1, (y - y0 - m1) + (-x + x0));
line(tang2, (y - y0 - m2) + (-x + x0));
display*[textplot*([[-c, 0, "F1"], [c, 0, "F2"], [coordinates(P)[], "P"]], align = {"above", 'right'}), draw*([el(color = red), P(color = black, symbol = solidcircle, symbolsize = 16), tang1(color = green), tang2(color = green), F1(color = blue, symbol = solidcircle, symbolsize = 16), F2(color = red, symbol = solidcircle, symbolsize = 16)], axes = none)];
 

sol := y = -3283/4253 - (3283*x)/4253, How can I determine the value of the coefficient of x?
How can I take the value of the coefficient of x? Thank you.

restart;
with(geometry);
with(plots);
_EnvHorizomtalName = 'x';
_EnvVerticalName = 'y';
_local(D);
line(delta, y = 1/3*x - 2, [x, y]);
line(deltap, y = (-1)/4*x + 1, [x, y]);
line(D, y = 3*x - 5, [x, y]);
point(S, 3, 0);
omega := Pi/3;
intersection(P, delta, D);
intersection(Pp, deltap, D);
projection(H, S, delta);
projection(K, S, deltap);
projection(M, S, D);
circle(c1, [H, K, M], 'centername' = O1);
display*[textplot*([[coordinates(S)[], "S"], [coordinates(P)[], "P"], [coordinates(H)[], "H"], [coordinates(K)[], "K"], [coordinates(M)[], "M"]], font = [times, bold, 16], align = [above, right]), draw*([delta(color = blue), deltap(color = blue), D(color = red), c1(color = black), S(color = black, symbol = solidcircle, symbolsize = 16), P(color = black, symbol = solidcircle, symbolsize = 16)], scaling = constrained, axes = none, view = [-15 .. 15, -15 .. 15])];
               [                /[             [9  -13     ]  
plots:-display [plots:-textplot |[[3, 0, "S"], [-, ---, "P"], 
               [                \[             [8   8      ]  

  [33  -9     ]  [52  4      ]  [9  2     ]]  
  [--, --, "H"], [--, --, "K"], [-, -, "M"]], 
  [10  10     ]  [17  17     ]  [5  5     ]]  

                                                  \              
  font = [times, bold, 16], align = [above, right]|, geometry:-d\
                                                  /              

  raw ([delta(color = blue), deltap(color = blue), 

  D(color = red), c1(color = black), 

  S(color = black, symbol = solidcircle, symbolsize = 16), 

  P(color = black, symbol = solidcircle, symbolsize = 16)], 

  scaling = constrained, axes = none, view = [-15 .. 15, -15 .. 15])];
  display*[textplot*([[3, 0, "S"], [9/8, -13/8, "P"], [33/10, -9/10, "H"], [52/17, 4/17, "K"], [9/5, 2/5, "M"]], font = [times, bold, 16], align = [above, right]), draw*([delta(color = blue), deltap(color = blue), D(color = red), c1(color = black), S(color = black, symbol = solidcircle, symbolsize = 16), P(color = black, symbol = solidcircle, symbolsize = 16)], scaling = constrained, axes = none, view = [-15 .. 15, -15 .. 15])]
Why the program will not run; Thank tou.

restart;
with(geometry);
with(plots);
_EnvHorizontalName := 'x';
_EnvVerticalName := 'y';
cb := color = blue;
t3 := thickness = 3;
l3 := linestyle = 3;
xA := 3;
yA := 0;
xB := -3;
yB := 0;
c := 6;
point(A, xA, yA);
point(B, xB, yB);
R := 5;
alpha := arctan(3/4);
evalf(%*180/Pi);
                          36.86989765

xH := (xA + xB)/2;
yH := (yA + yB)/2;
point(H, xH, yH);
xO1 := 0;
yO1 := 4;
point(O1, xO1, yO1);
xO2 := 0;
yO2 := -4;
point(O2, xO2, yO2);
segment(sAB, A, B);
segment(sHO, H, O1);
segment(sAO, A, O1);
segment(sBO, B, O1);
                              sAB

alpha1 := arctan((yO1 - yA)/(xO1 - xA));
beta := Pi + arctan((yO1 - yB)/(xO1 - xB));
AR1 := plottools[arc]([xO1, yO1], R, alpha1 .. beta, l3);
AR2 := plottools[arc]([xO2, yO2], R, -beta .. -alpha1, l3);
N := 80;
dt := (beta - alpha1)/(N - 1);
dr := draw({O1, O2, sAB, sHO, sAO(cb), sBO(cb)});
tex := textplot([[xA, yA - 0.5, "A"], [xB, yB - 0.5, "B"], [xO1, yO1 + 0.5, "O1"], [xO2, yO2 - 0.5, "O2"], [xH, yH - 0.5, "H"]]);
M1 := seq(plottools[disk]([R*cos(dt*t + alpha1) + xO1, R*sin(dt*t + alpha1) + yO1], 0.2, color = orange), t = 0 .. N);
M2 := seq(plottools[disk]([R*cos(dt*t + Pi + alpha1) + xO2, R*sin(dt*t + Pi + alpha1) + yO2], 0.2, color = orange), t = 0 .. N);
P1 := seq(plottools[polygon]([[R*cos(dt*t + alpha1) + xO1, R*sin(dt*t + alpha1) + yO1], [xA, yA], [xB, yB]], color = aquamarine, linestyle = dash), t = 0 .. N);
P2 := seq(plottools[polygon]([[R*cos(dt*t + Pi + alpha1) + xO1, R*sin(dt*t + Pi + alpha1) + yO1], [xA, yA], [xB, yB]], color = aquamarine, linestyle = dash), t = 0 .. N);
for t to N do
    E[t] := display(dr, tex, AR1, AR2, M1[t], P1[t]);
    F[t] := display(dr, tex, AR1, AR2, M2[t], P2[t]);
end do;
display([seq(E[t], t = 1 .. N), seq(F[t], t = 1 .. N)], insequence = true, axes = none, scaling = constrained, view = [-10 .. 10, -10 .. 10]);

NULL;
angle that does not turn on the lower arc. Thank you.

trouverTripletsDecroissants := proc(tripletInitial) local m, n, a, b, c, triplets, dernierTriplet; triplets := []; dernierTriplet := tripletInitial; m := dernierTriplet[3]; do m := m - 1; for n to m - 1 do if igcd(m, n) = 1 and (m - n) mod 2 = 1 then a := m^2 - n^2; b := 2*m*n; c := m^2 + n^2; if a^2 + b^2 = c^2 and a < dernierTriplet[1] and b < dernierTriplet[2] and c < dernierTriplet[3] then dernierTriplet := [a, b, c]; triplets := [op(triplets), dernierTriplet]; break; end if; end if; end do; break; if 0 < nops(triplets); end do; return triplets; end proc;
tripletInitial := [275, 252, 373];
tripletsDecroissants := trouverTripletsDecroissants(tripletInitial);
print(tripletsDecroissants);
               tripletInitial := [275, 252, 373]

           tripletsDecroissants := [[273, 136, 305]]

                       [[273, 136, 305]]

;
trouverTripletsDecroissants(275, 252, 373);
Error, (in trouverTripletsDecroissants) final value in for loop must be numeric or character
How to correct this error. Thank you.

p eriode:=proc(r::rational)""  local a,b,c,f,i,k,l,p,q,s:  s:=couvert(evalf(abs(r)-trunc(r),50),string):  if  s[1]="." then s:=s[2..length(s)] fi:  a:=numer(simplify(r)): b:=denom(simplify(r)):  if  igcd(b,10)=1 then       q:=0:       p:=1:      while (10^(p) mod b)<>1 do  p:=p+1 od:  else      f:=ifactors(b)[2]:      k:=0:l:=0:      for i  ;
to nops(f) do
        if f[i][1]=2 then k:=f[i][2] fi:
        if f[i][1]=5 then l:=f[i][2] fi: 
     od:
   c:=b/(2^k*5^l):
   q:=max(k,l):
   if c=1 then p:=0 else p:=op(2,periode(1/c)) fi:
fi:
[q,p,[seq(parse(s[k]),k=1+q..q+p)]]
end:
Error, improper op or subscript selector
Error, reserved word `fi` unexpected
I don't understand these errors. Thank you.

I'm working on simulating a triple pendulum. I have the numeric solution of the ODEs and some nice plots that indicate I'm on the right track. I would to animate these solutions. 
It seems like there is a way to plot the positions over time, as well as the lines between the points. If anyone can put me on the right track towards this, I would really appreciate it. 

I have a set of 3 differential equations called "odesys" that I need to solve. I input this line to solve them and I get the following error. 

sols := dsolve([odesys, ICs], {Theta[1](t), Theta[2](t), Theta[3](t)})
Error, (in dsolve) not a system with respect to the unknowns {Theta[1](t), Theta[2](t), Theta[3](t)}
 

I've also tried it using the exact syntax that the thetas are in the equation and I get this error. 

sols := dsolve([odesys, ICs], {Theta(t)[1], Theta(t)[2], Theta(t)[3]});
Error, (in dsolve) required an indication of the solving variables for the given system
 

I don't know even know which error message is better. The equations are second order and non-linear if that has any bearing on anything. 

How to show that any Pythagorean triplet can be obtained from <3,4,5> ? Thank you.
Can we simplify this program?

#Génération Géométrique et Algébrique des triplets Pythagoriciens
restart;
with(geometry);
with(LinearAlgebra);
_EnvHorizomtalName = 'x';
_EnvVerticalName = 'y';

with(plottools);
P := point([0, 0], color = black, symbol = cross, symbolsize = 25);
Oo := point([1/2, 1/2], color = black, symbol = cross, symbolsize = 25);
A := point([1, 1/2], color = black, symbol = cross, symbolsize = 25);
with(plots);
c1 := circle([1/2, 1/2], 1/2, color = blue);
NULL;
PA := line([0, 0], [1, 1/2], color = red);
eqC := (x - 1/2)^2 + (y - 1/2)^2 = 1/4;
eqPA := y = 1/2*x;
sol := solve({eqC, eqPA}, {x, y});

t1 := textplot([0, 0, 'typeset'("P"), font = [Times, Bold, 14]], 'align' = 'above');
t2 := textplot([1, 1/2, 'typeset'("A"), font = [Times, Bold, 14]], 'align' = 'right');
t3 := textplot([1/5, 1/10, 'typeset'("A'"), font = [Times, Bold, 14]], 'align' = 'above');
A1 := point([1 - 1/5, 1/10], color = black, symbol = cross, symbolsize = 25);
diff(A, x) := point([1/5, 1/10], color = black, symbol = cross, symbolsize = 25);
t4 := textplot([1 - 1/5, 1/10, 'typeset'("A1"), font = [Times, Bold, 14]], 'align' = 'right');
A2 := point([1 - 1/5, 1 - 1/10], color = black, symbol = cross, symbolsize = 25);
t5 := textplot([1 - 1/5, 1 - 1/10, 'typeset'("A2"), font = [Times, Bold, 14]], 'align' = 'right');
A3 := point([1/5, 1 - 1/10], color = black, symbol = cross, symbolsize = 25);
t6 := textplot([1/5, 1 - 1/10, 'typeset'("A3"), font = [Times, Bold, 14]], 'align' = 'right');

poly := Matrix([[1/5, 1/10], [1 - 1/5, 1/10], [1 - 1/5, 1 - 1/10], [1/5, 1 - 1/10]], datatype = float);
pol := polygonplot(poly, color = blue, transparency = 0.95);

display(c1, P, Oo, A, PA, seq(A || i, i = 1 .. 3), seq(t || i, i = 1 .. 6), pol, scaling = constrained, axes = none, size = [600, 600]);

R1 := Transpose(<<1, -2, 2> | <2, -1, 2> | <2, -2, 3>>);
R2 := Transpose(<<1, 2, 2> | <2, 1, 2> | <2, 2, 3>>);
R3 := Transpose(<<-1, 2, 2> | <-2, 1, 2> | <-2, 2, 3>>);
V := <3, 4, 5>;
MatrixVectorMultiply(R1, V);
MatrixVectorMultiply(R2, V);
MatrixVectorMultiply(R3, V);
t1 := <2225, 3648, 4273>;
t2 := MatrixVectorMultiply(1/R1, t1);
t3 := MatrixVectorMultiply(1/R3, t2);
t4 := MatrixVectorMultiply(1/R2, t3);
t5 := MatrixVectorMultiply(1/R1, t4);
t6 := MatrixVectorMultiply(1/R3, t5);
MatrixVectorMultiply(MatrixMatrixMultiply(MatrixMatrixMultiply(MatrixMatrixMultiply(MatrixMatrixMultiply(R1, R3), R2), R1), R3), V);
% - t1;
NULL;