Advanced Search

Rouben Rostamian

8676 Reputation

23 Badges

19 years, 340 days
Contact Rouben Rostamian

MaplePrimes Activity


These are answers submitted by Rouben Rostamian

Here's how I would do it...

Rouben Rostamian 8676
5 hours ago
3 1
> 

restart;
> 

with(plots):
> 

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

[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));

1, 5

0, 6

> 

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...

Rouben Rostamian 8676
May 02 2026
2 3
 

restart;
> 

with(plots):
> 

Z := n -> [Re, Im](8*(1/4 + 1/4*I*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)]);

proc (i) options operator, arrow; plottools:-polygon([[0, 0], Z(i-1), Z(i)]) end proc

> 

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...

Rouben Rostamian 8676
April 16 2026
3 3

This is an alternative to the solution presented by dharr.

> 

restart;
> 

Typesetting:-Settings(typesetprime=true):
> 

kernelopts(version);

`Maple 2024.2, X86 64 LINUX, Oct 29 2024, Build ID 1872373`

 (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 I = int(`∫D`*x, A).

 

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

y(0) = 0, y(q) = 0, (D(y))(q) = -infinity.

 

By Fubini's theorem, the integral I 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 lambda is the Lagrange multiplier that enforces the length constraint.

Letting F(x, y(x), diff(y(x), x)) = x*y(x)-lambda*sqrt(1+(diff(y(x), x))^2), we express J compactly as

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

> 

F := x*y(x) - lambda*sqrt(1 + diff(y(x),x)^2);

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

 (2)

NULLThe corresponding Euler-Lagrange equation is

"(&DifferentialD;)/(&DifferentialD; x)((&PartialD; F)/(&PartialD; y')) -(&PartialD; F)/(&PartialD; y)=0, "
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 -lambda 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 -1 and 1, and therefore
we may take it to be the sine of some angle theta, 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 diff(y(x), x) = tan*theta.  Let's make a note of it:NULL

> 

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

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

 (7)

In particular, we know that (D(y))(q) = -infinity, therefore theta(q) = -(1/2)*Pi, and consequently

sin(theta(q)) = -1. 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;

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

 (8)

Now we calculate the length of the boundary.  We need to integrate sqrt(1+(diff(y(x), x))^2)NULL.

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

the xi 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+(diff(y(x), x))^2)^(1/2)

  

(1+tan(theta)^2)^(1/2)

  

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 xi = 0 to xi = 1.

The result is a rather unpleasant-looking function of alpha:

> 

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

