Rouben Rostamian

Sliding with friction...

Answered: Rouben Rostamian 8813

21 hours ago

4 4

Earl, here is how I would solve the problem.

We have a curve (a track) in the vertical plane, and a bead that slides along the track. The coefficient of friction is μ. We wish to describe the bead's motion.

It's possible to derive the equation of motion entirely in Maple as you are attempting to do, however the reasoning is more transparent when done by the proverbial "paper and pencil" as in the attached sliding-with-friction.pdf. The calculation is generic, in the sense that it applies to any track shape and any choice of parameters.

In the attached Maple worksheet I have implemented the resulting equation, and then illustrated by applying it to the specific case of a circular track. The bead slides down the circle, goes past the bottom, and oscillates a few times before coming to a standstill abruptly in finite time due to friction. You may change the shape of the track and other parameters as you wish. They are all defined in the statement params := ... in the worksheet.

>

restart;

>	with(plots):

>	Typesetting:-Settings(typesetdot=true):

>	r := s -> < x(s), y(s) >;

(1)

The tangent vector

>	T := diff(r(s),s);

(2)

The length of the tangent vector

>	T_len := sqrt(T^+ . T);

((diff(x(s), s))^2+(diff(y(s), s))^2)^(1/2)

(3)

The unit tangent

>	T / T_len: convert(%, D); tau := unapply(%, s):

$Vector(2, {(1) = (D(x))(s)/sqrt((D(x))(s)^2+(D(y))(s)^2), (2) = (D(y))(s)/sqrt((D(x))(s)^2+(D(y))(s)^2)})$

(4)

The unit normal

>	<-tau(s)[2], tau(s)[1]>: convert(%, D); nu:= unapply(%, s):

$Vector(2, {(1) = -(D(y))(s)/sqrt((D(x))(s)^2+(D(y))(s)^2), (2) = (D(x))(s)/sqrt((D(x))(s)^2+(D(y))(s)^2)})$

(5)

The expression :

>	diff(r(s),s): convert(%, D): rp := subs(s=s(t ), %);

(6)

The expression :

>	diff(r(s),s,s): convert(%, D): rpp := subs(s=s(t), %);

(7)

The length of rp

>	rp_len := sqrt(rp^+ . rp);

(8)

The velocity

>	v := diff(r(s(t)),t);

(9)

The gravitational pull:

>	gvec := <0, -g>;

(10)

Here is Q (see documentation)

>	simplify(tau(s(t))^+ . v); Q := mu*signum(%);

(11)

Newton's Law of Motion: Mass times acceleration = applied forces

See the documentation for the derivation.

>	rp_lendiff(s(t),t,t) + rpp^+ . ( tau(s(t)) + Qnu(s(t)) ) * diff(s(t),t)^2 - gvec^+ . ( tau(s(t)) + Q*nu(s(t)) ): de := simplify(%);

(mu*(((D(x))(s(t))*((D@@2)(y))(s(t))-(D(y))(s(t))*((D@@2)(x))(s(t)))*(diff(s(t), t))^2+g*(D(x))(s(t)))*signum((diff(s(t), t))*((D(x))(s(t))^2+(D(y))(s(t))^2)^(1/2))+((D(x))(s(t))^2+(D(y))(s(t))^2)*(diff(diff(s(t), t), t))+((D(x))(s(t))*((D@@2)(x))(s(t))+(D(y))(s(t))*((D@@2)(y))(s(t)))*(diff(s(t), t))^2+(D(y))(s(t))*g)/((D(x))(s(t))^2+(D(y))(s(t))^2)^(1/2)

(12)

>

Example

Everything up to this point is generic, as we have not made any assumptions

on the shape of the curve or other parameters. Here we work out a specific

example.

>	params := g=1, mu=1/10, x = (s-> -cos(s) + 1), y = (s-> -sin(s) + 1) ;

(1.1)

Here is what the track and the starting point look like:

>	subs(params, [x(s), y(s), s=0..2*Pi]): track := display( plot(%, color="Blue", thickness=3), pointplot(eval([x(0),y(0)], {params}), symbol=solidbox, symbolsize=30, color="Green"), scaling=constrained);

and here is the differential equation for this particular instance:

>	subs(params, de): DE := simplify(%);

(1.2)

