Rouben Rostamian

MaplePrimes Activity

These are answers submitted by Rouben Rostamian

T := t -> (t+459.47) * 5/9;
seq(printf("%10d %10.2f\n", t, T(t)), t=100..110);
       100     310.82
       101     311.37
       102     311.93
       103     312.48
       104     313.04
       105     313.59
       106     314.15
       107     314.71
       108     315.26
       109     315.82
       110     316.37


Inversion in a sphere is much simpler than what you have in your worksheet.  Here is how I would do it.

Inversion in a sphere

Consider a sphere Sof radius R centered at the origin.  Let `#mover(mi("r"),mo("→"))` be the position vector of any

point Pin space. The inversion of P with respect to S is the point given by the position

vector `#mover(mi("ρ",fontstyle = "normal"),mo("→"))` = R^2*`#mover(mi("r"),mo("→"))`/abs(`#mover(mi("r"),mo("→"))`)^2  provided that `#mover(mi("r"),mo("→"))` is nonzero.  We see that inversion maps the outside

of the sphere into the inside, and vise versa.  Let's note that in Cartesian coordinates
`#mover(mi("r"),mo("&rarr;"))` ="<x,y,z>     implies    (rho)=(R^(2))/(x^(2)+y^(2)+z^(2))*<x,y,z>."


Here we illustrate this graphically by inverting a torus in a sphere.




The equation of a torus of radii a and b, displaced by d along the x axis in the Cartesian coordinates:

T := < d + (a + b*cos(s))*cos(t), (a + b*cos(s))*sin(t), b*sin(s) >;

Vector(3, {(1) = d+(a+b*cos(s))*cos(t), (2) = (a+b*cos(s))*sin(t), (3) = b*sin(s)})

The inversion of the torus in the sphere S

T__inverse := R^2/(T^+ . T) * T;

Vector(3, {(1) = R^2*(d+(a+b*cos(s))*cos(t))/((d+(a+b*cos(s))*cos(t))^2+(a+b*cos(s))^2*sin(t)^2+b^2*sin(s)^2), (2) = R^2*(a+b*cos(s))*sin(t)/((d+(a+b*cos(s))*cos(t))^2+(a+b*cos(s))^2*sin(t)^2+b^2*sin(s)^2), (3) = R^2*b*sin(s)/((d+(a+b*cos(s))*cos(t))^2+(a+b*cos(s))^2*sin(t)^2+b^2*sin(s)^2)})

Pick parameters for illustration:

params := { R=10, d=8, a=3, b=1 };

{R = 10, a = 3, b = 1, d = 8}

        sphere(eval(R, params), style=line),
        plot3d(eval(T, params), s=0..2*Pi, t=0..2*Pi, style=surface, color="Green"),
        plot3d(eval(T__inverse, params), s=0..2*Pi, t=0..2*Pi)
], scaling=constrained, viewpoint=[circleleft, frames=100], axes=none);





There are mnay ways of doing it, depending on what you want to do with the results.  Also it depends on how you represent quaternions.  You say that you want to represent a quaternion as a vector of length 4, but in your example your write qaa and qab as lists.  Lists and vectors are different types of objects in Maple.

Let's say you represent a quaternion as a list of length 4.  Then this will do what you want:


mult := proc(p::list,q::list)
                p[1]*q[1] - p[2]*q[2] - p[3]*q[3] - p[4]*q[4],
                p[1]*q[2] + p[2]*q[1] + p[3]*q[4] - p[4]*q[3],
                p[1]*q[3] - p[2]*q[4] + p[3]*q[1] + p[4]*q[2],
                p[1]*q[4] + p[2]*q[3] - p[3]*q[2] + p[4]*q[1]
end proc:



[a*e-b*f-c*g-d*h, a*f+b*e+c*h-d*g, a*g-b*h+c*e+d*f, a*h+b*g-c*f+d*e]

mult([8,4,6,2], [2,4,6,8]);

[-52, 76, 36, 68]

That said, I much prefer the [scalar, vector] pair representation of quaternions (see the Wikipedia article) rather than lists of length 4.  But that again depends on what you want to do with them.


de := a*diff(y(x),x,x) + b*diff(y(x),x) + c*y(x) = 0;
dsolve(de) assuming b^2 < 4*a*c;