sqrt(-1/(alpha-2))*EllipticK(sqrt(alpha/(alpha-2)))/sqrt(alpha)+piecewise(2 < alpha, sqrt(1/alpha)*EllipticF(sqrt((alpha-2)/alpha), sqrt(alpha/(alpha-2)))*(-sqrt(2)*sqrt(-2*alpha*sqrt((alpha-2)/alpha)/(alpha-2))*sqrt(-alpha*sqrt((alpha-2)/alpha))+2*sqrt(alpha*sqrt((alpha-2)/alpha)/(alpha-2))*sqrt(alpha*sqrt((alpha-2)/alpha)))/(2*sqrt(-alpha*sqrt((alpha-2)/alpha))*sqrt(alpha*sqrt((alpha-2)/alpha))), 0)

 (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 0 < alpha and alpha < 2.

From now on we limit alpha 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 = 2*q*(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 alpha and q.  We will need that later.  For now,

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

 Since diff(y(x), x) = tan*theta, and y(0) = 0, 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 xi = x/q,    We change the integration variable from x to xi 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 y at the right extreme, x = q, should be zero.  Setting the above expression

to zero gives us the value of alpha:

> 

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

1.652229532

 (15)

Once we have alpha, we may calculate q:

> 

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

.3915937454*L

 (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 y(0) = 0:

> 

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

0

 (19)

We substitute for the previously found numerical values of alpha and "q,"

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

> 

L := 100;

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):
2*int(x*%, x=0..my_q, numeric):
Re(%);

15468.55913

 (21)
> 

 

Download TearDrop.mw

 

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

Rouben Rostamian 8676
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...

Rouben Rostamian 8676
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...

Rouben Rostamian 8676
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:

Examle: Today's Julian Day is 2461135.   So we get

> 

invertjd(2461135);

2026, 4, 4

  

That's April 4, 2026.

 

Download julianday.mw

Here is what I get...

Rouben Rostamian 8676
March 15 2026
2 1
> 

restart;
> 

kernelopts(version); # sorry, can't stand the newfangled interface of 2025 and 2026

`Maple 2024.2, X86 64 LINUX, Oct 29 2024, Build ID 1872373`

> 

ode:=x*diff(y(x),x) = y(x)*cos(ln(y(x)/x));

x*(diff(y(x), x)) = y(x)*cos(ln(y(x)/x))

The general solution

> 

dsol := dsolve(ode) assuming positive;

y(x) = exp(2*arctan(1/(ln(x)+c__1)))*x

Singular solution:

> 

limit(dsol, c__1=infinity);

y(x) = x

Download mw.mw

 

A step-by-step method...

Rouben Rostamian 8676
December 18 2025
2 1

This does what you want in Maple 2024:

x := sqrt(sqrt(12)/(9 + sqrt(108)));
rationalize(x^2);
simplify(%);
sqrt(%);

Perhaps others can come up with a shorter way of doing that.

Missing data...

Rouben Rostamian 8676
November 08 2025
0 0

Execute the command indets([XDD], name) to see the names of unassigned parameters in your system of differential equations.  That command returns {g, q, r__rel, rho, t}.  The t is okay since that is the time varialble.  The other four items are not.  You need to assign values to them.

Additionally, you have gamma(t) as one of your unknowns, but gamma is a predefined constant in Maple.  Use some other name for it.

Correct your worksheet accordingly and write again if you run into other problems.

Use isolve()...

Rouben Rostamian 8676
October 25 2025
2 1

Let 19*ab be the birth year where a and b are single digits.

Then the age at the birthday party is 1968-1900-10*a-b,

which according to the hypothesis equals 1+9+a+b.

To find a and b, we do

> 

1 + 9 + a + b = 1968 - (1900 + 10* a + b);
sol := isolve(%);

10+a+b = 68-10*a-b

{a = 2*_Z1, b = 29-11*_Z1}

were _Z1 is an arbitrary integer.  We see that b will be a

positive single digit integer only if _Z1 is 2.  Therefore

> 

subs(_Z1=2, sol);

{a = 4, b = 7}

We conclude that the birth year is 1947, and therefore in

2025 the age of the "child" will be

> 

2025 - 1947;

78


 

Download mw.mw

Animation...

Rouben Rostamian 8676
August 31 2025
0 0

 

> 

restart;
> 

with(plots):
> 

with(plottools):
> 

F := (x, y) -> (3*surd(4 - (abs(x) - 2)^2, 2) + abs(x) - 20 - 4*y)*(3*surd(4 - (abs(x) - 6)^2, 2) + abs(x) - 20 - 4*y)*(x^2 + 4*y^2 - 100*sqrt(abs(7 - abs(2*y - 1))/(7 - abs(2*y - 1))))*(2*(abs(x) - 3)^2 - 9*y + 18*sqrt(abs(2 - abs(abs(x) - 4))/(2 - abs(abs(x) - 4))))*(-68*abs(abs(x) - 3/2) - 9*y + 54*sqrt(abs(43 - abs(136*abs(x) - 229))/(43 - abs(136*abs(x) - 229))))*(y - 5*sqrt(abs(1 - abs(x))/(1 - abs(x))));

proc (x, y) options operator, arrow; (3*surd(4-(abs(x)-2)^2, 2)+abs(x)-20-4*y)*(3*surd(4-(abs(x)-6)^2, 2)+abs(x)-20-4*y)*(x^2+4*y^2-100*sqrt(abs(7-abs(2*y-1))/(7-abs(2*y-1))))*(2*(abs(x)-3)^2-9*y+18*sqrt(abs(2-abs(abs(x)-4))/(2-abs(abs(x)-4))))*(-68*abs(abs(x)-3/2)-9*y+54*sqrt(abs(43-abs(136*abs(x)-229))/(43-abs(136*abs(x)-229))))*(y-5*sqrt(abs(1-abs(x))/(1-abs(x)))) end proc

> 

p := implicitplot(F(x, y) = 0, x = -10 .. 10, y = -10 .. 10, grid=[40,40], gridrefine=4):
> 

animate(scale, [p, 3/4 + 1/4*sin(t), 1], t=0..2*Pi, frames=15, scaling=constrained);

 

Download bat.mw

 

Try this...

Rouben Rostamian 8676
August 02 2025
1 1

The limitation may be imposed by the Unix shell, not Maple.  If your shell is bash, which seems to be the default in most Linux distributions, then "ulimit" controls the maximum available resources.  You may try the command "ulimit unlimited" in bash to remove caps on resources before starting Maple.  Look up ulimit in bash's man page for details. I assume that you are starting Maple from the command-line.

If your shell is tcsh (like mine), then the corresponding command is "unlimit".  Again, the man page has the details.

I assume you want something like this:pl...

Rouben Rostamian 8676
July 30 2025
0 0

I suppose you want something like this:

plots:-animate(plot, [[y1(1)+r1*cos(t),y1(2)+r1*sin(t),t=0..T]],
    T=0..2*Pi,color=blue, scaling=constrained);

If not, then you should be more explicit about what does not work.

Maple's solution...

Rouben Rostamian 8676
July 22 2025
3 6

@dharr has shown how to implement your textbook's approach to solving the PDE.  I just want to remark that a solution may be obtained directy in Maple through the traditional apprach.  Maple's solution appears to be computationally more efficient since it involves an integration over the interval [0,t] as opposed to the textbook's solution which integrates over [0,infinity].

Plotting the result at the default 25x25 grid size takes about 10 minutes on my computer.  On a 10x10 grid it takes about 20 seconds.

> 

restart;
> 

pde := diff(u(x,t),t) = diff(u(x,t),x,x);

diff(u(x, t), t) = diff(diff(u(x, t), x), x)

  
> 

ic := u(x,0) = exp(-x);

u(x, 0) = exp(-x)

  
> 

bc := u(0,t) = cos(t);

u(0, t) = cos(t)

  
> 

pdsol := pdsolve({pde,ic,bc}) assuming x > 0, t > 0;

u(x, t) = (1/2)*((1/4)*x*(Int(-sin(-t+zeta)*exp(-(1/4)*x^2/zeta)*(x^2-6*zeta)/zeta^(7/2), zeta = 0 .. t))+Pi^(1/2)*(erf((1/2)*(2*t+x)/t^(1/2))*exp(x+t)-erf((1/2)*(2*t-x)/t^(1/2))*exp(-x+t)-exp(x+t)+exp(-x+t)))/Pi^(1/2)

  
> 

plot3d(rhs(pdsol), x=1e-3..3, t=0..2*Pi, color="Orange");

 
 

 

Download heat-eq2.mw

 

Plotting arrowheads...

Rouben Rostamian 8676
June 28 2025
2 0

Here is a general proc for plotting arrowheads. 

> 

restart;
> 

with(plots):


Proc arrowhead()
        Draws an arrowhead at the point p in Cartesian coordinates.
        The arrowhead's reference point is its reentrant vertex.
Arguments:
        p: is the arrowhead's location specified as a list [a, b]  or a vector `<,>`(a, b) .
        angle: is the angle of the arrowhead (rotation) relative to the horizontal axis
        mag: (optional, default=1) is the magnification factor.  With the default mag = 1,
        the distant between the reentrant corner and the arrow tip is 1/10.

        Any additional arguments are assumed to be plotting options and are passed
        to plots:-display.

>

 arrowhead := proc(p::{[realcons,realcons],'Vector'(2,realcons)}, angle, {mag:=1})
        uses plots, plottools;
        polygon([[0,0], [-1/2,-sqrt(3)/2], [2,0], [-1/2,sqrt(3)/2]],
                style=polygon);
        rotate(%, angle);
        scale(%, mag/20, mag/20);
        translate(%, p[1], p[2]);
        display(%, _rest);
end proc:
   

 

Example

> 

arrowhead(<3,4>, Pi/3, mag=10, color=red, scaling=constrained);

 

Example

> 

A, B, C := <-1,0>, <0,0>, <1,0>;

Vector[column](%id = 36893614961805893492), Vector[column](%id = 36893614961805893628), Vector[column](%id = 36893614961805893748)

 (1)
> 

display(
        pointplot([A, B, C], symbol=solidcircle, symbolsize=25),
        plottools:-line([-2,0],[2,0]),
        arrowhead((A+B)/2,   Pi, mag=1.2),
        arrowhead((B+C)/2,   0,  mag=1.2),
        arrowhead(A-<1/2,0>, 0,  mag=1.2),
        arrowhead(C+<1/2,0>, Pi, mag=1.2),
scaling=constrained, tickmarks=[0,0], size=[800,100]);
 
   

 

Download arrowhead.mw

1 2 3 4 5 6 7 Last Page 1 of 58