Solve the DE:

>	dsol := dsolve({DE, s(0)=0, D(s)(0)=0}, numeric, output=operator);

(1.3)

Plot s(t):

>	odeplot(dsol, [t,s(t)], t=0..20);

Warning, cannot evaluate the solution further right of 16.400954, maxfun limit exceeded (see ?dsolve,maxfun for details)

The solution extends up to the time where the bob stops moving.

We set tmax according to Maple's warning message above;

>	tmax := 16.400954;

(1.4)

This proc displays the bob at a time t:

>	frame := proc(t) local S; S := subs(dsol, s); subs(params, [x(S(t)), y(S(t))]); pointplot(%, symbol=solidcircle, symbolsize=30, color="Red"); end proc:

Animate the bob's motion. Friction stops the bob after a few oscillations.

>	animate(frame, [t], t=0..tmax, background=track, scaling=constrained, frames=60);

>

Download sliding-with-friction.mw

Unfortunately not possible...

Answered: Rouben Rostamian 8813

August 05 2026

1 1

Maple's export to PDF has been broken for a long time, and despite my occasional prodding, it remains useless.

Try, for instance:

restart;
with(plottools):
plots:-display(disk([0, 0], 1), disk([1, 0], 1), color = [red, green], transparency = 0.8, scaling = constrained);

Right-lick on the image and select Export -> PDF, and save.

Viewing the saved result we see two glaring shortcomings:

The transparency is lost. That's despite the fact that the image is encoded as a PDF 1.4 which has native support for transparency.
The graphics is embedded in a 612 x 792 pts (US letter size) "paper". For the life of me, I have no idea how graphics embedded in a letter-size page can be of use to anyone. Shouldn't the graphics be saved in dimensions compatible with its bounding box? I seem to recall that Maple 11 used to do that properly.

As an alternative to PDF, you may want to try saving your graphics in the Encapsulated PostScript format. That provides faithful vector graphics, and respects the image's bounding box. It does not handle transparency because PostScript doesn't.

By the way, the export menu spells "Postscript" with lowercase "S". The proper spelling should be with an uppercase S as that is a trademark. That should be corrected.

Simple answer by ChatGPT...

Answered: Rouben Rostamian 8813

July 23 2026

1 8

For whatever it's worth, ChatCPT derives a very simple answer through elementary trigonometry. I don't know how to show ChatCPT's derivation here, so I will just show the final result:

f := (n,x) -> n*x + 2*sum((n-k)/k*sin(k*x), k=1..n-1);

In case you are interested, I obtained that by asking ChatCPT:

integrate 1-cos(n*x))/(1-cos(x)) where n is an integer

An outline of the solution...

Answered: Rouben Rostamian 8813

July 18 2026

6 2

@C_R The following worksheet outlines the technique for solving the initial boundary value problem corresponding to the boundary condition BC3 = u(a,z,t)=0. The insulated boundary condition may be treated similarly.

The implementation of the technique is tedious but straightforward if you are familiar with separation of variables and eigenfunction expansion of solutions of PDEs. janhardo has provided AI-generated details of some concrete cases.

Download heat-pde-outline.mw <-- contains typos, don't use!

Edit:

I spotted a view careless typos in the originally posted worksheet. These have been corrected in this version 2.

An outline of the procedure for

solving heat conduction in a cylinder

>

restart;

We wish to solve the following initial boundary value problem (IBVP)

on a cylinder, with a nonhomogeneous boundary condition at

the cylinder's bottom:

>	PDE := diff(u(r,z,t),t) = alpha(diff(u(r,z,t),r,r) + diff(u(r,z,t),r)/r + diff(u(r,z,t),z,z)); BC1 := -kD[2](u)(r, 0, t) = q(r,t); # bottom BC2 := u(r,h,t) = 0; # top BC3 := u(a,z,t) = 0; # lateral surface IC := u(r,z,0) = 0;

diff(u(r, z, t), t) = alpha*(diff(diff(u(r, z, t), r), r)+(diff(u(r, z, t), r))/r+diff(diff(u(r, z, t), z), z))

We are going to change the dependent variable from to

in such a way that satisfies homogeneous boundary conditions all around
For the sake of consistency of the boundary conditions, we need to have

at all .

