Advanced Search

Rouben Rostamian

7711 Reputation

22 Badges

19 years, 2 days
Contact Rouben Rostamian

MaplePrimes Activity


These are answers submitted by Rouben Rostamian

Take y to be a function of x.  Then...

Rouben Rostamian 7711
June 15 2019
5 0

Take y to be a function of x. Then differentiate the equation once to get rid of B. Then isolate A, and then differentiate again to get rid of A. The resulting differential equation will not have A or B in it. Here is how.

> 

restart;
> 

t1 := x^2 + A*y(x)^2 = B;

x^2+A*y(x)^2 = B

> 

t2 := diff(t1, x);

2*x+2*A*y(x)*(diff(y(x), x)) = 0

> 

t3 := isolate(t2, A);

A = -x/(y(x)*(diff(y(x), x)))

> 

t4 := diff(t3, x);

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

That's the differential equation you are looking for, but let's simplify it:

> 

de := factor(isolate(t4, diff(y(x),x,x)));

diff(diff(y(x), x), x) = -(x*(diff(y(x), x))-y(x))*(diff(y(x), x))/(y(x)*x)

OK.  Let's see if Maple can recover the original family of curves from this:

> 

dsolve(de, y(x), implicit);

_C1*x^2-_C2+(1/2)*y(x)^2 = 0

Yes!  That's equivalent to the original equation, with some re-named constants.


 

Download mw.mw

PS: The English verb for computing a derivative is "to differentiate", not "to derive", but don't feel bad for having mixed up the two.  Some native English speakers make the same mistake.

 

Geodesics...

Rouben Rostamian 7711
June 12 2019
4 3

Some time ago when you asked about a related question, I made worksheets for calculating geodesics on all five regular polyhedra but didn't post them because I felt that there ought to be a better approach than what I have right now.  Now that you are asking again, let me show how I do geodesics on a cube.  I have attempted to document the approach in cube.mw in detail, but let me know if you have questions.  I have not documented the work on the other four regular polyhedra so I am not posting those worksheets.  I do, however, post a sample of a geodesic on a dodecahedron just to show that I have done it.

 

Download worksheet: cube.mw

 

Animation of the Geneva Drive...

Rouben Rostamian 7711
June 07 2019
3 0

I saw this question a few days ago and was unable to read the original undecipherable code, so I wrote my own code from scratch which I did not post then since I was waiting for an update from the original poster.

Checking the thread today I see that there has been quite a bit of discussion in the meantime, so perhaps this is a good time to post my code which is entirely independent of that presented in the discussion.

Here are a few sample animations for your enjoyment!  See Geneva Drive in Wikipedia for more info.

 

Download geneva.mw

 

Missing colon...

Rouben Rostamian 7711
June 04 2019
1 3

Insert missing colon in
if flag=1 then  dumbindex:=dumbindex+1  xdumb[dumbindex]:=x[i]: ...
to make it
if flag=1 then  dumbindex:=dumbindex+1:  xdumb[dumbindex]:=x[i]: ...

Solution with fixes...

Rouben Rostamian 7711
June 02 2019
2 6
> 

restart;
> 

with(plots):
> 

with(CurveFitting):
> 

