Advanced Search

Rouben Rostamian

8586 Reputation

23 Badges

19 years, 140 days
Contact Rouben Rostamian

MaplePrimes Activity


These are answers submitted by Rouben Rostamian

T*...

Rouben Rostamian 8586
January 15 2020
2 0

Here is how to construct T*.  I don't quite understand the definitions of the Ybar matrices. Perhaps you can explain with an example.

> 

restart;
> 

with(LinearAlgebra):
> 

N := 3;

3

> 

A := IdentityMatrix(N+1);

Matrix(4, 4, {(1, 1) = 1, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = 1, (2, 3) = 0, (2, 4) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 1, (3, 4) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 1})

> 

B := Matrix([t^k$k=0..N]);

Vector[row](4, {(1) = 1, (2) = t, (3) = t^2, (4) = t^3})

> 

interface(rtablesize=(N+1)^2):
> 

T_star := KroneckerProduct(A,B);

Matrix(4, 16, {(1, 1) = 1, (1, 2) = t, (1, 3) = t^2, (1, 4) = t^3, (1, 5) = 0, (1, 6) = 0, (1, 7) = 0, (1, 8) = 0, (1, 9) = 0, (1, 10) = 0, (1, 11) = 0, (1, 12) = 0, (1, 13) = 0, (1, 14) = 0, (1, 15) = 0, (1, 16) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 0, (2, 5) = 1, (2, 6) = t, (2, 7) = t^2, (2, 8) = t^3, (2, 9) = 0, (2, 10) = 0, (2, 11) = 0, (2, 12) = 0, (2, 13) = 0, (2, 14) = 0, (2, 15) = 0, (2, 16) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 0, (3, 4) = 0, (3, 5) = 0, (3, 6) = 0, (3, 7) = 0, (3, 8) = 0, (3, 9) = 1, (3, 10) = t, (3, 11) = t^2, (3, 12) = t^3, (3, 13) = 0, (3, 14) = 0, (3, 15) = 0, (3, 16) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 0, (4, 5) = 0, (4, 6) = 0, (4, 7) = 0, (4, 8) = 0, (4, 9) = 0, (4, 10) = 0, (4, 11) = 0, (4, 12) = 0, (4, 13) = 1, (4, 14) = t, (4, 15) = t^2, (4, 16) = t^3})

 

 

 

Download mw.mw

PS: The following interprets the notation in the defintion of Ybar in a way that leads to the solution that yiou are looking for.

> 

restart;
> 

N := 3:
> 

Y1 := <y[1,k]$k=0..N>;

Vector(4, {(1) = y[1, 0], (2) = y[1, 1], (3) = y[1, 2], (4) = y[1, 3]})

> 

Y2 := <y[2,k]$k=0..N>;

Vector(4, {(1) = y[2, 0], (2) = y[2, 1], (3) = y[2, 2], (4) = y[2, 3]})

> 

Y1_bar := <(y[2,k]*Y1)$k=0..N>;

_rtable[18446883978168624542]

> 

Y2_bar := <(y[1,k]*Y2)$k=0..N>;

_rtable[18446883978035869566]

 

Download mw2.mw

 

Cost function...

Rouben Rostamian 8586
January 15 2020
1 2

Suppose the const function is f(x) = x^2.  Then you would plot it with:
plot(x^2, x=0..4);
If you have some other cost function in mind, then you should tell us what it is.

 

Banded storage...

Rouben Rostamian 8586
January 15 2020
1 1

If the zero entries of B are not going to change in your calculations, then you should use banded storage for B.  That way, only the nonzero entries of B are stored, thus saving a ton of memory if your N is large.

> 

restart;
> 

N:=5;

5

> 

B := LinearAlgebra:-BandMatrix([[0],[$1..N]],0,N+1);

Matrix(6, 6, {(1, 1) = 0, (1, 2) = 1, (1, 3) = 0, (1, 4) = 0, (1, 5) = 0, (1, 6) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 2, (2, 4) = 0, (2, 5) = 0, (2, 6) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 0, (3, 4) = 3, (3, 5) = 0, (3, 6) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 0, (4, 5) = 4, (4, 6) = 0, (5, 1) = 0, (5, 2) = 0, (5, 3) = 0, (5, 4) = 0, (5, 5) = 0, (5, 6) = 5, (6, 1) = 0, (6, 2) = 0, (6, 3) = 0, (6, 4) = 0, (6, 5) = 0, (6, 6) = 0})

You can change the entries on the main diagonal or on the superdiagonal if needed:

