Rouben Rostamian

Simple answer by ChatGPT...

Answered: Rouben Rostamian 8776

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 8776

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 8776

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 8776

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 8776

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 8776

July 04 2026

1 0

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

A partial answer...

Answered: Rouben Rostamian 8776

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 8776

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 8776

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 8776

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 8776

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 8776

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 8776

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?

Julian Day conversion...

Answered: Rouben Rostamian 8776

April 04 2026

2 1

This is a Maple version of a C code that I obtained from
Ian Miller <ian_m@cix.compulink.co.uk>

some 30 years ago.

Proc receives a Julian Day and determines the corresponding

year, month, day

>

invertjd := proc(jd::posint)
        local tmp1, tmp2, y, m, d;
        tmp1 := jd + 68569;
        tmp2 := iquo(4*tmp1,146097);
        tmp1 := tmp1 - iquo(146097 * tmp2 + 3, 4);
        y := iquo( 4000 * (tmp1 + 1) , 1461001);
        tmp1 := tmp1 - iquo(1461 * y, 4) + 31;
        m := iquo(80 * tmp1 , 2447);
        d := tmp1 - iquo(2447 * m , 80);
        tmp1 := iquo(m , 11);
        m := m + 2 - 12 * tmp1;
        y := 100 * (tmp2 - 49) + y + tmp1;
        return (y, m, d);
end proc: