Robert Israel

6577 Reputation

21 Badges

18 years, 209 days
University of British Columbia
Associate Professor Emeritus
North York, Ontario, Canada

MaplePrimes Activity


These are answers submitted by Robert Israel

As a general rule, never use fewer than the default 10 Digits for floating-point evaluation of a complicated expression.  Roundoff error will kill you.  If you want to round to fewer digits, do that afterwards.

> evalf(U . Sig . Vt);
map(fnormal, %, 9);
simplify(%, zero);

                         [1.     3.     -2.]
                         [                 ]
                         [2.     7.     5. ]
                         [                 ]
                         [-2.    -3.    4. ]
                         [                 ]
                         [5.     -3.    -2.]


The first question is, what equation do you really want to solve?  What you wrote is a vector, not an equation.  I suspect that what you really want is your vector to be equal to a given vector, let's say b.  Maple doesn't solve vector equations directly, but you can convert to lists or sets and solve for the components.

 

> with(VectorCalculus):
W:= Vector(3,symbol=w):
V:= Vector(3,symbol=v);
B:= Vector(3, symbol=b):
eqs:= convert((W &x V)/((W &x V).(W &x V)) - B, list);

[(w[2] v[3] - w[3] v[2])//                       2
[                        \(w[2] v[3] - w[3] v[2])

                            2                          2\         
   + (w[3] v[1] - w[1] v[3])  + (w[1] v[2] - w[2] v[1]) / - b[1],

  (w[3] v[1] - w[1] v[3])//                       2
                          \(w[2] v[3] - w[3] v[2])

                            2                          2\         
   + (w[3] v[1] - w[1] v[3])  + (w[1] v[2] - w[2] v[1]) / - b[2],

  (w[1] v[2] - w[2] v[1])//                       2
                          \(w[2] v[3] - w[3] v[2])

                            2                          2\       ]
   + (w[3] v[1] - w[1] v[3])  + (w[1] v[2] - w[2] v[1]) / - b[3]]

> solve(eqs, [w[1], w[2], w[3]]);

[ ]



Maple couldn't find a solution.  The reason is that in general there isn't a solution: your vector is always orthogonal to V, so if B is not orthogonal to V there is indeed no solution.  However:

 

> solve(eqs, [w[1], w[2], v[3]]);

[[                                1                        /     
[[w[1] = - ----------------------------------------------- \b[3]
[[                                 /    2       2       2\       
[[         (b[2] v[2] + b[1] v[1]) \b[1]  + b[3]  + b[2] /       

  /              2       2                      2            \\  
  \v[1] w[3] b[2]  + b[1]  v[1] w[3] + v[1] b[3]  w[3] - b[2]//,

                                  1                        //
  w[2] = - ----------------------------------------------- \\
                                   /    2       2       2\   
           (b[2] v[2] + b[1] v[1]) \b[1]  + b[3]  + b[2] /   

      2                 2                      2            \    
  b[1]  w[3] v[2] + b[3]  w[3] v[2] + w[3] b[2]  v[2] + b[1]/ b[3

   \           b[2] v[2] + b[1] v[1]]]
  ]/, v[3] = - ---------------------]]
                       b[3]         ]]
                                    ]]

Your answer is correct for the one-sided limit as x -> 0+.  In Maple's notation, that would be

limit(x/sqrt(1-cos(x)^2), x = 0, right)

But the one-sided limit as x -> 0- is -sqrt(2), so Maple is correct that the (two-sided) limit is undefined.

Hint: to find the intersection of a line (given parametrically) and a plane (given by an equation in x, y, z), plug the parametric expressions for x, y, z on the line into the equation of the plane, and solve for the parameter t.

Remember that variables are complex by default.  Simplification is tricky where logs of complex quantities are involved, and many of the identities you're used to are not necessarily true.  This case is interesting because it looks like the additional assumptions prevent a desired simplification from taking place, rather than allowing it to take place: without the assumptions on mu and sigma I get (in Maple 13)

mu*((sigma^2+mu^2)^(1/2)/mu)^(-rho), with them I get 

