## 6071 Reputation

15 years, 215 days

## Will this do?...

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

## The easy way...

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 of radius  centered at the origin.  Let  be the position vector of any

point in space. The inversion of  with respect to  is the point given by the position

vector   provided that  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
=

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

 > restart;
 > with(plots):
 > with(plottools):

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

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

The inversion of the torus in the sphere

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

Pick parameters for illustration:

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

 > plots:-display([         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);

## One possible solution...

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:

 > restart;
 > 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:
 >
 > mult([a,b,c,d],[e,f,g,h]);

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

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.

## This will get you sines and cosines...

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

## Solution...

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

 > plot(P);

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                         break;                 end if;         end do;         if i = nops(P) then                 error "there is no second gradient in the given P";         else                 Q := [ P[i], P[i+1] ];         end if;         return Q; end proc:
 > select_second(P);

## Calculations in Maple...

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

 > restart;

In your hand calculations the symbol  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  for the exponent to remove the confusion.

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

Below, we work with expr/, 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;

We want this equation to hold identically for all .  Therefore the

exponents of these quantities should be zero, that is,

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

 > solve(sys);

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.

## Emulating Matlab...

How do you like this?

```restart;

# 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]);
```

## Simple by hand...

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

Try this:

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

## Scale the determinant...

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.

## x versus x(t)...

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.

 > restart;
 > with(DEtools):
 > 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

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

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

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

 > sqrt(3/2)*u*(1+u)*q: F3 := unapply(%, [x,y,u]);          # the right-hand side of ode3

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

 > indets({ode1,ode2,ode3});

 > 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

## Animation...

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.

 > restart;
 > with(geometry):

Define the ellipse through its semi-major axis length  and focal length

 > a, c := 3, 2;

The eccentricity

 > e := c/a;

The foci:

 > point(B,0,0): point(C,2*c,0):

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:

Animate:

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

## Reduction to an ODE...

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.

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

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

Substitute a special form of  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, %));

Let's get rid of the irrelevant multiplicative factor:

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

Simplify and then introduce :

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

That's a fourth order ODE in the unknown .  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
, etc.

## Try this...

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.

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

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

 > ic := 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;

 > pde_sys;

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