pts := <0, 3.016315859, 5.753536424, 8.489713643, 11.22797756, 13.96519812, 16.09414745, 18.83136802, 21.26445296, 24.00167353, 26.73889409, 29.47611466, 31.60506399, 34.34228455, 37.07950512, 39.51259007, 41.6415394, 44.37875996, 46.81184491, 49.54906547, 51.6780148, 54.11109975, 56.84832032, 58.97726964, 62.32276145, 64.45171078, 66.88479572, 69.01374505, 71.44578665, 74.1814422, 76.91814109, 79.0455254, 81.47704533, 84.21165752, 86.33695514, 88.76691004, 90.89116431, 93.3195542, 96.04894966, 98.17059556, 100.5958554, 103.3231642, 105.4432451, 107.8679832, 110.5947703, 113.3236441, 116.0514746, 118.1746855, 120.9066893, 125.4640285, 127.8960701, 130.634334, 132.763805, 135.1994984, 137.9382839, 140.06932, 142.5034483, 145.2411905, 147.976846, 150.7119799, 153.4450271, 156.1770309, 158.9069481, 161.9399575, 164.6693529, 167.0951345, 169.2178237, 171.9472192, 174.3756091, 176.4998633, 179.2334322, 181.9680444, 184.4011293, 186.531122, 189.2719943, 192.0149532, 194.4542983, 197.2014306, 201.7822452, 203.921628, 207.5889923, 210.3397763, 213.699875, 217.3651526, 220.1133283, 222.859939, 224.9961918, 227.7412375, 230.1795391, 232.3131835, 234.7504419, 237.4913141, 240.5368437, 242.9720153, 245.7108009, 247.8313153, 247.8413153, 247.8513153, 250>,
<-0.041150436, -0.041150436, -0.04113729, -0.038722776, -0.041110998, -0.041097851, -0.041087627, -0.04107448, -0.041062795, -0.041049649, -0.041036503, -0.041023356, -0.041013132, -0.040999985, -0.040986839, -0.040975154, -0.040964929, -0.040951783, -0.040940097, -0.040926951, -0.040916726, -0.040905041, -0.040891895, -0.04088167, -0.040865602, -0.040855378, -0.040843692, -0.040833467, -0.038420415, -0.034805218, -0.033591388, -0.029979112, -0.026365376, -0.020348812, -0.011933802, -0.004718015, 0.006098363, 0.0169162, 0.034939601, 0.051759396, 0.069781335, 0.09260747, 0.113029316, 0.132251939, 0.156278757, 0.175502841, 0.197128292, 0.210346036, 0.222366019, 0.233194081, 0.235607134, 0.233218913, 0.232028454, 0.226036722, 0.222447817, 0.217655307, 0.215265626, 0.214078088, 0.217693285, 0.222509166, 0.232127781, 0.244147763, 0.26097048, 0.280196024, 0.298219424, 0.31504068, 0.329459108, 0.347482508, 0.358300346, 0.369116723, 0.377534655, 0.383551219, 0.383562904, 0.381171762, 0.372780123, 0.35958575, 0.345189232, 0.32238939, 0.27918669, 0.255183243, 0.214377529, 0.183172901, 0.149569828, 0.113566848, 0.088365639, 0.06676648, 0.049967135, 0.031970027, 0.019974876, 0.009178949, -0.000414835, -0.008806474, -0.018397336, -0.023188385, -0.026777289, -0.030369115, -0.030369115, -0.030369115, -0.030369115>;

"pts:=[[[[[0],[3.016315859],[5.753536424],[8.489713643],[11.22797756],[13.96519812],[16.09414745],[18.83136802],[21.26445296],[24.00167353],[&vellip;]]]],["99 element Vector[column]"]],[[[[[-0.041150436],[-0.041150436],[-0.04113729],[-0.038722776],[-0.041110998],[-0.041097851],[-0.041087627],[-0.04107448],[-0.041062795],[-0.041049649],[&vellip;]]]],["99 element Vector[column]"]]"

> 

A := unapply(Spline(pts[1], pts[2], x, degree = 3), x):
> 

plot(A(x), x=0..250);

> 

PDE := diff(C(x,t),t) = diff(C(x,t),x,x) - diff(C(x,t),x) + 5*diff(C(x,t), x$4);

diff(C(x, t), t) = diff(diff(C(x, t), x), x)-(diff(C(x, t), x))+5*(diff(diff(diff(diff(C(x, t), x), x), x), x))

> 

IBC := C(x,0)=A(x), D[1](C)(0,t)=0, D[1](C)(250,t)=0, D[1,1](C)(0,t)=0, D[1,1](C)(250,t)=0:
> 

pdsol := pdsolve(PDE, {IBC}, numeric);

_m140347094656320

> 

pdsol:-plot3d(x=0..250, t=0..600, orientation=[120,65,0]);

In the long run the solution converges to a constant.
Here is what it looks like at t=600.

> 

pdsol:-plot(t=600, view=-0.05..0, gridlines=false);


 

Download mw.mw

Another way...

Rouben Rostamian 7711
May 30 2019
1 0

As an alternative to what tomleslie has suggested, you may also want to try:

restart;
with(plots):
P, Q := <1,1>, <2,2>:
frame := t -> pointplot([P, (1-t)*P + t*Q], connect, color="Green"):
nf := 50:
display([seq(frame(i/(nf-1)), i=0..nf-1)], insequence);

 

Solution via finite differences...

Rouben Rostamian 7711
May 24 2019
0 2

As far as I can tell, Maple's numerical solver is not equipped to solve a heat conduction problem with discontinuous conductivity (but I will be happy to be corrected on this.)  For now, I have written a finite differences solver with a specific goal of handling such cases which you will find in
https://www.mapleprimes.com/posts/210623-A-Finite-Difference-Scheme-For-The-Heat