P := [[1,1], [6,6], [7,12], [8,15], [10,17]];

[[1, 1], [6, 6], [7, 12], [8, 15], [10, 17]]


Returns the endpoint coordinates of the segment of the "second gradient".

select_second := proc(P::list)
        local i, s, Q;

        if nops(P) < 3 then
                error "the argument P should have at least 3 points";
        end if;

        # calculate the sequence of slopes
        seq(P[i+1]-P[i], i=1..nops(P)-1);
        s := map(p -> p[2]/p[1], [%]);

        # find the first index where consecutive slops are different
        for i from 2 to nops(P)-1 do
                if s[i] <> s[i-1] then
                end if;
        end do;
        if i = nops(P) then
                error "there is no second gradient in the given P";
                Q := [ P[i], P[i+1] ];
        end if;

        return Q;
end proc:


[[6, 6], [7, 12]]




Okay, your hand-written calculations are good.  Here I will do those in Maple.


In your hand calculations the symbol m is used with two

different meanings, first as an unknown exponent, and

second, as an abbreviation for "meter".  That's not a good.

I will write k for the exponent to remove the confusion.

expr := v = P/L * eta^n * R^k;

v = P*eta^n*R^k/L

Below, we work with expr/v, in order to make the left-hand side

of the expression into the dimensionless "1".

subs(v=m/s, P=N/m^2, L=m, eta=N*s/m^2, R=m, expr/v);
eq := simplify(%) assuming m>0, N>0, s>0;

1 = s*N*(N*s/m^2)^n*m^k/m^4

1 = m^(-4-2*n+k)*s^(1+n)*N^(1+n)

We want this equation to hold identically for all m, s, N.  Therefore the

exponents of these quantities should be zero, that is,

sys := { -4 - 2*n + k = 0, 1 + n = 0, 1 + n = 0 };

{1+n = 0, -4-2*n+k = 0}


{k = 2, n = -1}

This answer agrees with your hand calculations.


Note:  I formed the expression "sys" above by typing the three

exponents by hand.  It is possible to have Maple extract those

exponents programmatically but that would go far beyond the

scope of the beginner user of Maple.




How do you like this?


# macro emulates Matlab's "jet" colorscheme:
macro(jet=["zgradient", ["DarkBlue", "Blue", "DeepSkyBlue", "Cyan", "LawnGreen", "Yellow", "Orange", "Red", "DarkRed"]]):

plot3d((1-xi^2)*(1+tau), tau=-1..1, xi=-1..1,
    colorscheme=jet, orientation=[-120,72,0], axes=framed);

plot3d((1-xi^2)*(1+tau), tau=-1..1, xi=-1..1,
    colorscheme=jet, gridstyle=triangular, orientation=[-120,72,0]);

As Kitonum has noted, this is not an infinite fraction.  The integrand is z = ..  Then we see that z^2 = x*z, and therefore z = x.  Then the integral of z is x^2/2 + constant.  QED.

Forget about solve.  After you define psi, just do:

plots:-contourplot(psi, x=-1..1, eta=-2..2);


Try this:

plot(x^2, x=-1..1, labels=[typeset(NFE), typeset(log__10(`Max Err`))]);

Here is a slightly modified version of your worksheet.  The determinant grows too fast for efficient numerical calculations.  Therefore I have scaled it by dividing by cosh(...) which does not affect the roots since cosh is never zero.  Then you see that you can calculate the first five roots without any fuss.



If you look closely at your equations, you will see that at some places you have x, and at some other places you have x(t).  They should all be x(t).  Same goes for y(t) and u(t).

I have fixed these in the attached worksheet.  Note the technique—to define a differential equation dx/dt = F(x,y,u), we begin with defining the function (not an expression!) F, and then express the equation as dx(t)/dt = F(x(t), y(t), u(t)).

But note that the solution to your system is not defined in the entire length of the time interval that you have requested.  Since you have x, y, and u in the denominators, the solution stops when any of these variables reaches zero.



p := 2*u*x*y^2;
q := (x-1)/x;



doth := simplify(-(sqrt(6)*x^3*(y-1)^2+u*((3+sqrt(6))*x^3*(y-1)^2+3*y^2-9*x*y^2+3*x^2*(y^2+2*y-1)))/p);


