tomleslie

13876 Reputation

20 Badges

15 years, 181 days

MaplePrimes Activity


These are replies submitted by tomleslie

@torabi 

I'm more comfortable with MATLAB software (maybe true, I wouldn't know. Begs the obvious questions: why did you write in Maple in the first place? why do you have so much trouble rewriting in Matlab?) and MATLAB  in graphical aspect  is more powerful than maple to modify the charts (no it isn't, try to justify this statement somehow, and you will get buried)  and also has a higher quality of the colors (no it doesn't, last time I checked, both handle 24-bit RGB, which is more than enough for any sensible person).

is possible to convert final data (for plotting graph)  to matlab file? Yes, this is trivial to do. Add the command

dataMat:=plottools:-getdata(p1)[3];

to the end of your original Maple worksheet. This will provide a matrix of the data which is actually being plotted. You can then use the ExportMatrix() command to 'export' this matrix to a file in either Mallab binary or ascii formats, which can then easily be read into Matlab

 

You changed the observation period, T, for the function.  Nonetheless, I made those changes in my worksheet & my results concur with yours for CASE 1, but not CASE 2.  So I am still in the same predicament regarding this discrepancy.

Your original worksheet *seemed* to have multiple definitions on this issue. You defined 'T=1'. Sometimes you used [-T/2, T/2] as the period of interest, sometimes [0,T], and (at least once) [-T,T]. So your original actually has three different "definitions" of the period: either [-1/2,1/2], [0,1], or [-1,1]. I decided to use T=1, with period of interest [-T/2, T/2], ie [-1/2, 1/2]. You have to specify what the period of the function is, because I can't do it for you!! I can easily handle the case where the period is [-1,1] because all you have to do in my original worksheet is set T=2 (rather than T=1). Then all relevant caculations will be performed over the time range [-1,1] rather than [-1/2, 1/2]. I have done this in the attached below. Numerical values change (a bit), but Fourier series still agrees with original function over the revised range and Parseval's theeorm still holds

  restart;
#
# OP may not have the OrthogonalExpansions loaded
#
# It is available from
#
# https://www.maplesoft.com/applications/view.aspx?sid=33406
#
  with(OrthogonalExpansions):
#
# Define some parameters
#
  T:=2:
  alpha:=1:
  tau:=1:
#
# Define the "test" function
#
  fx:=piecewise( x<-T/2,
                 undefined,
                 x<=0,
                -(1+x/(alpha*tau))*exp(x/(alpha*tau)),
                 x<=T/2,
                 (1-x/(alpha*tau))*exp(-x/(alpha*tau)),
                 undefined
               );
#
# Compute the 'general' fourier series for the
# function fx
#
  Fseries:= FourierSeries
            ( fx,
              x = -T/2 .. T/2,
              -n..n,
              'Coefficients'
            );
#
# Extract the coefficients (as functions) for this
# general Fourier series. Express as functions, for
# convenient later use
#
  a:=j->evalf(eval(op([2,1,1], Coefficients), i=j));
  b:=j->evalf(eval(op([2,1,2], Coefficients), i=j));

fx := piecewise(x < -1, undefined, x <= 0, -(1+x)*exp(x), x <= 1, (1-x)*exp(-x), undefined)

 

Sum(2*Pi*i*(Pi^2*i^2+2*(-1)^i*exp(-1)-1)*sin(Pi*i*x)/(Pi^4*i^4+2*Pi^2*i^2+1), i = -n .. n)

 

proc (j) options operator, arrow; evalf(eval(op([2, 1, 1], Coefficients), i = j)) end proc

 

proc (j) options operator, arrow; evalf(eval(op([2, 1, 2], Coefficients), i = j)) end proc

(1)