But before that, we rearrange BC1 into a more practical equivalent form:

>	BC1_alt := D[2](u)(r,0,t) + 1/k*q(r,t) = 0;

Here is the change of variable from to :

>	trans := u(r,z,t) = v(r,z,t) - 1/kq(r,t)(z-h);

We substitute this expression for into the IBVP and arrive at an IBVP for .

That can be done by hand easily, but here we do it in Maple.

>	PDETools:-dchange(trans, PDE, {v(r,z,t)}, known=q): isolate(%, diff(v(r,z,t),t)): pde := collect(%, [diff(v(r,z,t),r,r), diff(v(r,z,t),r), diff(v(r,z,t),z,z)], simplify);

diff(v(r, z, t), t) = alpha*(diff(diff(v(r, z, t), r), r))+alpha*(diff(v(r, z, t), r))/r+alpha*(diff(diff(v(r, z, t), z), z))+(-z+h)*((diff(diff(q(r, t), r), r))*alpha*r+(diff(q(r, t), r))*alpha-(diff(q(r, t), t))*r)/(k*r)

>	bc1 := PDETools:-dchange(trans, BC1_alt, {v(r,z,t)}, known=q); bc2 := PDETools:-dchange(trans, BC2, {v(r,z,t)}); PDETools:-dchange(trans, BC3, {v(r,z,t)}, known=q): bc3 := subs(q(a,t)=0, %); PDETools:-dchange(trans, IC, {v(r,z,t)}, known=q): ic := isolate(%, v(r,z,0));

We see that all boundary conditions are homogeneous now, which is good.

The price we pay is that the PDE and the initial conditions are nonhomogeneous,

but that's okay. We can handle that.

To simplify this writeup, we write for the rightmost term in the PDE

and therefore write the PDE as

>	pde_alt := diff(v(r,z,t),t) = alpha*(diff(v(r,z,t),r,r) + diff(v(r,z,t),r)/r + diff(v(r,z,t),z,z)) + P(r,z,t);

diff(v(r, z, t), t) = alpha*(diff(diff(v(r, z, t), r), r)+(diff(v(r, z, t), r))/r+diff(diff(v(r, z, t), z), z))+P(r, z, t)

Similarly, we write for the right-hand side of the initial condition, and thus

.

We look for the eigenvalues and eigenfunctions associated with the

partial differential operator

and the associated boundary conditions. The separation of variables leads to

the eigenfunctions where
diff(R(x), x, x)+(diff(R(x), x))/r+(mu/a)^2*R(x) = 0, R(a) = 0
diff(Z(x), x, x)+(gamma/h)^2*Z(x) = 0, (D(Z))(0) = 0, Z(h) = 0 .

each of which has infinitely many solutions:

where is the zeroth order Bessel function, is the th root of ,

and is the th root of the cosine, that is, .

Next, we expand the solution , the "forcing term" , and
and in terms of the eigenfunctions;

>	s1 := v(r,z,t) = Sum(Sum(V__mn(t)R__m(r)Z__n(z), m=1..infinity), n=1..infinity); s2 := P(r,z,t) = Sum(Sum(Psi__mn(t)R__m(r)Z__n(z), m=1..infinity), n=1..infinity); s3 := Q(r,z) = Sum(Sum(Phi__mnR__m(r)Z__n(z), m=1..infinity), n=1..infinity);

v(r, z, t) = Sum(Sum(V__mn(t)*R__m(r)*Z__n(z), m = 1 .. infinity), n = 1 .. infinity)

P(r, z, t) = Sum(Sum(Psi__mn(t)*R__m(r)*Z__n(z), m = 1 .. infinity), n = 1 .. infinity)

Q(r, z) = Sum(Sum(Phi__mn*R__m(r)*Z__n(z), m = 1 .. infinity), n = 1 .. infinity)

Since and are known, the coefficients and may be calculated

via the orthogonality properties of the eigenfunction. We then substitute the

expression for into the PDE and arrive at

>	subs(s1, s2, pde_alt);

Recalling that d^2*Z__n/dz^2 = -(`γ__n`/h)^2*Z__n and d^2*R__m/dr^2+dR__m/(r*dr) = -(`μ__m`/a)^2*R__m , this simplifies to