ddotphi := simplify((y-1)*(sqrt(6)*x*y^2*(x-1)+u*(3*y^2-(6+sqrt(6))*x*y^2+x^2*((3+sqrt(6))*y^2+6*y-3)))/p);


x*(1-x)*doth - x^2*ddotphi/sqrt(6):
F1 := unapply(%, [x,y,u]);    # the right-hand side of ode1

proc (x, y, u) options operator, arrow; (1/2)*(1-x)*(-x^3*(y-1)^2*(u+1)*6^(1/2)-3*(x-1)*((y-1)^2*x^2+2*x*y^2-y^2)*u)/(u*y^2)-(1/12)*x*(y-1)*(x*y^2*(x-1)*(u+1)*6^(1/2)+3*((x-1)^2*y^2+2*x^2*y-x^2)*u)*6^(1/2)/(u*y^2) end proc

ode1 := diff(x(t),t) = F1(x(t),y(t),u(t));

diff(x(t), t) = (1/2)*(1-x(t))*(-x(t)^3*(y(t)-1)^2*(u(t)+1)*6^(1/2)-3*(x(t)-1)*((y(t)-1)^2*x(t)^2+2*x(t)*y(t)^2-y(t)^2)*u(t))/(u(t)*y(t)^2)-(1/12)*x(t)*(y(t)-1)*(x(t)*y(t)^2*(x(t)-1)*(u(t)+1)*6^(1/2)+3*((x(t)-1)^2*y(t)^2+2*x(t)^2*y(t)-x(t)^2)*u(t))*6^(1/2)/(u(t)*y(t)^2)

F2 := unapply(%, [x,y,u]);        # the right-hand side of ode2

proc (x, y, u) options operator, arrow; (1/2)*(y-1)*(-x^3*(y-1)^2*(u+1)*6^(1/2)-3*(x-1)*((y-1)^2*x^2+2*x*y^2-y^2)*u)/(y*u*x) end proc

ode2 := diff(y(t), t) = F2(x(t),y(t),u(t));

diff(y(t), t) = (1/2)*(y(t)-1)*(-x(t)^3*(y(t)-1)^2*(u(t)+1)*6^(1/2)-3*(x(t)-1)*((y(t)-1)^2*x(t)^2+2*x(t)*y(t)^2-y(t)^2)*u(t))/(y(t)*u(t)*x(t))

F3 := unapply(%, [x,y,u]);          # the right-hand side of ode3

proc (x, y, u) options operator, arrow; (1/2)*6^(1/2)*u*(u+1)*(x-1)/x end proc

ode3 := diff(u(t), t) = F3(x(t),y(t),u(t));

diff(u(t), t) = (1/2)*6^(1/2)*u(t)*(u(t)+1)*(x(t)-1)/x(t)


{t, diff(u(t), t), diff(x(t), t), diff(y(t), t), u(t), x(t), y(t)}

DEplot3d([ode1, ode2, ode3], [x(t), y(t), u(t)], t = -2 .. 2, x = 0 .. 1, y = 0 .. 1, u = -1 .. 0, [[x(0) = .9, y(0) = .25, u(0) = -.95], [x(0) = .9, y(0) = .25, u(0) = -.5], [x(0) = .85, y(0) = .2, u(0) = -.15], [x(0) = .25, y(0) = .5, u(0) = -.5], [x(0) = .7, y(0) = .5, u(0) = -.6]], linecolor = [red, green, black, navy, maroon], thickness = 2, axes = boxed, labels = [x(t), y(t), u(t)], stepsize = 0.1e-1, orientation = [30, 65]);

Warning, plot may be incomplete, the following errors(s) were issued:
   cannot evaluate the solution further right of .16790195, probably a singularity
Warning, plot may be incomplete, the following errors(s) were issued:
   cannot evaluate the solution further right of .51702724e-1, probably a singularity
Warning, plot may be incomplete, the following errors(s) were issued:
   cannot evaluate the solution further right of .57085115e-2, probably a singularity
Warning, plot may be incomplete, the following errors(s) were issued:
   cannot evaluate the solution further right of .54099962, probably a singularity




The English term for "exinscrit circle" in the second sentence is  "excircle".