#
# Check Parseval's identity for the function fx,
# using increasing numbers of terms in the Fourier
# series.
#
# Takes ~20secs on my machine, so just wait
#
# Agreement gets pretty good if the number of terms
# is sufficient
#
  evalf((2/T)*int(fx^2, x=-T/2..T/2));
  nterms:= [  10,    20,   50,  100,    200,   500,
             1000, 2000, 5000, 10000, 20000, 50000, 100000]:
  seq( (1/2)*add( a(k)^2+b(k)^2, k=-p..p), p in nterms);
 

.4323323584

 

.3938358159, .4125759754, .4243078225, .4282997890, .4303110026, .4315225996, .4319272762, .4321297667, .4322513096, .4322918320, .4323120948, .4323242528, .4323283056

(2)

#
# As a check that everything has been
# defined more-or-less correctly.
#
# Compute the finite fourier series for
# different numbers of terms, and plot it
# along with the original function.
#
  tr:=[ seq
        ( FourierSeries
          ( fx,
            x = -T/2 .. T/2,
            1..k
          ),
          k=5..20,5
        )
     ]:
   plot([tr[],fx], x=-T/2..T/2);

 

NULL

NULL

Download parseval2.mw

Nonetheless, I made those changes in my worksheet & my results concur with yours for CASE 1, but not CASE 2.  So I am still in the same predicament regarding this discrepancy.

I don't really understand either "CASE1" or "CASE2" in your worksheet. You appear to "type in" various quantities (S1, S2, S3, S4) without any reference to the the fourier series you are trying to compute. Where did these expressions come from? Are they correct? On the other hand the worksheet(s) I have supplied, follow a simple, logical process whihc can be summarised as

  1. Define the (piecewise)  time domain function
  2. Compute its Fourier series representation
  3. Extract the Fourier coefficients from this representation
  4. Use these coefficients to "demonstrate" Parseval's equation
  5. Graph the original function together with its Fourier series representation for various numbers of terms

Now in your original worksheet, try answering the following

  1. Where did S1 come from? It has no *obvious* relation to the Fourier series representation of fx! It just "appears" from nowhere
  2. Where did S2 come from? See (1) above
  3. Do you expect me to believe these expressions are correct? Not unless you prove it!!
  4. Why do you have/need two ways of working out the Fourier series representation? Particularly when they give two different answers. Logically, only one of them can be correct. So which of you calculations is correct, and which one did you screw up?

Last Point

Some time ago I attempted to use the OrthogonalExpansions toolbox per someone else's suggestion.  I am not sure it works with MAPLE 12.  It does look like a nice package to use if I could get it to work.

The OrthogonalExpansions package was created in Maple 13, but (on a quick examination of its code) I would be very, very surprised if it didn't work in Maple 12. It would take all of 10 minutes to download/install and test a few of the examples from its help pages. If they work then you are almost certainly "good to go". Then you can actually run the worksheet(s) I supplied. This seems like a such a simple process, I am completely at a loss to know why you haven't  tried it!

@panke 

The Threads:-Sleep() command was introduced in Maple 15 - and ince you are using Maple 13........

In any case, this would only work if you are deliberately using multithreaded code, whih I'm guessing you aren't

There is no direct equivalent of the Matlab pause() command - I'm not sure why you would want one!.

or anywhere close. Let's take the first few lines of your "Matlab code", with my comments (added in red)

fdsolve := proc({alpha:=NULL, t0:=NULL, t1=NULL, x0=NULL, y0=NULL,  N=NULL}, params)

Matlab does not have a proc() command - try function(). You can use the Matlab command validateAttributes() to check for various poperties of the passed arguments, but I don't think that that you can designate default values using anything similar to the '=NULL' command. You can of course use a variable number of arguments using nargincheck() and varargin()
    local t, h, h1, h2, c, b, x, y, L, n, l, X, Y, f, g;

Matlab does not have a local variable definition. It does however support local and global variables. As a general rule, if you use a name on the left-hand side of an expression, it will be regarded as 'local' to the function.
    eval(F(t,x,y), params);