>	Sum(Sum(diff(V__mn(t),t)R__m(r)Z__n(z), m=1..infinity), n=1..infinity) = alphaSum(Sum((-(mu__m/a)^2V__mn(t)R__m(r)Z__n(z) - (gamma__n/h)^2V__mn(t)R__m(r)Z__n(z)), m=1..infinity), n=1..infinity) + Psi__mn(t)R__m(r)*Z__n(z);

Sum(Sum((diff(V__mn(t), t))*R__m(r)*Z__n(z), m = 1 .. infinity), n = 1 .. infinity) = alpha*(Sum(Sum(-mu__m^2*V__mn(t)*R__m(r)*Z__n(z)/a^2-gamma__n^2*V__mn(t)*R__m(r)*Z__n(z)/h^2, m = 1 .. infinity), n = 1 .. infinity))+Psi__mn(t)*R__m(r)*Z__n(z)

Then, matching the coefficients of the eigenfunctions, we obtain a doubly-infinite set of ODEs

in the unknowns ;

>	odes := diff(V__mn(t),t) = - alpha((mu__m/a)^2 + (gamma__n/h)^2)V__mn(t) + Psi__mn(t);

diff(V__mn(t), t) = -alpha*(mu__m^2/a^2+gamma__n^2/h^2)*V__mn(t)+Psi__mn(t)

With each ODE there is an associated initial condition, obtained by setting

in s1 and equating the result to s3, whereby we obtain

>	ics := V__mn(0) = Phi__mn;

The quantities and are know as pointed out before, and therefore

the odes + ics may be solved for each and . That may done analytically if

(and therefore ) is an elementary function, else it would be done

numerically.

With thus computed, we obtain from s1, and then

from ,

as introduced near the beginning of this worksheet.

Download heat-pde-outline-ver2.mw

Perhaps this is all you need...

Answered: Rouben Rostamian 8813

July 15 2026

3 1

A thorough answer will depend on what you intend to do with the result. Based on what you have shown so far, the following should serve your needs.

Download mw.mw

Edit: Alternatively, you may simply do

restart;
e := Matrix(3,3, (i,j) -> 1/2*(diff(u[i](x[1],x[2],x[3]),x[j])
                               + diff(u[j](x[1],x[2],x[3]),x[i])));

Specify the file name as ./data.xlsx...

Answered: Rouben Rostamian 8813

July 12 2026

1 1

To do what you want, specify the file name as ./data.xlsx, as in

restart;
with(ExcelTools):
Import("./data.xlsx");

This was tested and works as expected in Linux. It should work in other operating systems as well.

The shooting method...

Answered: Rouben Rostamian 8813

July 05 2026

2 0

You may solve a set of ODEs interfacing over any number of adjacent intervals through the method of shooting. The idea is that you have the value of the solution at the left end, you leave the value of the slope there as an arbitrary parameter, and then apply fsolve() to determine the value of the parameter so that the solution hits a desired value at the right end.

The following worksheet illustrates how to solve such a problem over two adjacent intervals. It may be extended in the obvious way to solve the interface problem over any number of adjacent intervals.

See if you can adapt this to solve your problem.

Solving an interface problem through

the shooting method

>

restart;

>	with(plots):

We solve two ODEs over adjacent intervals. The valules of the

solutions are given at the leftmost and rightmost boundaries.

At the interface of the two intervals the solution and its derivative

are continuous.

With obvious changes, this worksheet may be extended to

solve such interface problems over any number of intervals.

>	de1 := diff(u(x),x,x) + sin(u(x)) = 1;

diff(diff(u(x), x), x)+sin(u(x)) = 1

>	de2 := diff(v(x),x,x) + diff(v(x),x) + v(x) = 2;

diff(diff(v(x), x), x)+diff(v(x), x)+v(x) = 2

The domain consists of the adjacent intervals C[1]..C[2] and C[2]..C[3].
Change as needed.

>	C := [1,2,4];

Boundary conditions at the leftmost and rightmost boundaries.

Change as needed.

>	bc_left := 1; bc_right := -1;

The shooting procedure.

1. Shoot from the left boundary with and with a

given initial slope .

2. Solve the first ODE over C[1]..C[2], find the values of and at C[2],

and use those values as initial data to solve the second ODE over C[2]..C[3].