> 

B[1,2] := 100;

100

But you cannot change the entries elsewhere because there is no memory set aside
for them:

> 

B[2,1] := 100;

Error, attempt to assign a value outside Matrix bands

> 

 
 

 

 

Solution...

Rouben Rostamian 8586
January 14 2020
0 1

Here is one way of doing it.  Make the field "undefined" outside of the desired region.

> 

restart;
> 

f := map2(piecewise,  x^2+y^2<=1, [y, -sin(x)-(1/10)*y], undefined);

f := [piecewise(x^2+y^2 <= 1, y, undefined), piecewise(x^2+y^2 <= 1, -sin(x)-(1/10)*y, undefined)]

> 

plots:-fieldplot(f, x = -1 .. 1, y = 0 .. 1, scaling=constrained);

 

Download mw.mw

 

No differentiation necessary...

Rouben Rostamian 8586
January 11 2020
1 3
 

You have:
u = exp(-lambda*z*sqrt(x^2+y^2)).

 

You want to get


v = z*exp(-lambda*z*sqrt(x^2+y^2))/(1+2*sqrt(x^2+y^2)*lambda).

 

No differentiation is necessary. You may produce v through

"v=(z u)/(1-(2 ln(u))/z)."
Here is how we verify this in Maple.

> 

restart;
> 

u = exp(-lambda*z*sqrt(x^2 + y^2));
v = z*u/(1 - 2*ln(u)/z);
simplify(%, {%%}) assuming positive;

u = exp(-lambda*z*(x^2+y^2)^(1/2))

v = z*u/(1-2*ln(u)/z)

v = z*exp(-lambda*z*(x^2+y^2)^(1/2))/(2*(x^2+y^2)^(1/2)*lambda+1)

 

Download mw.mw

 

Another way...

Rouben Rostamian 8586
January 05 2020
1 1

 [f(k) $k=0..5];  

Arbitrary length strings...

Rouben Rostamian 8586
January 01 2020
3 0

Some years ago I wrote this worksheet to make strings of arbitrary length from STL letters.  See if it is useful to you.

3D text

2015-05-31

2015-06-17 revised (see below)

2017-02-05 added a "thickening effect" to produce an "Egyptian" look.

The idea for this worksheet comes from a post by Acer in MaplePrimes:

    http://www.mapleprimes.com/questions/204335-Can-Rotate-3d-Text-Like-This-Be-Done-In-Maple

The goal is to produce text with 3D letters.  It requires the letter shapes to be available in the STL format.

I downloaded a set of such letter shapes in a file named All_Alphabet_Letters_A-Z.zip  from

     http://www.thingiverse.com/thing:15198

and saved it in the /usr/local/share/char-shapes-stl directory.

 

This worksheet defines two procs named string_to_plot3d() and string_to_plot3d_cylinder() which take a string as the first argument and a number called "gap" as the second argument.  These procs extract the 3D letters corresponding to the characters in the string, arrange them in the 3D space with "gap" distance between the characters, and return the corresponding PLOT3D structures.  The first proc arranges the characters in a linear fashion.  The second proc wraps the characters around a cylinder.

> 

restart;
> 

with(plottools): with(plots): with(StringTools):
> 

letters_dir := "/usr/local/share/char-shapes-stl/":

This table translates single character names to their corresponding STL file names.

> 

STL_shapes_table := table([
  seq(c = sprintf("Letter_%s.stl", c), c in "ABCDEFGHIJKLMNOPQRSTUVWXYZ"),
  "a" = "Letter_a_lower_case.stl",
  "b" = "Letter_b_lower_case.stl",
  "c" = "Letter_c_lower_case.stl",
  "1" = "one.stl",
  "2" = "two.stl",
  "3" = "three.stl",
  "4" = "four.stl",
  "5" = "five.stl",
  "6" = "six.stl",
  "7" = "seven.stl",
  "8" = "eight.stl",
  "9" = "nine.stl",
  "0" = "zero.stl",
  "@" = "Letter_at.stl",
  "&" = "and_sign.stl",
  #"(" = "bracket_open_close.stl",  # this gives a composite "()" shape, not
  #")" = "bracket_open_close.stl",  # very useful
  "$" = "dollar_sign.stl",
  "€" = "euro_sign.stl",
  "!" = "exclamation_mark.stl",
  "-" = "minus_sign.stl",
  "#" = "number_sign.stl",
  "+" = "plus_sign.stl",
  "£" = "pound_sign.stl",
  "?" = "question_mark.stl",
  "§" = "section_sign.stl",
  "¥" = "yen_sign.stl",
NULL]):