Matlab does have an 'eval' command, but it bears no relation to the Maple eval() command. Just because the name is the same in the two applications. Don't assume that the commands will do the same thing
     f= unapply(%, [t,x,y]);

Matlab has no 'unapply()' function

So the first four lines of your "Matlab" code contain multiple errors in both in syntax and general understanding. I didn't have the time, inclination or energy to check the rest

I'm still curious that you haven't answered the question I asked in my original response - what is the prupose of this code translation??? Your original code works in Maple: why not just use it

See the attached where I have used the "toy" example where the xlist contains the integers 1..10, and the fxlist contains their squares. This is just so that you can keep track of what is happening. Any two list of the same length would work

  restart;
#
# Create a toy list of x-values and f(x)-values
#
  xlist:= [1, 2, 3,  4,  5,  6 , 7,  8,  9,  10]:
  fxlist:=[1, 4, 9, 16, 25, 36, 49, 64, 81, 100]:
#
# A function which checks whether the supplied
# argument exists in xlist. If it does, output
# the corresponding entry in fxlist. If it doesn't
# then return NULL
#
  f:=p-> `if`( member(p, xlist, 'r'),
               fxlist[r],
               NULL
             ):
#
# Trial values
#
  B:=f(11);
  B:=f(5);
  B:=f(-3);
  B:=f(2.5);

 

25

 

 

(1)

 


 

Download memcheck.mw

as in

A1 := convert(A, Vector); #note capital 'V'

A2 := GivensRotationMatrix(A1, 1, 2);

Try setting Digits:=50 at the top of your worksheet. This will perform all floating-point calculations to 50-Digit precision

A drawback is that, it will also display the result of all calculations to 50 digits, which can be a bit much, so you can set interface(displayprecision=) to something a little more sensible, like 10 maybe, as in the attached. The latter does not affect the number of digits used in cacluclation, just the number used for subsequent display

restart;

Digits:=50:
interface(displayprecision=10);
f := proc (y) options operator, arrow; -(1-tanh(y-.25)^2)^2+2*tanh(y-.25)^2*(1-tanh(y-.25)^2) end proc;
g := proc (y) options operator, arrow; (2*cosh(y-.25)^2-3)/cosh(y-.25)^4 end proc;
f(10);
g(10);

20

 

proc (y) options operator, arrow; -(1-tanh(y-.25)^2)^2+2*tanh(y-.25)^2*(1-tanh(y-.25)^2) end proc

 

proc (y) options operator, arrow; (2*cosh(y-.25)^2-3)/cosh(y-.25)^4 end proc

 

0.27186141816874231204810584779132022747218897494844e-7

 

0.27186141816874231204810584779132022747218915387174e-7

(1)

 

Download precProb.mw

what format your data is in. Just get it into something machine-readable,