3. Return the value of the at =C[3].

>

right_val := proc(a)
        local uval, udval;
        if not type (a, realcons) then return 'procname'(args) end if;
        dsolve({de1, u(C[1])=bc_left, D(u)(C[1])=a}, numeric);
        uval, udval := rhs~(%(C[2]))[2..3][]; # u and u' at C[2]
        dsolve({de2, v(C[2])=uval, D(v)(C[2])=udval}, numeric);
        subs(%(C[3]), v(x));
end proc:

Then the following (for illustration only) finds the value of the initial slope that
shoots to hit the given boundary condition at C[3].

>	fsolve(right_val(a) = bc_right);

Procedure finds the initial slope and then solves the ODEs, and plots.

>

do_plot := proc()
        local s, dsol1, dsol2, uval, udval;
        s := fsolve(right_val(a) = bc_right);
        dsol1 := dsolve({de1, u(C[1])=bc_left, D(u)(C[1])=s}, numeric);
        uval, udval := rhs~(dsol1(C[2]))[2..3][];
        dsol2 := dsolve({de2, v(C[2])=uval, D(v)(C[2])=udval}, numeric);
        display(
                [odeplot(dsol1, x=C[1]..C[2]),
                 odeplot(dsol2, x=C[2]..C[3])],
                color=[red,blue]);
end proc:

>	do_plot();

Download shooting.mw

Or slightly shorter...

Answered: Rouben Rostamian 8813

July 04 2026

1 0

ithprime(b) $b=1..19;

A partial answer...

Answered: Rouben Rostamian 8813

June 21 2026

2 2

After a bit of algebra it is possible to reduce the original integral to the following but I don't know how to evaluate the result symbolically. If anyone can figure out how, then I will post my calculations that lead to that result, otherwise it's not worth the trouble.

Download mw.mw

An animation...

Answered: Rouben Rostamian 8813

June 17 2026

2 3

Here an implementation of the equations of motion derived in
https://www.youtube.com/watch?v=1h4zlDRjV0c
but you may want to have a look at the older Maple demos that the other posters have noted. They are more illuminating.

>

restart;

>	with(plots):

>	with(plottools):

Position and velocity of the small mass

>	A := <a[1](t), a[2](t), a[3](t)>; U := <D(a[1])(t), D(a[2])(t), D(a[3])(t)>;

Position and velocity of the large mass

>	B := <b[1](t), b[2](t), b[3](t)>: V := <D(b[1])(t), D(b[2])(t), D(b[3])(t)>:

Differential equations of motion

>	de1 := diff(A,t,t) = -1/R^2(U^+.U)A - 2m[2]/(m[1]+m[2])/R^2(U^+.V)B: de2 := diff(B,t,t) = -1/R^2(V^+.V)B - 2m[1]/(m[1]+m[2])/R^2(U^+.V)A:

Initial conditions:

>	e := 0.01: A =~ < 0, e, sqrt(R^2 - e^2)>, U =~ < 0.01, 0, 0 >, B =~ < R, 0, 0>, V =~ < 0, 1, 0 >: subs(t=0, [%]): ic := map(seq, %)[];

a[1](0) = 0, a[2](0) = 0.1e-1, a[3](0) = (R^2-0.1e-3)^(1/2), (D(a[1]))(0) = 0.1e-1, (D(a[2]))(0) = 0, (D(a[3]))(0) = 0, b[1](0) = R, b[2](0) = 0, b[3](0) = 0, (D(b[1]))(0) = 0, (D(b[2]))(0) = 1, (D(b[3]))(0) = 0

Set the parameter values and solve the system of differential equations

>	subs(R=1, m[1]=2, m[2]=3, {de1, de2, ic}): dsol := dsolve(%, numeric, output=operator):

Proc to produce animation frame at time t:

>

frame := proc(t)
        local A, B;
        if not type (t, realcons) then return 'procname'(args) end if;
        A := subs(dsol, [a[1], a[2], a[3]]);
        B := subs(dsol, [b[1], b[2], b[3]]);
        display(
                line(-A(t), A(t), color="Red"),
                line(-B(t), B(t), color="Green"),
                pointplot3d([A(t),-A(t)], symbol=solidsphere,
                        symbolsize=25, color="Red"),
                pointplot3d([B(t),-B(t)], symbol=solidsphere,
                        symbolsize=40, color="Green"),
                spacecurve(A(s), s=0..t, linestyle=dot, color="Red"),
                spacecurve(B(s), s=0..t, linestyle=dot, color="Green"),
        scaling=constrained, thickness=4);