mu*exp(-(1/2)*rho*(ln(sigma^2+mu^2)+2*ln(1/mu))).  The original expression to be simplified is

exp(-ln((sigma^2+mu^2)^(1/2)/mu)*rho+ln((sigma^2+mu^2)^(1/2)/mu)+ln(mu^2/(sigma^2+mu^2)^(1/2)))

I'm not sure exactly what happens here, but I suspect that the assumptions allow some intermediate steps so that the simplification to go along a different path which would not be justified without the assumptions, but ends up with a more complicated expression.

While the is function could certainly be improved, it will never be infallible: there will always be cases where it can't determine whether an inequality is true.  In fact there are theorems that say there is no algorithm to determine whether a function is always positive (even in just one variable, and under some rather restrictive conditions).  See  e.g.  <http://www.mapleprimes.com/questions/36680-Verify-Problem#comment63227>.

It appears to be NULL.

The equation you're trying to solve is

.8778418191 = (int(5.309250137*10^18*sqrt(E)*exp(-.5*E^2/KT^2)*ln(.1006036217*E)/KT^(3/2), 
E = 9.94 .. infinity))/(int(4.409090909*10^18*E^(3/2)*exp(-.5*E^2/KT^2)*ln(0.8771929825e-1*E)
/(KT^(3/2)*(0.233e-1+E)), E = 11.4 .. infinity));

There's no real chance of a symbolic solution, since the integrals can't even be done in closed form. You could try for numerical solutions. But in this case, it looks to me like the right side may always be greater than 1.

I don't know if there's a symbolic solution, but here's how to do it if a[0] ... a[q] and b[m] are given.  I'll try q=3, and let Q = C C*.

>  with(LinearAlgebra):
    q:= 3: Values:= {a[0]=3, a[1]=2, a[2]=1, a[3] = 0, b[m] = -1};
    A:= subs(Values,Matrix(q,q, proc(i,j) if i >= j then a[i-j] else b[m] end if end proc));
B:= subs(Values,Matrix(q,q,proc(i,j) if i > j then b[m] else a[q+i-j] end if end proc));
   Q:= Matrix(q,q,symbol=x, shape=symmetric);
eq:= Q + A . B^%T . Q^(-1) . B . A^%T - (A . A^%T + B . B^%T);
S:= fsolve({seq(seq(eq[i,j],j=i..q),i=1..q)});

  S := {x[1, 1] = 14.50730379, x[1, 2] = 4.953701411,

        x[1, 3] = -1.739291904, x[2, 2] = 13.64941817,

        x[2, 3] = 4.567714781, x[3, 3] = 11.20425848}

> QS:= subs(S, Q);

               [14.50730379     4.953701411    -1.739291904]
               [                                           ]
         QS := [4.953701411     13.64941817    4.567714781 ]
               [                                           ]
               [-1.739291904    4.567714781    11.20425848 ]

Now if the singular value decomposition of C is C = U Sigma V* (where U and V are unitary matrices and Sigma is diagonal and positive semidefinite), C C* = U Sigma^2  U*.  So what we need is that U Sigma^2 U* is the diagonalization of QS.

> E,V:= Eigenvectors(QS);

          [5.39969606970105]
          [                ]
  E, U := [14.4000049895441],
          [                ]
          [19.5612793807549]

        [0.460441223326331 , -0.630480111273029 , 0.624890957808561]

        [-0.625326394755049 , 0.269265793840270 , 0.732436230869434]

        [0.630048236156887 , 0.728004643911506 , 0.270274783431376]
> Sigma:= DiagonalMatrix(sqrt~(E));
                          [2.32372461141613 , 0 , 0]
                          [                        ]
                 Sigma := [0 , 3.79473384963216 , 0]
                          [                        ]
                          [0 , 0 , 4.42281351412818]

V is an arbitrary unitary matrix.  For example, if we take the identity matrix, here's one solution:

> C1:= U . Sigma;
  C1 :=

        [1.06993860275395 , -2.39250421976761 , 2.76377617305221]

        [-1.45308633366043 , 1.02179202243375 , 3.23942886012644]

        [1.46405859273708 , 2.76258386494040 , 1.19537496468836]