You state " it’s all on paper ". Well the bad news is that you are going to have to type it somewhere. It really doesn't matter where.

  1. You could put the data into a txt file, Easiest way to handle this is also the most person-readable, two values per line, 'x' and 'f(x)', separated by a comma and a (few) space(s) if you feel like it (aka CSV). But you could use almost anything else!!!
    1. eg type all the 'x-values', one per line and then all the f(x)-values, one per  line
    2. or type all the x-values in a single line, followed by all the f(x) values in a single line
    3. Don't really care how you do it! If you can explain the layout of the data in the text file then it can be read and formatted for use in Maple.
  2. you could put the data in an excel spreadsheet. Obvious choices  would be one column for x-values and another column for f(x) values. Or one row for x-values and another for f(x values. It really doesn't matter - it can be read!
  3. If the data only exists "on paper" and you are never going to want it anywhere else, then you might as well just type it into Maple rather than a text file which is subsequently imported
    1. Simplest way to do this is probably just to use the Matrix() command and type Matrix([[x, f{x}], [x,f(x)], [x, f(x)],....])
    2. but if you wanted to you could type a list of x-values, like xlist:={xval1, xval2, xval3,  ] and further list of f(x) values, like flist:= [f(xval1), f(xval2), f(xval3),...]
    3. Again, it really doesn't matter what format you choose
    4. If the data exists at all it can be restructured/reformatted into anything required but any subsequent piece of code

The addition operator '+'; is commutative and associative

The multiplication operator '*'; is commutative and associative

Commutativity and associativity are properties of the individual operators: but if you use both in the same expression - then you have to worry about 'precedence'

a+b*c wil l give (with default precedence) a+bc, because '*' has higher precedence than '+'. This is so "obvious" that no-one gives it a second thought. But there might be a 'parallel universe' whose inhabitants have decided that '+' has a higher precedence than '*', in which case a+b*c would evaluate to (a+b)*c.

An analogous argument applies to logical operators: they do have a precedence (and it is largely a matter of "convention"!)

 

 

I finally stopped reading "algorithms" on Wiki and elsewhere, applied wetware and came to the conclusion you are correct for both the application of Cayley-Menger and Tartaglia approach. Checking facial  triangle inequalities is absolutely necessary.

I have added a "do not use" warning to my original response.

I suppose I now have to determine exactly why your 'elegant' response actually works. I've no doubt that it does - but it ain't obvious why (to me at least). Is a geometrical explanation possible?

Apologies for my obtuseness

@vv 

  1. I can find no indication in McCrea's book (reference in revious post), that his method for computing the volume of tetrahedron, requires checking the triangle inequality for faces separately
  2. Similarly in any of the relevant Wikipedia articles on Cayley-Menger determinants/Tartaglia's method can I find any reference to the fact that tthe triangle inequality for faces needs checking separately. See the following references
    1. https://en.wikipedia.org/wiki/Distance_geometry_problem#Cayley%E2%80%93Menger_determinants
    2. https://en.wikipedia.org/wiki/Tetrahedron#Volume
    3. https://en.wikipedia.org/wiki/Niccol%C3%B2_Fontana_Tartaglia
  3. The only mistake I did make in my second post was to assume that a positive determinant required a positive-definite matrix: positive-definiteness is sufficient but not necessary: mea culpa, it was late and I was tired!

For your example test:=[1,1,3,5,1,3], the attached checks for positive tetrahedroal volumes (determinant>0) using both McCrea's matrix and the Cayley-Menger/Tartaglia matrix. Both of these indicate that some tetrahedra can be formed with this list of edges.

Interestingly(?), unless I screwed up somewhere, none of the associated matrices are positive-definite

  restart;
  with(LinearAlgebra):
  with(combinat):
#
# VV test example
#
  test:=[1,1,3,5,1,3]:
#
# Function which computes McCrea determinant
#
  fM1:= ( a, b, c, abar, bbar,cbar )
         ->
         Matrix
         ( 3,
           3,
           [ [      2*a^2,     a^2+b^2-cbar^2, a^2+c^2-bbar^2],
             [ a^2+b^2-cbar^2,      2*b^2,     b^2+c^2-abar^2],
             [ a^2+c^2-bbar^2, b^2+c^2-abar^2,     2*c^2     ]
           ]
         ):

#
# Function for Tartaglia's approach, basically returns
# Cayley-Menger determinant
#
  fM2:= (a, b, c, d, e, f)
        ->
        Matrix
        ( 5,
          5,
          [ [  0,  a^2, b^2, c^2, 1 ],
            [ a^2,  0,  d^2, e^2, 1 ],
            [ b^2, d^2,  0,  f^2, 1 ],
            [ c^2, e^2, f^2,  0,  1 ],
            [  1,   1,   1,   1,  0 ]
          ]
        ):

  g1:= (L, f)-> `if`( Determinant( f(L[]))>0,
                 L,
                 NULL
               ):
  g2:= (L, f)->`if`( IsDefinite( f(L[]),
                                  query=`positive_definite`
                                ),
                     L,
                     NULL
                   ):

#
# Check when the determinant of the McCrea
# matrix is positive

  check1:= seq( g1(k, fM1),
                 k in permute(test, 6)
              );
#
# Check whether the McCrea matrix is positive-
# definite
#
  check2:= seq( g2(k, fM1),
                k in permute(test, 6)
              );
#
# Check when the determinant of the Cayley-
# Menger matrix is positive
#

  check3:= seq( g1(k, fM2),
                 k in permute(test, 6)
              );
#
# Check whether the Cayley-Menger matrix is
# positive definite
#
  check4:= seq( g2(k, fM2),
                k in permute(test, 6)
              );

[1, 1, 3, 5, 1, 3], [1, 1, 3, 1, 5, 3], [1, 3, 1, 5, 3, 1], [1, 3, 1, 1, 3, 5], [1, 3, 5, 1, 3, 1], [1, 5, 3, 1, 1, 3], [3, 1, 1, 3, 5, 1], [3, 1, 1, 3, 1, 5], [3, 1, 5, 3, 1, 1], [3, 5, 1, 3, 1, 1], [5, 1, 3, 1, 1, 3], [5, 3, 1, 1, 3, 1]

 

 

[1, 1, 3, 3, 5, 1], [1, 1, 3, 3, 1, 5], [1, 3, 1, 5, 3, 1], [1, 3, 1, 1, 3, 5], [1, 3, 5, 1, 3, 1], [1, 5, 3, 3, 1, 1], [3, 1, 1, 5, 1, 3], [3, 1, 1, 1, 5, 3], [3, 1, 5, 1, 1, 3], [3, 5, 1, 1, 1, 3], [5, 1, 3, 3, 1, 1], [5, 3, 1, 1, 3, 1]

 

(1)

 

 

Download tetra3.mw

@Ritti 

Your worksheet runs perfectly for me, so I would assume that this is a Maple version issue.

I've tried this on Maple 18 (released in 2014)  and other versions up to and including the current issue, ie Maple 2018.0.

It works in all of them. What version are you running?

@vv 

So let's go through the thought process

  1. there is no indication from the original question that the list of six numbers is in any particular order. Employing  Tartaglia's formula (https://en.wikipedia.org/wiki/Niccol%C3%B2_Fontana_Tartaglia) as you have done, implies an order. You have to designate the vertices of the tetrahedron as 1..4 (you have 24 choices). Having performed this designation, you have to populate the relevant matrix with the side lengths d11 (=0), d12, d13, d14 etc. But you don't know which side length is which - and you really don't want d12 and d34 in the same row/column. So how do you populate the upper triangle of your symmetric matrix given a list of six numbers in no particular order? I decided that in fact there are 120 ways to do this - six possible choices for d12, five for d13, four for d23 and so on. So for a random list of six numbers you have 120 population possibilities -obviously less if two of these numbers happen to be the same, but this gets really "icky" to check for (I know because I tried)
  2. Thus for any list of six numbers, you have to check (worst case) 120 Tartaglia determinants
  3. The determinant of the Tartaglia matrix (reference above) actually computes 288*V^2 where V is the volume of the tetrahedron. So any positive determinant represents a meaningful volume and hence a realisable tetrahedron
  4. The Tartaglia matrix is an extension of the Cayley-Menger matrix (https://en.wikipedia.org/wiki/Distance_geometry_problem#Cayley%E2%80%93Menger_determinants) and if I were using the Cayley-Menger matrix, I accept that I would have to establish the facial triangle inequalities. However I decided not to use this approach, but rather to use the McCrea matrix (see Analytical Geometry of Three Dimensions, by William McCrea, Dover Books, somewhere around p20). This computes the volume of a tetrahedron as 288*V^2=8*determinant_of_McCrea_matrix, without any reference to the facial triangular inequalities. Like the Tartaglia approach above, provided the determinant is positive, the tetrahedron is realiisable
  5. Why did I decide to go the McCrea route rather than Tartaglia? No particullar reasons, other than
    1. in either case I would have to compute 720 determinants for an input list of length six in no known order
    2. the Tartaglia approach requires computing the determinant of 4x4 matrices whereas the McCrea approach requires computing the determinant of 3x3 matrices, and I thought the latter *might* be more efficient, although I doubt if there is much in it
  6. Should I have used IsDefinite(...query='positive_definite') rather than simply checking Determinant()>0 - maybe. If this bothers you excessively, then you can use the revised version attached

  restart;
  with(LinearAlgebra):
  with(combinat):
#
# List to be tested
#
  test:=[2,2,2,2,2,1]:
#
# Function which computes McCrea determinant
#
#  fM:= ( a, b, c, abar, bbar,cbar )
#       ->
#       Determinant
#       ( Matrix
#         ( 3,
#           3,
#           [ [      2*a^2,     a^2+b^2-cbar^2, a^2+c^2-bbar^2],
#             [ a^2+b^2-cbar^2,      2*b^2,     b^2+c^2-abar^2],
#             [ a^2+c^2-bbar^2, b^2+c^2-abar^2,     2*c^2     ]
#           ]
#         )
#       ):
#
# revised function checking directly for positive-definite
#
  fM:= ( a, b, c, abar, bbar,cbar )
       ->
       IsDefinite
       ( Matrix
         ( 3,
           3,
           [ [      2*a^2,     a^2+b^2-cbar^2, a^2+c^2-bbar^2],
             [ a^2+b^2-cbar^2,      2*b^2,     b^2+c^2-abar^2],
             [ a^2+c^2-bbar^2, b^2+c^2-abar^2,     2*c^2     ]
           ]
         ),
         query='positive_definite'
       ):
#
# Check the McCrea determinant for every
# combination of the supplied list, output
# all cases where this is >0
#
# Revised 'if' condition to accommodate the IsDefinite()
# command
#
  seq( `if`( fM(k[]), k, NULL ),
        k in permute(test, 6)
     );

[2, 2, 2, 2, 2, 1], [2, 2, 2, 2, 1, 2], [2, 2, 2, 1, 2, 2], [2, 2, 1, 2, 2, 2], [2, 1, 2, 2, 2, 2], [1, 2, 2, 2, 2, 2]

(1)

 

Download tetra2.mw

Z-differences the wrong way round???

My, unposted response ( cos I took too long) was

restart;
with(LinearAlgebra):
K:=4:
Testna:= RandomMatrix(K-1, 12);
Z:= Vector[row]( K+1,
                              i-> z-2*(i-1)
                          );
CC:= Matrix( K-1,
                      K-1,
                      (i,j) -> add( Testna[ i,
                                                       j+(K-1)*(k-1)
                                                    ]
                                          *
                                          ( Z[k]-Z[k+1] ),
                                          k=1..K
                                       )
                     );

It was only later I noiced that this did not, in fact give the same answers as your code. Discrepancy arises from the successive differences in the Z-vector. The above uses Z[k]-Z[k+1]; the code you posted essentially computes  Z[k+1]-Z[k]. Based on the OP's original question, I *think* the former is correct. So your code(??) should read

K:= 4:
Z:= Vector[row](K+1, n-> z-2*(n-1));
# DZ:= Z[2..] - Z[..-2]; #Delta Z
DZ:= Z[..-2] - Z[2..]; # swap round
CC:= Matrix((K-1)$2, (i,j)-> add(Testna[i, j+(K-1)*(n-1)]*DZ[n], n= 1..K));

 

 

The (incorrect) equation

sum(x^n,n=0..infinity) = 1/(1-x)

Suppose x=2 - in other words

sum(2^n, n=0..infinity) = -1

Now do you believe that?

First 93 94 95 96 97 98 99 Last Page 95 of 207