end proc:

Animate:

>	st := time(): animate(frame, [t], t=0..50, frames=300); time() - st;

Download: animation.mw

Here's how I would do it...

Answered: Rouben Rostamian 8813

June 01 2026

3 4

>

restart;

>	with(plots):

>	points := [3,6], [1,0], [4,2], [5,1];

Find the min and max of the coordinates

>	xmin, xmax := min(seq(p[1], p in points)), max(seq(p[1], p in points)); ymin, ymax := min(seq(p[2], p in points)), max(seq(p[2], p in points));

>	display( plottools:-polygon([points], color="Honeydew", thickness=2), pointplot([seq(seq([i,j], i=xmin..xmax), j=ymin..ymax)], symbol=circle), pointplot([points], symbol=solidcircle, color="Red"), symbolsize=20, scaling=constrained, axes=none);

Download grid.mw

This may be what you want...

Answered: Rouben Rostamian 8813

May 02 2026

2 3

restart;

>	with(plots):

>	Z := n -> [Re, Im](8(1/4 + 1/4I*sqrt(3))^n);

proc (n) options operator, arrow; ([Re, Im])(8*(1/4+((1/4)*I)*sqrt(3))^n) end proc

>	Triangle := i -> plottools:-polygon([[0,0], Z(i-1), Z(i)]);

>	frame := proc(n::posint) local i, colors; colors := seq(cat("OldPlots ", i), i=1..6); seq(Triangle(i), i=1..n); display(%, color=[colors, colors]); end proc:

>	display(seq(frame(i), i=1..6), insequence, scaling=constrained, size=[800,400]);

mw.mw

An alternative solution...

Answered: Rouben Rostamian 8813

April 16 2026

3 3

This is an alternative to the solution presented by dharr.

>

restart;

>	Typesetting:-Settings(typesetprime=true):

>	kernelopts(version);

(1)

Let D be a domain enclosed by the closed curve C in the Cartesian xy plane.

The curve goes through the origin and its length L is prescribed. We wish

to find D that maximizes the integral .

We expect C to be tear-drop shaped, symmetric about the x axis, with a pointy

vertex at the origin and smooth otherwise. The curve C meets the x axis at

the far right extreme of D at some point x=q, where q is to be determined.

Because of the presumed symmetry and smoothness, the slope of C at x=q

will be infinite.

Let the upper half of the boundary C be the graph of the function

y(x), 0 < x < q. We have

.

By Fubini's theorem, the integral may be expressed as
I = 2*(int(x*y(x), x = 0 .. q)) .

Moreover, since the overall length of C is L, we have

2*(int(sqrt(1+(diff(y(x), x))^2), x = 0 .. q)) = L .

Thus, finding the maximizing y(x) amounts to finding the extrema of

the functional

J(y, q) = 2*(int(x*y(x)-lambda*sqrt(1+(diff(y(x), x))^2), x = 0 .. q)) ,

where is the Lagrange multiplier that enforces the length constraint.

Letting , we express compactly as

"J(y,q) = 2 (∫)[0]^(q)F(x,y,y') &DifferentialD;x."

>	F := xy(x) - lambdasqrt(1 + diff(y(x),x)^2);

(2)

The corresponding Euler-Lagrange equation is

which evaluates to

>	Diff(Physics:-diff(F, diff(y(x),x)), x) - Physics:-diff(F, y(x)) = 0;

Diff(-lambda*(diff(y(x), x))/(1+(diff(y(x), x))^2)^(1/2), x)-x = 0

(3)

We divide the equation by and rearrange:

>	Diff(diff(y(x), x)/sqrt(1 + diff(y(x), x)^2), x) = - x/lambda;

Diff((diff(y(x), x))/(1+(diff(y(x), x))^2)^(1/2), x) = -x/lambda

(4)

Then, we integrate with respect to x

>	int(lhs((4)), x) = C + int(rhs((4)), x);