Do you mean you want to plot points [x, f(x)] for integers x from a to b?  Make a and b the bounds for the x axis, choose Options, set Style to Point, remove the checkmark from Adaptive Plotting, and set Number of Points to b-a+1.  This doesn't really make x an integer, though: rather, the x values will be floats that happen to be equal to integers.  The distinction is important if f requires integers rather than floats.  In that case, don't use the Plot Builder, use pointplot in the plots package.

If L is a list of n elements with all L[j] >= 0 and sum(L[j], j=1..n) = 1, empirical[L] is a discrete distribution with
Probability(X = j) = L[j] for j in [1,2,...,n], 0 otherwise.  The equivalent in Statistics is ProbabilityTable (not Empirical).  What this Maple code is doing is simulating a Markov chain with transition matrix P and states labelled S,C,R.  I might implement it in a more modern Maple (and in a somewhat more efficient way) like this:


with(Statistics):
P1:= [[.6,.3,.1],[.4,.5,.1],[.3,.4,.3]]:
X:= [seq(Sample(ProbabilityTable(P1[i]),51),i=1..3)]:
R:= Array(1..51):
R[1]:= 1:
for j from 2 to 51 do R[j]:= X[R[j-1]][j] end do:
S:= cat(op(convert(subs(1="S",2="C",3="R",R),list)));

implicitplot with filledregions=true works with a single curve (say f(x,y) = 0), where it plots the regions where f(x,y) > 0 and f(x,y) < 0 in different colours.  For two curves f(x,y) = 0 and g(x,y) = 0, note that f(x,y) < 0 and g(x,y) < 0 is equivalent to max(f(x,y), g(x,y)) < 0.   Thus:

> plots[implicitplot](max(x^2 - y, y - (x + 5)), x = -10 .. 10, y = -10 .. 10, 
filledregions=true, gridrefine=3);

But if you want to use plot, my favourite technique is this.  First determine where the curves intersect.

> solve(x^2 = x + 5);

1/2+1/2*21^(1/2), 1/2 - sqrt(21)/2

> plottools[transform]((x,y) -> [x,y + x + 5])(plot(x^2-x-5,x=1/2 - sqrt(21)/2 .. 1/2 + sqrt(21)/2, 
filled=true));

If you want an implicit function, you shouldn't be assigning something to R.  Instead you want to define an equation.  I'm not sure whether your (m/2R) is supposed to be m/(2*R) or (m/2)*R.  I'll take it as the former.  Note that hirnyk make the other choice.

> eq:= R = exp(L/2 - m/(2*R)); 

Now the question is what you want to do with this.  For example, you could solve explicitly for R

> solve(eq, R, AllSolutions);

Given numeric values for m and L, you could solve numerically for R:

> fsolve(eval(eq, {L=2,m=1}), R=0..1);
   fsolve(eval(eq, {L=2,m=1}), R=1..10);

You could find the derivative of R with respect to L or m using implicit differentiation

> implicitdiff(eq, R, L);

You could plot the surface defined by the equation

> plots[implicitplot3d](eq, R=0.01 ..3,L=0..3,m=0..3,axes=box,grid=[30,30,30],style=patchnogrid,axes=box);


All that the help page says is that it is a CNF, not which particular one it is.  In fact, it seems to be the formula obtained as follows:

(1) negate the formula

(2) convert to DNF by expanding &not and distributing &and over &or

(3) negate again, and expand &not.

You have only one data point (x0, y0, z0).  That's not enough to determine your polynomial, and Maple is telling you that (i.e. for any A,B,C,D,E you can find F so the point is on the curve).  In order to determine all 6 parameters you need 6 data points.  By the way, instead of setting up and solving the equations you might use PolynomialInterpolation from the CurveFitting package.  Also useful is Fit from the Statistics package, for cases where you have more points than parameters and want to find the best fit (rather than a perfect fit).

First 37 38 39 40 41 42 43 Last Page 39 of 138