Here is an animation of the excircle tangent to the side AC.  We see that it goes through the ellipse's vertex, as asserted.  I haven't attempted a proof but let me add an extra observation: It appears that the locus of the excircle's center is a vertical line.



Define the ellipse through its semi-major axis length a and focal length "c." 

a, c := 3, 2;

3, 2

The eccentricity

e := c/a;


The foci:


Proc to draw one animation frame

frame := proc(t)
        local r, A, T, E, obj, circ;
        r := a*(1-e^2)/(1 - e*cos(t));  # equation of ellipse in polar coords
        point(A, r*cos(t), r*sin(t));
        triangle(T, [B,C,A]);
        ellipse(E, [foci=[B,C], distance=2*a]);
        excircle(obj, T);
        circ := %[1];
        draw([T,E,circ], color=[red,"Green",blue], printtext=true, axes=none);
end proc:


nframes := 30:
plots:-display([seq(frame(i/nframes*Pi/2.0), i=1..nframes)], insequence, scaling=constrained);








I don't like the idea of using both psi and psi[2] in the same Maple worksheet.  In this particular case that's harmless but often it is not.  So I have changed psi[2] to w.


Do := f -> diff(f, r $ 2) + diff(f, phi $ 2)/r^2 - cos(phi)*diff(f, phi)/(r^2*sin(phi));

proc (f) options operator, arrow; diff(f, `$`(r, 2))+(diff(f, `$`(phi, 2)))/r^2-cos(phi)*(diff(f, phi))/(r^2*sin(phi)) end proc

pde := Do(Do(psi(t, r, phi)) - diff(psi(t, r, phi), t)/nu);

diff(diff(diff(diff(psi(t, r, phi), r), r), r), r)+(diff(diff(diff(diff(psi(t, r, phi), phi), phi), r), r))/r^2-4*(diff(diff(diff(psi(t, r, phi), phi), phi), r))/r^3+6*(diff(diff(psi(t, r, phi), phi), phi))/r^4-cos(phi)*(diff(diff(diff(psi(t, r, phi), phi), r), r))/(r^2*sin(phi))+4*cos(phi)*(diff(diff(psi(t, r, phi), phi), r))/(r^3*sin(phi))-6*cos(phi)*(diff(psi(t, r, phi), phi))/(r^4*sin(phi))-(diff(diff(diff(psi(t, r, phi), r), r), t))/nu+(diff(diff(diff(diff(psi(t, r, phi), phi), phi), r), r)+(diff(diff(diff(diff(psi(t, r, phi), phi), phi), phi), phi))/r^2+2*(diff(diff(psi(t, r, phi), phi), phi))/r^2-cos(phi)*(diff(diff(diff(psi(t, r, phi), phi), phi), phi))/(r^2*sin(phi))+2*cos(phi)^2*(diff(diff(psi(t, r, phi), phi), phi))/(r^2*sin(phi)^2)-2*cos(phi)*(diff(psi(t, r, phi), phi))/(r^2*sin(phi))-2*cos(phi)^3*(diff(psi(t, r, phi), phi))/(r^2*sin(phi)^3)-(diff(diff(diff(psi(t, r, phi), phi), phi), t))/nu)/r^2-cos(phi)*(diff(diff(diff(psi(t, r, phi), phi), r), r)+(diff(diff(diff(psi(t, r, phi), phi), phi), phi))/r^2+(diff(psi(t, r, phi), phi))/r^2-cos(phi)*(diff(diff(psi(t, r, phi), phi), phi))/(r^2*sin(phi))+cos(phi)^2*(diff(psi(t, r, phi), phi))/(r^2*sin(phi)^2)-(diff(diff(psi(t, r, phi), phi), t))/nu)/(r^2*sin(phi))

Substitute a special form of psi into the PDE.  This introduces a new unknown w:

exp(-_c[1]*t)*sin(phi)*r*(_C2*BesselJ(1,sqrt(_c[1]/nu)*sin(phi)*r) + _C3*BesselY(1,sqrt(_c[1]/nu)*sin(phi)*r)):
psi(t,r,phi) = %*w(r*sin(phi));
expr_tmp := simplify(eval(pde, %));

