Robert Israel

6542 Reputation

21 Badges

17 years, 166 days
University of British Columbia
Associate Professor Emeritus
North York, Ontario, Canada

MaplePrimes Activity

These are answers submitted by Robert Israel

You really have two different systems, one for gamma3(t) > c and the other for gamma3(t) < c.

I suggest you start with the system for gamma3(t) > c (since the initial condition has gamma3(0) > c), and integrate until gamma3(t) hits c (the events mechanism in dsolve provides a convenient way to do that).  Then switch to the system for gamma3(t) < c.

The matrix A is almost a rotation matrix.  I suspect you want to find the Euler angles.  There are various conventions for these, but one of them is

Matrix(3, 3, {(1, 1) = cos(alpha), (1, 2) = sin(alpha), (2, 1) = -sin(alpha), (2, 2) = cos(alpha), (3, 3) = 1}) Matrix(3, 3, {(1, 1) = cos(beta), (1, 3) = -sin(beta), (2, 2) = 1, (3, 1) = sin(beta), (3, 3) = cos(beta)}) Matrix(3, 3, {(1, 1) = 1, (2, 2) = cos(gamma), (2, 3) = sin(gamma), (3, 2) = -sin(gamma), (3, 3) = cos(gamma)})

You won't get that exactly, but we can try to get a good fit using LSSolve in the Optimization package.  There's a slight complication: gamma is not a good name for a variable in LSSolve, because gamma is actually the name of a constant in Maple.

After entering the matrix above as B and your matrix as A:

> Optimization[LSSolve](subs(gamma=c,   [seq(seq(A[i,j]-B[i,j],i=1..3),j=1..3)]),alpha=0..2*Pi,beta=0..2*Pi,c=0..2*Pi);

The first number being so small indicates an extremely good fit.  Here's the matrix with these values of the angles.

> subs(gamma=c, %[2], B);

1) Use the parametric form of plot, and ArcLength in the Student[Calculus1] package.

2) Use polarplot in the plots package.  I think you'll have to set up the area integral yourself: use int for integration.

3) Use implicitplot3d in the plots package.  Again you'll want to set up a multiple integral yourself.

Put a semicolon ; after the third ::nonnegint or just remove it altogether (you don't have to specify a return type for a procedure)

Take out that local (you don't have any local variables). 
Also, you have one Ackermann( too many in the else clause


if m=0 then n+1
elif n=0 then Ackermann(m-1,1)
else Ackermann(m-1,Ackermann(m,n-1))
end if
end proc;

The sum command, like most Maple commands, starts by evaluating its arguments.
So F(0,i) is called with i as a symbolic variable, and the result is 0 (this is in fact correct for positive integers i).  Unfortunately, you want it at i=0, so you get an incorrect result.  Maple's symbolic sum and product commands, when given a range for the index such as a .. b where b is symbolic, may implicitly make the assumption that b > a in order to get a closed-form formula.   A cure for this would be to use add instead of sum.  This has special evaluation rules so it will only evaluate its first argument with numerical values substituted for the index variable.  And by the way, you shouldn't use the global variable i for the index variable in F.  Thus:

> F:= proc(a,b) local i; simplify(product(2^(a+1-i)-1, i=1..b)) end proc;
    add(F(0,i), i = 0 .. 0);


This is what functions are for.  Try:

> A:= (r,z) -> int( f(r,z, teta), teta=0..0.78);

Then, for example,


I think the singularity in your PDE at r=0 could be causing trouble.  You might try a change of variables to remove the singularity.  If you take r^2 = s:

> with(PDEtools):

I'm not sure what you'll do about a boundary condition at s=0, though: the change of variable has made D[1](Phi)(0,t) = 0 automatically true if B is differentiable.

The point is that 1 can appear in unexpected places in Maple's internal representations of things, so anything of the form subs(1 = ..., ...) is likely to give unexpected results.  In this case it seems that at some level a/b is represented as a^1 * b^(-1).  Reasonable enough, I guess.

Have you tried implementing the Gale-Shapley algorithm?  See the pseudocode at

> max(A,B) assuming A::RealRange(0,1), B::RealRange(2,3);


Maple can calculate the inverse of a Matrix A as A^(-1).  If you do want the adjoint, that is Adjoint(A) (in the LinearAlgebra package).  You shouldn't have to worry about pivot elements.  What you do have to worry about, however, is that the inverse of a 6 x 6 symbolic matrix will be horrendously complicated.  In general (if there is no special simplification) each entry will be a fraction whose denominator is Determinant(A), which is the sum of 6! = 720 terms, each being the product of 6 entries of A.  For example, the inverse of

has a length, as reported by the length command, of 2246858.  If the matrix entries are more complicated, the length will increase even more.  Further computations with such a monster are likely to be very difficult.  You should ask yourself whether you really need an explicit expression for A-1.  What do you intend to do with it?  Could you reformulate things so that an explicit inverse is not needed?

A simpler way is to embed your expression in a one-parameter family, solve for the parameter, and differentiate.  The result is a first-order differential equation that has the family as its general solution.  For example:

> family := y(x) = C/x + sin(x);
   sol:= solve(family, C);
   de1:= diff(sol, x) = 0;


Of course there are many ways to put in a parameter, corresponding to lots of differential equations. 

> F:= N -> Matrix(N,N, 
proc(i,j) if 2+j-i < 0 then 0 else (((1+j-i)*N+j)/(2+j-i)!) *z[2+j-i] end if end proc);

Then, for example,

> F(5);

Something like this, perhaps?

> for ro from 1 to 10 do
      n1:= 2*ro - 1:
      print('A'^n1 = A[1]^n1, 'A'^(n1+1) = A[1]^(n1+1))
    end do:

I'm not sure what you mean by "the equation which y(x) satisfies".  One equation that y(x) satisfies is y(x) = 1/x + sin(x).  You could transform this in a number of ways.  Or do you want a differential equation?  Perhaps something like


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