(diff(y(x), x))/(1+(diff(y(x), x))^2)^(1/2) = C-(1/2)*x^2/lambda

(5)

where C is the integration constant.

To further simplify, observe that the the fraction (diff(y(x), x))/sqrt(1+(diff(y(x), x))^2) takes values between and , and therefore
we may take it to be the sine of some angle , as in sin(theta) = (diff(y(x), x))/sqrt(1+(diff(y(x), x))^2), `and`(-(1/2)*Pi <= theta, theta <= (1/2)*Pi)

Then from (5) we get

>	sin(theta) = C - x^2/(2*lambda);

sin(theta) = C-(1/2)*x^2/lambda

(6)

Moreover, sin*theta = (diff(y(x), x))/sqrt(1+(diff(y(x), x))^2) implies that . Let's make a note of it:

>	diff(y(x),x) = tan(theta);

(7)

In particular, we know that , therefore , and consequently

Plugging this into (6) we get -1 = C-q^2/(2*lambda) , and therefore

C = q^2/(2*lambda)-1 . We plug this back into (6) and arrive at
sin*theta = q^2/(2*lambda)-1-x^2/(2*lambda) and q^2/(2*lambda)-1-x^2/(2*lambda) = q^2*(1-(x/q)^2)/(2*lambda)-1. We introduce the abbreviations
xi = x/q, alpha = q^2/(2*lambda) , whereby

>	sin(theta) = alpha*(1-xi^2) - 1;

(8)

Now we calculate the length of the boundary. We need to integrate .

We apply (7) to change the variable from to and then apply (8) to change to

the variable;

>	sqrt(1 + diff(y(x),x)^2); subs((7), %); simplify(%) assuming theta > -Pi/2, theta < Pi/2: convert(%, sincos): subs(cos(theta) = sqrt(1 - sin(theta)^2), %); subs((8), %);

1/(1-sin(theta)^2)^(1/2)

1/(1-(alpha*(-xi^2+1)-1)^2)^(1/2)

(9)

Calculating the boundary's length requires integrating (9) from to .

The result is a rather unpleasant-looking function of :

>	int((9), xi=0..1);

(10)

Let's plot (10) to see what we have;

>	plot((10), alpha=-4..4);

We see that (10) takes on real values only where .

From now on we limit to that range. Within that range, the expression

simplifies considerably:

>	simplify((10)) assuming alpha > 0, alpha < 2;

(1/2)*2^(1/2)*EllipticK((1/2)*2^(1/2)*alpha^(1/2))/alpha^(1/2)

(11)

The total length L of the domain's boundary is given by
L = 2*(int(sqrt(1+(diff(y(x), x))^2), x = 0 .. q)) and 2*(int(sqrt(1+(diff(y(x), x))^2), x = 0 .. q)) = 2*q*(int(1/(1-(alpha*(-xi^2+1)-1)^2)^(1/2), xi = 0 .. 1)) , therefore

>	L = 2q(11); isolate(%, q);

L = q*2^(1/2)*EllipticK((1/2)*2^(1/2)*alpha^(1/2))/alpha^(1/2)

q = (1/2)*L*2^(1/2)*alpha^(1/2)/EllipticK((1/2)*2^(1/2)*alpha^(1/2))

(12)

That establishes a relationship between and . We will need that later. For now,

let's calculate the value of y(x) at x=q, for a given. We call that .

Since , and , we have

>	y(q,alpha) = Int(tan(theta), x=0..q); subs(tan(theta)=sin(theta)/sqrt(1-sin(theta)^2), %); subs((8), %);

y(q, alpha) = Int(tan(theta), x = 0 .. q)

y(q, alpha) = Int(sin(theta)/(1-sin(theta)^2)^(1/2), x = 0 .. q)

y(q, alpha) = Int((alpha*(-xi^2+1)-1)/(1-(alpha*(-xi^2+1)-1)^2)^(1/2), x = 0 .. q)

(13)

where , We change the integration variable from to and
and let Maple do the integration:

>	y(q,alpha) = q*subs(q=1, x=xi, rhs((13))); value(%) assuming alpha>0, alpha<2: simplify(%): lhs(%) = subs((12), rhs(%));

y(q, alpha) = q*(Int((alpha*(-xi^2+1)-1)/(1-(alpha*(-xi^2+1)-1)^2)^(1/2), xi = 0 .. 1))