psi(t, r, phi) = exp(-_c[1]*t)*sin(phi)*r*(_C2*BesselJ(1, (_c[1]/nu)^(1/2)*sin(phi)*r)+_C3*BesselY(1, (_c[1]/nu)^(1/2)*sin(phi)*r))*w(r*sin(phi))

2*exp(-_c[1]*t)*((5/2)*_C2*((-(1/5)*cos(phi)^2*r^2*_c[1]+(1/5)*r^2*_c[1]-(3/5)*nu)*(D(w))(r*sin(phi))+((2/5)*nu*r*(cos(phi)-1)*(cos(phi)+1)*((D@@3)(w))(r*sin(phi))+((cos(phi)^2*r^2*_c[1]-r^2*_c[1]+(3/5)*nu)*((D@@2)(w))(r*sin(phi))-(1/5)*((D@@4)(w))(r*sin(phi))*nu*r^2*(cos(phi)-1)*(cos(phi)+1))*sin(phi))*r)*BesselJ(1, (_c[1]/nu)^(1/2)*sin(phi)*r)+(5/2)*_C3*((-(1/5)*cos(phi)^2*r^2*_c[1]+(1/5)*r^2*_c[1]-(3/5)*nu)*(D(w))(r*sin(phi))+((2/5)*nu*r*(cos(phi)-1)*(cos(phi)+1)*((D@@3)(w))(r*sin(phi))+((cos(phi)^2*r^2*_c[1]-r^2*_c[1]+(3/5)*nu)*((D@@2)(w))(r*sin(phi))-(1/5)*((D@@4)(w))(r*sin(phi))*nu*r^2*(cos(phi)-1)*(cos(phi)+1))*sin(phi))*r)*BesselY(1, (_c[1]/nu)^(1/2)*sin(phi)*r)+r^3*(_c[1]/nu)^(1/2)*sin(phi)*(cos(phi)-1)*(cos(phi)+1)*(BesselJ(0, (_c[1]/nu)^(1/2)*sin(phi)*r)*_C2+BesselY(0, (_c[1]/nu)^(1/2)*sin(phi)*r)*_C3)*((D(w))(r*sin(phi))*_c[1]-2*((D@@3)(w))(r*sin(phi))*nu))/(r^2*nu*sin(phi)^2)

Let's get rid of the irrelevant multiplicative factor:

expr := r^2*nu*sin(phi)^2 / exp(-_c[1]*t) * expr_tmp;

5*_C2*((-(1/5)*cos(phi)^2*r^2*_c[1]+(1/5)*r^2*_c[1]-(3/5)*nu)*(D(w))(r*sin(phi))+((2/5)*nu*r*(cos(phi)-1)*(cos(phi)+1)*((D@@3)(w))(r*sin(phi))+((cos(phi)^2*r^2*_c[1]-r^2*_c[1]+(3/5)*nu)*((D@@2)(w))(r*sin(phi))-(1/5)*((D@@4)(w))(r*sin(phi))*nu*r^2*(cos(phi)-1)*(cos(phi)+1))*sin(phi))*r)*BesselJ(1, (_c[1]/nu)^(1/2)*sin(phi)*r)+5*_C3*((-(1/5)*cos(phi)^2*r^2*_c[1]+(1/5)*r^2*_c[1]-(3/5)*nu)*(D(w))(r*sin(phi))+((2/5)*nu*r*(cos(phi)-1)*(cos(phi)+1)*((D@@3)(w))(r*sin(phi))+((cos(phi)^2*r^2*_c[1]-r^2*_c[1]+(3/5)*nu)*((D@@2)(w))(r*sin(phi))-(1/5)*((D@@4)(w))(r*sin(phi))*nu*r^2*(cos(phi)-1)*(cos(phi)+1))*sin(phi))*r)*BesselY(1, (_c[1]/nu)^(1/2)*sin(phi)*r)+2*r^3*(_c[1]/nu)^(1/2)*sin(phi)*(cos(phi)-1)*(cos(phi)+1)*(BesselJ(0, (_c[1]/nu)^(1/2)*sin(phi)*r)*_C2+BesselY(0, (_c[1]/nu)^(1/2)*sin(phi)*r)*_C3)*((D(w))(r*sin(phi))*_c[1]-2*((D@@3)(w))(r*sin(phi))*nu)