Get the worksheet from there and then solve your problem by entering:
k := x -> piecewise(x < 10, 10, 20);
sol := heat_solve(c=k, a=0, b=20, alpha=20, beta=10, u__0=0, T=15, n=50, m=50):

Then animate the solution:
display(convert(sol, list), insequence);
 

How about this?...

Rouben Rostamian 7711
May 20 2019
2 3
> 

restart;

The equation of heat condition in a medium of
variable conductivity k(x) is
"(&PartialD; u)/(&PartialD; t)=(&PartialD;)/(&PartialD; x)(k(x) (&PartialD; u)/(&PartialD; x))."

In your case the conductivity is a piecewise constant
function which jumps from k__1 to k__2 at the interval's

midpoint. In principle we should be able to take

"k(x) = `k__1`+(`k__2`-`k__1`) Heaviside(x-L),"

or equivalently

"k(x)=`k__1`+(`k__2`-`k__1`) piecewise(x<L, `k__1`, `k__2`)."

I am unable to make either of these representations
work with Maple's numeric solver, so I do the next
best thing which is to approximate the Heaviside function
with a smooth function such as

> 

myHeaviside := x -> (1+erf(4*x))/2;

proc (x) options operator, arrow; 1/2+(1/2)*erf(4*x) end proc

Let's compare myHeaviside with the real thing:

> 

plot([Heaviside(x), myHeaviside(x)], x=-10..10,
  color=[red,blue], axes=frame, gridlines=false,
  legend=[Heavisie, myHeaviside]);

OK, then, let's start:

> 

pde := Diff(u(x,t),t) = Diff(k(x)*Diff(u(x,t),x), x);

Diff(u(x, t), t) = Diff(k(x)*(Diff(u(x, t), x)), x)

> 

ibc := u(x,0)=0, u(0,t)=v1,  u(2*L,t)=v2;

u(x, 0) = 0, u(0, t) = v1, u(2*L, t) = v2

> 

k := x -> k1 + (k2-k1)*myHeaviside(x-L);

proc (x) options operator, arrow; k1+(k2-k1)*myHeaviside(x-L) end proc

> 

L := 10: v1 := 20: v2 := 10: k1 := 10: k2 := 20:

To handle the sharp change in conductivity, we refine the spacestep size from

its default value of 2*L/20 by dividing it by 5.

> 

pdsol := pdsolve(pde, {ibc}, numeric, spacestep=(2*L/20)/5);

_m139639715851360

> 

pdsol:-animate(t=15, frames=50, labels=[x,u(x,t)],
  title="time = %f");

Download mw.mw

PS:  Looking at the animation closely, we see that the graph wobbles a bit near the initial time. That's an artifact due to our insufficiently small step size. Reducing the spacestep from (2*L)/20/5 to (2*L)/20/20 makes the wobble go away.

Geodesics...

Rouben Rostamian 7711
May 08 2019
5 5

Note: Figures in this worksheet involve semi-transparent surfaces.  These do not display properly on the web.  You will need to look in the worksheet itself to see them.  For reference, this stand-alone animation shows a representative case.

Geodesics on a tetrahedron

This worksheet provides tools for computing and

plotting geodesics on a tetrahedron.  The tetrahedron

need not be regular, although we take a regular tetrahedron

for the illustrations here.

> 

restart;
> 

with(plots):

Vertices of a regular tetrahedron

> 

pi[1], pi[2], pi[3], pi[4] := <0,0,0>, <1,0,0>, <1/2,sqrt(3)/2,0>, <1/2, sqrt(3)/6, sqrt(6)/3>:

Plot the tetrahedron with semi-transparent faces.

> 

tetra := plottools:-tetrahedron(convert~([pi[1],pi[2],pi[3],pi[4]], list)):
T := display([
  display(tetra, style=line, color=black),
  display(tetra,color="PeachPuff", transparency=0.2, style=surface)
], scaling=constrained):

A face of the tetrahedron is specified through the indices u, v, w of its vertices.

E.g., u, v, w = 1, 4, 3corresponds to the face whose vertices are `&pi;__1`, `&pi;__4`, `&pi;__3`.

Along the edge connecting the vertices u and w we take a point P given by
P = a*`&pi;__u`+(1-a)*`&pi;__w`.

Along the edge connecting the vertices v and w we take a point Q given by

Q = b*`&pi;__v`+(1-b)*`&pi;__w`.

If we move from P to Q along the directed line segment PQ, and then

continue on the surface of the tetrahedron along a geodesic path, we

leave the previous face and enter a new face with vertices u__new, v__new, `w__new `

within which the geodesic traces a line segment QR.  The following proc