The Import() command, introduced in Maple2015, reads an STL shape and returns a PLOT3D structure:

> 

Import(cat(letters_dir, STL_shapes_table["@"]), title="");

 

This is a front-end to  Import().  It throws an error if the requested character is unavailable.

> 

getletter := proc(c::string)
  if not assigned(STL_shapes_table[c]) then
    error "shape %1 not available", c
  else
    Import(cat(letters_dir, STL_shapes_table[c]), title="");
  end if;
end proc:
> 

getletter("A");

 
> 

string_to_plot3d := proc(str, gap)
  local n, i, shape, a, b, p, x;
  n := length(str);
  p := NULL;
  x := 0;
  for i from 1 to n do
    if str[i] = " " then
      x := x + gap;
      next;
    end if;
    shape := getletter(str[i]);
    a := op(1, getdata(shape, 'rangesonly')[1]);
    b := op(2, getdata(shape, 'rangesonly')[1]);
    translate(shape, x-a, 0, 0);
    p := p, select(has, [op(%)], POLYGONS)[];
    x := x + b - a + gap;
  end do;
  return PLOT3D(p);
end proc:

Test 1:

> 

string_to_plot3d("MAPS", 5):
display(%, scaling=constrained, color="Chocolate",
    style=patchnogrid, axes=none, orientation=[-80,30,0]);

 

Test 2:

> 

string_to_plot3d("THIS IS A TEST", 10):
display(%, scaling=constrained, color="Chocolate",
   style=patchnogrid, axes=none, orientation=[-110,-30,0]);

 

Here we wrap the string around a cylinder.  Additionally, we apply a transform
that thickens the letters toward their bases, giving them a sort of "Egyptian" look.

> 

string_to_plot3d_cylinder := proc(str, gap)
  local n, i, shape, a, b, p, x, L, R;
  n := length(str);
  p := NULL;
  x := 0;
  for i from 1 to n do
    if str[i] = " " then
      x := x + gap;
      next;
    end if;
    shape := getletter(str[i]);
    a := op(1, getdata(shape, 'rangesonly')[1]);
    b := op(2, getdata(shape, 'rangesonly')[1]);
    translate(shape, x-a, 0, 0);
    p := p, select(has, [op(%)], POLYGONS)[];
    x := x + b - a + gap;
  end do;
  p := PLOT3D(p);
  # the "2.5" below is the thickening factor.  Letters are 2.5 times thicker at the base than at the top.
  p := transform((x,y,z)->[x,y,z*(1-(2.5-1)*y/40)])(p);
  L := x;
  R := L/(2*Pi);
  return transform((x,y,z) -> [(R+z)*cos(x/R), (R+z)*sin(x/R), y])(p);
end proc:

Test 3:

> 

string_to_plot3d_cylinder(" THIS IS A TEST ", 10):
display(%, color="Chocolate", scaling=constrained, axes=none, style=patchnogrid);

 
> 

 


 

Download 3d-letters.mw

Animation...

Rouben Rostamian 8586
December 30 2019
4 3
> 

restart;
> 

with(plots,animate):

Let the x axis be vertical and point up.  A rope of length L is attached to the x axis

at the point x = L.  The other end at at x = 0 is loose.

 

When the rope is set into motion, let y(x, t) be the horizontal displacement
at time tof the point whose vertical coordinate is "x."

 

Through applying Newton's law of motion we find that y(x, t) solves the
partial differential equation
    "(&PartialD;)/(&PartialD; x)(rho g x (&PartialD; y)/(&PartialD; x))=rho((&PartialD;)^2 y)/((&PartialD;)^( )t( )^2)"
where rho (not necessarily a constant) is the rope's linear density, and g is

the acceleration due to gravity.  This is derived under the assumption that

the displacement y(x, t)is small.

 

The details of the derivation and a solution through separation of

variables may be found at:

https://www.phys.uconn.edu/~rozman/Courses/m3410_18s/downloads/hanging-chain.pdf

Here I will go ahead and use that solution in order to produce
an animation in Maple.

The solution is actually expressed as an infinite sum but here I will take

only the first few terms for the sake of simplicity, like this:

> 

y := Sum(a[n]*BesselJ(0, 2*omega(n)*sqrt(x/g))*cos(omega(n)*t), n=1..4);