Simplify and then introduce eta = r*sin(phi):

simplify(expr, {cos(phi)^2 = 1 - sin(phi)^2}):
algsubs(r*sin(phi)=eta, %):
ode := collect(%, [(D@@4)(w)(eta), (D@@3)(w)(eta), (D@@2)(w)(eta), (D@@2)(w)(eta)], simplify);

eta^3*nu*(BesselJ(1, (_c[1]/nu)^(1/2)*eta)*_C2+BesselY(1, (_c[1]/nu)^(1/2)*eta)*_C3)*((D@@4)(w))(eta)-2*(BesselJ(1, (_c[1]/nu)^(1/2)*eta)*_C2+BesselY(1, (_c[1]/nu)^(1/2)*eta)*_C3-2*eta*(_c[1]/nu)^(1/2)*(BesselJ(0, (_c[1]/nu)^(1/2)*eta)*_C2+BesselY(0, (_c[1]/nu)^(1/2)*eta)*_C3))*eta^2*nu*((D@@3)(w))(eta)+3*eta*(-(5/3)*_c[1]*eta^2+nu)*(BesselJ(1, (_c[1]/nu)^(1/2)*eta)*_C2+BesselY(1, (_c[1]/nu)^(1/2)*eta)*_C3)*((D@@2)(w))(eta)-3*(_C2*(-(1/3)*_c[1]*eta^2+nu)*BesselJ(1, (_c[1]/nu)^(1/2)*eta)+_C3*(-(1/3)*_c[1]*eta^2+nu)*BesselY(1, (_c[1]/nu)^(1/2)*eta)+(2/3)*(_c[1]/nu)^(1/2)*eta^3*_c[1]*(BesselJ(0, (_c[1]/nu)^(1/2)*eta)*_C2+BesselY(0, (_c[1]/nu)^(1/2)*eta)*_C3))*(D(w))(eta)

That's a fourth order ODE in the unknown w(eta).  Maple is unable to obtain the general solution of

that ODE.  Perhaps solutions can be obtained if you specify the values of some of the constants
_c__1, _C__3, etc.




Look at this post and the ensuing discussion on ways to handle discontinuous coefficients in the heat equation.   Converting the heat equation into a system of two first order equations (suggested in that discussion by vv) works in your case.  Here it is.


pde_sys :=
        k(x)*diff(T(x,t),x) = Q(x,t),
        diff(Q(x,t),x) = c(x)*rho(x)*diff(T(x,t),t);

k(x)*(diff(T(x, t), x)) = Q(x, t), diff(Q(x, t), x) = c(x)*rho(x)*(diff(T(x, t), t))

bc := D[1](T)(0,t) = q0, T(1,t) = t0;

(D[1](T))(0, t) = q0, T(1, t) = t0

ic := T(x,0) = f0(x);

T(x, 0) = f0(x)

k        := x->piecewise(x<=0.5,3-1*x,x>0.5,2.5);
c        := x->piecewise(x<=0.5,1.2-0.4*x,x>0.5,1);
rho      := x->piecewise(x<=0.5,20-10*x,x>0.5,15);
f0       := x->t0*x;
t0       := 20;
q0       := -2;

k := proc (x) options operator, arrow; piecewise(x <= .5, 3-x, .5 < x, 2.5) end proc

c := proc (x) options operator, arrow; piecewise(x <= .5, 1.2+(-1)*.4*x, .5 < x, 1) end proc

proc (x) options operator, arrow; piecewise(x <= .5, 20-10*x, .5 < x, 15) end proc

proc (x) options operator, arrow; t0*x end proc




piecewise(x <= .5, 3-x, .5 < x, 2.5)*(diff(T(x, t), x)) = Q(x, t), diff(Q(x, t), x) = piecewise(x <= .5, 1.2-.4*x, .5 < x, 1)*piecewise(x <= .5, 20-10*x, .5 < x, 15)*(diff(T(x, t), t))

pdsol := pdsolve({pde_sys}, {ic, bc}, numeric, spacestep=0.01);


tmax := 20;


pdsol:-animate(T(x,t), t=0..tmax, frames=100);

pdsol:-animate(Q(x,t), t=0..tmax, frames=100);







4 5 6 7 8 9 10 Last Page 6 of 44