calculates and returns the indices u__new, v__new, w__new of the new vertices,

and the values of a and bcorresponding to the points Q and R on the

new face.

Later on, we construct a geodesic on the tetrahedron's surface through

repeated applications of this proc.

> 

doit := proc(u, v, w, a, b)
  local u_new, v_new, w_new, a_new, b_new,
        fourth_vertex := ({1,2,3,4} minus {u,v,w})[],
        U := pi[u],
        V := pi[v],
        W := pi[w],
        P := a*U + (1-a)*W,
        Q := b*V + (1-b)*W,
        T := V+W-U,
        s1 := b*a/evalf((a-b)),
        t1 := -(b-1)*a/evalf(b);
  if t1 <= 1 then
    w_new := v;
    u_new := w;
    v_new := fourth_vertex;
    a_new := 1 - b;
    b_new := t1;
  else
    w_new := w;
    u_new := v;
    v_new := fourth_vertex;
    a_new := b;
    b_new := s1;
  end if;
  return u_new, v_new, w_new, a_new, b_new, Q;
end proc:
> 

 

Examples

 

In all of the following examples we begin with

u, v, w = 1, 2, 3, which means that the first line

segment of the geodesic lies in the `&pi;__1`, `&pi;__2`, `&pi;__3` plane

(which is the tetrahedron's base in our case).

 

Consider the vectors A = `&pi;__2`-`&pi;__1` and B = `&pi;__3`-`&pi;__1` as a basis

of the vectors in the plane of the base of the tetrahedron..

Let C = A*m+B*n.  For the first segment of the geodesic

we take a line within the base which parallels the vector C.

If the ratio m/n is rational, then the geodesic will be a closed

path, otherwise the geodesic curve will be dense in the

tetrahedron's surface.

 

The parameter a that appears in the computations below

may be assigned a value between 0 and 1. It specifies a

parallel displacement of the initial line segment

(which is parallel to the vector C as noted above.)  If a = 1,

the segment goes through the vertex `&pi;__1`.  If "a=0," it goes
through the vertex `&pi;__3`.  Don't change the parameter "b."

> 

 

Example 1:  m=1, n=0
 

We construct and plot 100 segments of the geodesic

but the result is a quadrilateral.  The plot is folded

over itself 25 times.

> 

u:=1: v:=2: w:=3: m:=1: n:=0: a:=0.7: b := m/(m+n)*a:
pts := a*pi[u] + (1-a)*pi[w]:  # starting point
for i from 1 to 100 do
  (u,v,w,a,b,q) := doit(u,v,w,a,b);
  pts := pts, q;
end do:
display([T, pointplot3d([pts], connect, color=red)],
  axes=none, orientation=[160,75,0]);

 
> 

 

> 

 

Example 2:  m=1, n=1
 
> 

u:=1: v:=2: w:=3: m:=1: n:=1: a:=0.7: b := m/(m+n)*a:
pts := a*pi[u] + (1-a)*pi[w]:  # starting point
for i from 1 to 100 do
  (u,v,w,a,b,q) := doit(u,v,w,a,b);
  pts := pts, q;
end do:
display([T, pointplot3d([pts], connect, color=red)],
  axes=none, orientation=[160,75,0]);

 
> 

 

> 

 

Example 3:  m=2, n=3
 
> 

u:=1: v:=2: w:=3: m:=2: n:=3: a:=0.7: b := m/(m+n)*a:
pts := a*pi[u] + (1-a)*pi[w]:  # starting point
for i from 1 to 100 do
  (u,v,w,a,b,q) := doit(u,v,w,a,b);
  pts := pts, q;
end do:
display([T, pointplot3d([pts], connect, color=red)],
  axes=none, orientation=[160,75,0]);

 

Example 4:  m=sqrt(2), n=1
 

The ratio m/n is irrational, therefore the geodesic is dense

on the surface of the torus.

> 

u:=1: v:=2: w:=3: m:=sqrt(2.0): n:=1: a:=0.7: b := m/(m+n)*a:
pts := a*pi[u] + (1-a)*pi[w]:  # starting point
for i from 1 to 100 do
  (u,v,w,a,b,q) := doit(u,v,w,a,b);
  pts := pts, q;
end do:
display([T, pointplot3d([pts], connect, color=red)],
  axes=none, orientation=[160,75,0]);

 
> 

 


 

Download geodesics.mw

 

Solution...

Rouben Rostamian 7711
May 03 2019
6 2

Any sequence a__0, a__1, a__2 .. ()generated according to

(*)         `a__i+2` < (1/4)*`a__i+1`^2/a__ifor all i = 0, 1, 2, () .. ()

has the desired property.  For instance, let's take

> 

restart;
> 

recur := a(0)=1, a(1)=1/5, a(i+2)=1/5*a(i+1)^2/a(i);

a(0) = 1, a(1) = 1/5, a(i+2) = (1/5)*a(i+1)^2/a(i)

> 

rsolve({recur}, {a(i)}):
ans := unapply(%, i);

proc (i) options operator, arrow; 5^(-(1/2)*i*(i+1)) end proc

This sequence converges to zero more slowly than that suggested by vv, enabling us to

extend the test further down in the sequence.  Applying vv's method of testing, we have

> 

P:= n -> add(ans(k)*x^k, k=0..n):
seq( [n, nops([fsolve(P(n))])], n=1..40);

[1, 1], [2, 2], [3, 3], [4, 4], [5, 5], [6, 6], [7, 7], [8, 8], [9, 9], [10, 10], [11, 11], [12, 12], [13, 13], [14, 14], [15, 15], [16, 16], [17, 17], [18, 18], [19, 19], [20, 20], [21, 21], [22, 22], [23, 23], [24, 24], [25, 25], [26, 26], [27, 27], [28, 28], [29, 29], [30, 30], [31, 31], [32, 32], [33, 33], [34, 34], [35, 35], [36, 36], [37, 37], [38, 38], [39, 39], [40, 40]

 

Explanation:

The inequality (*) comes from the paper:

 

David C. Kurtz, A Sufficient Condition for All the Roots of a Polynomial To Be Real
The American Mathematical Monthly, Vol. 99, No. 3 (Mar., 1992), pp. 259-263

 

where it is shown that (*) is sufficient to make all roots

of a polynomial real and distinct.

 

 

Download real-roots.mw

Brute-force calculation...

Rouben Rostamian 7711
April 27 2019
0 3

If n is not too large, then a brute-force approach (exhaustive search) will do. It appears that there are solutions provided that the value of a is small. Here is an illustration.

> 

restart;
> 

n := 10;  a := 1/30;

10

1/30

Generate a random vector with integer entries in the 0..99 range

> 

X := LinearAlgebra:-RandomVector[row](n, generator=rand(100));

Vector[row](10, {(1) = 89, (2) = 33, (3) = 17, (4) = 51, (5) = 55, (6) = 42, (7) = 82, (8) = 24, (9) = 89, (10) = 92})

> 

S := a*n*add(X);

574/3

> 

for k from 2 to n do
  for h from 1 to k-1 do
    if S/((n+1-k)*(n+k-2*h))
         <= min( X[k]-X[k-1], (X[h+1]-X[h])/(n+1-k) ) then
      printf("h=%d, k=%d\n", h, k);
    end if;
  end do;
end do;  

h=3, k=4
h=3, k=5
h=3, k=7
h=6, k=7
h=3, k=9
h=6, k=9
h=8, k=9

 

 

 

Download mw.mw

Coordplot...

Rouben Rostamian 7711
April 25 2019
2 0

plots:-coordplot(elliptic);

Animation...

Rouben Rostamian 7711
April 25 2019
3 3

Animation:

Download worksheet: conic-sections.mw

Alternatively...

Rouben Rostamian 7711
April 23 2019
2 0

Alternatively,

Eval(diff(f(r),r), r=C);

This may help...

Rouben Rostamian 7711
April 22 2019
4 0

I have Maple 2018 and I see the problem that you are facing. Here is how I would address the issue.

I don't know much about exportplot().  I use plotsetup() which seems to be better documented.

Put this in your .mpl file:

plotsetup(cps,
	plotoutput="test.eps",
	plotoptions="portrait, noborder, width=500, height=500");
plot3d(cos(x)*cos(y), x = -Pi .. Pi, y = -Pi .. Pi,
	axes=boxed, orientation=[55,75,0]);
plotsetup(default);

The resulting test.eps looks pretty much like what you get from exporting from a worksheet but it can benefit from some minor adjustments.

First, the bounding box is too tight—there are no margins around the figure.  To fix that we want to scale down the figure a little.  So we look into the test.eps file and find the line
0.1 0.1 scale

and change that to
0.09 0.09 scale

Now the size of the graph is good but it is not centered in the window.  To fix that, we translate the plot a little before scaling, as in
15 20 translate
0.09 0.09 scale
The result looks good to me. You may wish to play around with those adjustments to get something that looks good to you.

First 24 25 26 27 28 29 30 Last Page 26 of 53