Sum(a[n]*BesselJ(0, 2*omega(n)*(x/g)^(1/2))*cos(omega(n)*t), n = 1 .. 4)

where

> 

omega := n -> 1/2*sqrt(g/L)*BesselJZeros(0,n);

proc (n) options operator, arrow; (1/2)*sqrt(g/L)*BesselJZeros(0, n) end proc

and a__n are arbitrary constant coefficients.  In practice the a__ns are determined

in terms of the initial conditions but I will not go into that.  Instead, I will just

pick arbitrary values for them.  To get realistic solutions, we need the a__n to be

reasonably small so that y is small as it is stipulated in the derivation of the

equation of motion.

Here are the parameter values.  Change them if you want.

> 

params := g=1, L=1;

g = 1, L = 1

Experiment: Let's take a__1 = .2, a__2 = .1 and the rest of the a__ns equal to zero.  Then

we get the following animation:

> 

eval(value(y), {a[1]=0.2, a[2]=0.1, a[3]=0, a[4]=0}):
%animate(plot, [[%, x, x=0..L], thickness=5, color="Indigo"],
         t=0..36.4, frames=100, scaling=constrained):
subs(params, %):
value(%);

 

Download hanging-rope.mw

 

Answer...

Rouben Rostamian 8586
December 27 2019
0 0

We know that the derivative of a constant is zero.  As far as Maple knows, u[i] does not depend on x and therefore diff(u[i],x)=0.  You need to tell Maple that u[i] is actually a function of x, like this:

n := 0;

for k from 0 to n do
#A[k] := sum((Diff(u[i](x), x))*(Diff(u[k-i](x), x)), i = 0 .. k);
A[k] := sum((diff(u[i](x), x))*(diff(u[k-i](x), x)), i = 0 .. k)
end do;

 

Animation...

Rouben Rostamian 8586
December 27 2019
3 0

Here's an animation to go with vv's illustration.

Worksheet: classic-puzzle.mw

In this variant the pieces move simultaneously:

Worksheet: classic-puzzle-alt.mw

A solution...

Rouben Rostamian 8586
December 14 2019
0 0

This will get you [−1,0]. There is no particularly useful way to get (−1,0) since Maple will strip away the parentheses.

a := solve({-1 < x, 0 < 2*x/(x^2 - 1), x < 1}, {x});
[(op~(a) minus {x})[]];

 

Reference...

Rouben Rostamian 8586
December 14 2019
2 1

Earl, have a look at the (lengthy!) discussion in 
https://www.mapleprimes.com/questions/125273-Dsolve-Events-How-To-Control-For-A-Sign-Change
regarding events.

In particular this response from Allan Wittkopf
https://www.mapleprimes.com/questions/125273-Dsolve-Events-How-To-Control-For-A-Sign-Change#comment125426
shows how to use the evenfired option.  I think he is the author of Maple's events handler, but I can be wrong.

 

See if you can figure this out...

Rouben Rostamian 8586
December 14 2019
1 1 
with(Student[Calculus1], SurfaceOfRevolution):
plots:-display([
  SurfaceOfRevolution(r, r=0..5, axis=vertical, output=plot, transparency=0.5),
  SurfaceOfRevolution(12-r^2, r=0..3, axis=vertical, output=plot)
], scaling=constrained, color=["Gold", "DarkRed"], caption="");

 

Will this do?...

Rouben Rostamian 8586
December 10 2019
0 2

restart

expr := Y = G__B*D1*G__A*G__f/(G__A*G__B*G__M*G__c+G__A*G__B*G__R+1)+G__B*G__A*G__c*Y__sp/(G__A*G__B*G__M*G__c+G__A*G__B*G__R+1)+G__B*D1*G__d/(G__A*G__B*G__M*G__c+G__A*G__B*G__R+1)

expand(expr/D1)

Y/D1 = G__B*G__A*G__f/(G__A*G__B*G__M*G__c+G__A*G__B*G__R+1)+G__B*G__A*G__c*Y__sp/(D1*(G__A*G__B*G__M*G__c+G__A*G__B*G__R+1))+G__B*G__d/(G__A*G__B*G__M*G__c+G__A*G__B*G__R+1)

``

 

Download mw.mw

 

Custom-made graphics...

Rouben Rostamian 8586
December 07 2019
3 2

This custom-made graphics imposes a color grid over the domain in order to visualize which parts get mapped where.

Worksheet: mw.mw


The original domain is on the left.  The mapped domain is on the right.

 

First 22 23 24 25 26 27 28 Last Page 24 of 58