y(q, alpha) = (1/2)*L*(EllipticK((1/2)*2^(1/2)*alpha^(1/2))-2*EllipticE((1/2)*2^(1/2)*alpha^(1/2)))/EllipticK((1/2)*2^(1/2)*alpha^(1/2))

(14)

The value of at the right extreme, , should be zero. Setting the above expression

to zero gives us the value of :

>	my_alpha := fsolve(rhs((14))/L);

(15)

Once we have , we may calculate :

>	my_q := evalf(subs(alpha=my_alpha, rhs((12))));

(16)

Now let us calculate the shape of the curve. We have diff(y(x), x) = tan*theta and tan*theta = sin*theta*(1/(cos*theta)) and sin*theta*(1/(cos*theta)) = sin*theta/sqrt(-sin^2*theta+1) ,

which then may be evaluated according to (8). We get

>	sin(theta)/sqrt(1-sin(theta)^2); subs((8), %): subs(xi=x/q, %): diff(y(x),x) = %;

sin(theta)/(1-sin(theta)^2)^(1/2)

diff(y(x), x) = (alpha*(-x^2/q^2+1)-1)/(1-(alpha*(-x^2/q^2+1)-1)^2)^(1/2)

(17)

We are almost there. Finding the shape of the upper boundary, y(x), is

a matter of integrating the above:

>	int(rhs((17)), x) assuming q>0, alpha>0, alpha<2: ans := simplify(%);

(((q^2-x^2)*alpha-2*q^2)/(q^2*(alpha-2)))^(1/2)*((alpha-2)*EllipticE(x/q, alpha^(1/2)/(alpha-2)^(1/2))+EllipticF(x/q, alpha^(1/2)/(alpha-2)^(1/2)))*q^3*((q^2-x^2)/q^2)^(1/2)/(-(q^2-x^2)^2*alpha^2+(2*q^4-2*q^2*x^2)*alpha)^(1/2)

(18)

Let's verify that the solution determined above satisfies the initial condition :

>	simplify(subs(x=0, ans));

(19)

We substitute for the previously found numerical values of and

and the desired value for the overall length of the boundary, and plot the result

>

L := 100;

(20)

>	subs(alpha=my_alpha, q=my_q, ans): plot([%,-%], x=0..my_q, scaling=constrained);

Finally, we find the maximum value of objective function. That's a matter of

evaluating 2*(int(x*y(x), x = 0 .. q)) . Unfortunately there seems to be no symbolic

expression for that, so we do it numerically:

>	subs(alpha=my_alpha, q=my_q, ans): 2int(x%, x=0..my_q, numeric): Re(%);

(21)

>

Download TearDrop.mw

This ought to do it;plot3d(1, t = -Pi .....

Answered: Rouben Rostamian 8813

April 10 2026

1 8

This ought to do it;

plot3d(1, t = -Pi .. Pi, p = 0 .. Pi, coords = spherical, grid=[25,11]);

You wrote: ...want output of all numbers...

Answered: Rouben Rostamian 8813

April 04 2026

1 0

You wrote:
...want output of all numbers less than 234 that are divisible by eleven.

This will to do that:

seq(11*n, n=1..21);

The output is
11, 22, 33, 44, 55, 66, 77, 88, 99, 110, 121, 132, 143, 154, 165, 176, 187, 198, 209, 220, 231

Aside: Your subject line refers to "prime number set". But that does not have anything to do with the question. Perhaps you meant to ask something else?

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

Ask a Question

Create a Post

8813 Reputation

23 Badges

MaplePrimes Activity

These are answers submitted by Rouben Rostamian

Sliding with friction...

Unfortunately not possible...

Simple answer by ChatGPT...

An outline of the solution...

Perhaps this is all you need...

Specify the file name as ./data.xlsx...

The shooting method...

Or slightly shorter...

A partial answer...

An animation...

Here's how I would do it...

This may be what you want...

An alternative solution...

This ought to do it;plot3d(1, t = -Pi .....

You wrote: ...want output of all numbers...

Save this setting as your default sorting preference?

Ask a Question

Create a Post

Generating PDF…

Save this setting as your default sorting preference?
Note: You can change your preference any time in your account settings
Don't show this again