tomleslie

13876 Reputation

20 Badges

15 years, 163 days

MaplePrimes Activity


These are answers submitted by tomleslie

and it is a nasty one, which I'd like others to confirm before I submit an SCR

Consider the (slightly revsed) version of Kitonum's code, shown below

  restart:
  local D:

  solProc:=proc( P )
                 local Q, f, np, nq, M, N, Sys1, Sys2, Sys3;
                 Q:=x-y-4=0:
                 f:= v-> < coeff(lhs(v), x),coeff(lhs(v), y), coeff(lhs(v), z)>:
                 np:=f(P);
                 nq:=f(Q);
                 M:=[1,0,0]:
                 N:=[0,0,-1]:
                 Sys1:=eval(P,[x,y,z]=~M);
                 Sys2:=eval(P,[x,y,z]=~N);
                 Sys3:=(np.nq)=sqrt(np.np)*sqrt(nq.nq)*cos(Pi/4) assuming real;
                 solve([Sys1, Sys2, Sys3]);
           end proc:

  solProc( A*x+B*y+C*z+D=0);
  solProc( a*x+b*y+c*z+d=0);

Running this code in Maple 2018 returns

{A = 0, B = B, C = 0, D = 0}, {A = -2 B, B = B, C = 2 B, D = 2 B}
{a = 0, b = b, c = 0, d = 0}, {a = -2 b, b = b, c = 2 b, d = 2 b}

But running this code in Maple 2019, Maple 2020 and Maple 2021 returns

                  {A = 0, B = B, C = 0, D = 0}

{a = 0, b = b, c = 0, d = 0}, {a = -2 b, b = b, c = 2 b, d = 2 b}

so for these three (latest) Maple releases, the number of solutions returned by the solve() command depends on whether the parameters in the equations are uppercase or lowercase - now that is nasty

from 'linalg' to 'LinearAlgebra' is shown in the attached.

Apart from the discrepancies noted in the comments in the attached, the two worksheets appear to do the same thing

  restart;
  with(plots):
  unprotect(gamma, D);
  interface(rtablesize = 10): _EnvHorizontalName := 'x': _EnvVerticalName := 'y':
  f := (x, y) -> 4*x^2 + 4*y*x + y^2 - 8*x + 16*y - 17: #(for instance)

  Fg := proc(P::polynom, v::set, V::list, N::list)
           local C, M, i, j;
           C := coeffs(f(x, y), v, M);
           seq(`if`(member(op(i, [M]), N, 'j'), op(j, V) = op(i, [C]), NULL), i = 1 .. nops([M]));
        end proc:
  Fg(f(x, y), {x, y}, [A, B, C, D, E, F], [x^2, y*x, y^2, x, y, 1]):
  assign(%);
  Delta := -4*A*C + B^2:
  var := [x, y]:


  with(LinearAlgebra):
  AA := Matrix([seq([seq(diff(f(x, y), var[i], var[j])/2, j = 1 .. 2)], i = 1 .. 2)]):
  vp := sort([entries(Eigenvalues(AA), 'nolist')]):
  print(`Valeur propres de AA ` = vp):
  P11, DD:=JordanForm(AA, output=['Q', 'J']):
  print(`Matrice diagonale semblable à AA: DD` = DD):
  G := map(Normalize, GramSchmidt([Column(P11, 1 .. 2)]), 2):
  PP := map(simplify, `<|>`(op(G))):
  print(`Matrice de passage orthogonale:   PP` = PP);
  print(`Directions principales de la conique:`):
#
# In the original worksheet with the deprecated linalg()
# package, the col() command returns the relevant columns
# but converts them to rows - no idea why. The corresponding
# Column() command from the LinearAlgebra() package returns
# columns as columns, which seems sensible. Since these are
# only used in the following print() statement, I don't see
# why anyone would worry about the difference. However, if
# this discrepancy is significant, then one could just use
# the Transpose() command, as in the second print() statement
# below
#
  print(`I1 ` = Column(PP, 1), ` J ` = Column(PP, 2)):
  print(`I1 ` = Transpose(Column(PP, 1)), ` J ` = Transpose(Column(PP, 2))):
  alpha := 1/2*arctan(B/(A - C)):
  print('alpha' = evalf(%)):                 

  M1 := Matrix(1, 2, [X, Y]);
  M2 := Matrix(2, 1, [X, Y]);
  M1.AA.M2;
  N := Matrix(1, 2, [D, E]);
#
# In the original worksheet with the deprecated linalg()
# package, the '+' operator doesn't seem to do anything??
# One can circumvent this in the original by using an
# elementwise add (ie +~ ), but in any case, using the
# LinearAlgebra() package, the '+' operator is functional
# so the additions are completed
#
  M1.AA.M2 + N.M2;
  M1.Transpose(PP).AA.PP.M2 + N.PP.M2;
 

`Valeur propres de AA ` = [0, 5]

 

`Matrice diagonale semblable à AA: DD` = (Matrix(2, 2, {(1, 1) = 0, (1, 2) = 0, (2, 1) = 0, (2, 2) = 5}))

 

`Matrice de passage orthogonale:   PP` = (Matrix(2, 2, {(1, 1) = (1/5)*sqrt(5), (1, 2) = (2/5)*sqrt(5), (2, 1) = -(2/5)*sqrt(5), (2, 2) = (1/5)*sqrt(5)}))

 

`Directions principales de la conique:`

 

`I1 ` = (Vector(2, {(1) = (1/5)*sqrt(5), (2) = -(2/5)*sqrt(5)})), ` J ` = (Vector(2, {(1) = (2/5)*sqrt(5), (2) = (1/5)*sqrt(5)}))

 

`I1 ` = (Vector[row](2, {(1) = (1/5)*sqrt(5), (2) = -(2/5)*sqrt(5)})), ` J ` = (Vector[row](2, {(1) = (2/5)*sqrt(5), (2) = (1/5)*sqrt(5)}))

 

alpha = .4636476090

 

Vector[row](2, {(1) = X, (2) = Y})

 

Vector(2, {(1) = X, (2) = Y})

 

Vector(1, {(1) = (4*X+2*Y)*X+(2*X+Y)*Y})

 

Vector[row](2, {(1) = -8, (2) = 16})

 

Vector(1, {(1) = (4*X+2*Y)*X+(2*X+Y)*Y-8*X+16*Y})

 

Matrix(%id = 36893488148080909356)

(1)

 

Download withLA.mw

 

just a lot *simpler* (and works!).

You should use add() not sum() unless you are seeking a "formula" for the summation. I think all of the uneval quotes which you have are redundant.

  restart;

  p:=(i,j)-> local k;
             add( cat(A,__,i,j,k), k=1..3);
  q:= i -> local j;
            add(p(i, j), j = 1 .. 3);
  add(q(i), i = 1 .. 3);

proc (i, j) local k; options operator, arrow; add(cat(A, __, i, j, k), k = 1 .. 3) end proc

 

proc (i) local j; options operator, arrow; add(p(i, j), j = 1 .. 3) end proc

 

A__111+A__112+A__113+A__121+A__122+A__123+A__131+A__132+A__133+A__211+A__212+A__213+A__221+A__222+A__223+A__231+A__232+A__233+A__311+A__312+A__313+A__321+A__322+A__323+A__331+A__332+A__333

(1)

  p:=(i,j)-> local k;
             2*add( cat(A,__,i,j,k), k=1..3);
  add(q(i), i = 1 .. 3);

proc (i, j) local k; options operator, arrow; 2*add(cat(A, __, i, j, k), k = 1 .. 3) end proc

 

2*A__111+2*A__112+2*A__113+2*A__121+2*A__122+2*A__123+2*A__131+2*A__132+2*A__133+2*A__211+2*A__212+2*A__213+2*A__221+2*A__222+2*A__223+2*A__231+2*A__232+2*A__233+2*A__311+2*A__312+2*A__313+2*A__321+2*A__322+2*A__323+2*A__331+2*A__332+2*A__333

(2)

 

Download sums.mw

I don't understand why you are generating ~1000 fits, n=2..3 in 0.001 steps.

Why not just treat 'n' as another fitting parameter restricted to the range 2..3. You then just have one fit to do - see the attached

  restart;
  with(Statistics):
  with(plots):

  j := Matrix(20, 3, [[70, 8, 8.468140006],
                      [70, 10, 4.105515432],
                      [70, 12, 2.36261199],
                      [70, 14, 1.422093923],
                      [70, 16, 0.9],
                      [100, 8, 20.47249229],
                      [100, 10, 9.618450629],
                      [100, 12, 5.360869165],
                      [100, 14, 3.399312905],
                      [100, 16, 2.640399788],
                      [130, 8, 35.90466304],
                      [130, 10, 17.62958097],
                      [130, 12, 9.828362586],
                      [130, 14, 5.866863694],
                      [130, 16, 3.799262645],
                      [160, 8, 57.31814648],
                      [160, 10, 34.49692774],
                      [160, 12, 15.39340528],
                      [160, 14, 9.991012951],
                      [160, 16, 6.049343013]
                     ]
           ):
  expr:= u*h^n*exp(-b*d)/(-c*d + h):
#
# Fit with parameter 'n' restricted to 2 .. 3
#
  ans:= NonlinearFit
         ( expr,
           j,
           [h, d],
           initialvalues = [b = 0.1, c = -30, u = 0.1, n=2.0],
           parameterranges=[n=2..3],
           output = [parametervalues, residualsumofsquares],
           iterationlimit=10000
         );
#
# Plot data points and fit for a "visual check"
#
  display( [ pointplot3d
             ( j,
               symbol=solidsphere,
               symbolsize=16,
               color=black
             ),
             plot3d
             ( eval(expr, ans[1]),
               h=70..160,
               d=8..16,
               color=red,
               style=surface
             )
           ]
        );

[[b = HFloat(0.24364537414299728), c = HFloat(-20.866766308291584), n = HFloat(2.7593674556587215), u = HFloat(0.11018019598352247)], 20.29141052]

 

 

 

Download NLFit.mw

No idea where you got the idea that the intersection point was   [249/100, 471/100, 0].

If I put forward an alternative solution to that already provided by Kitonum, I get the same answer (s)he did

  with(geom3d):
  line(l1, [57/5 + 20*t, -21/5 - 40*t, 20*t], t);
  line(l2, [249/100 - 120*t, 471/100 - 80*t, 200*t], t);
  intersection( p, l1, l2);
  coordinates(p);

l1

 

l2

 

p

 

[933/160, 111/16, -891/160]

(1)

 

Download insect.mw

 

Something wrong with

ExportMattrix( somePathName/someFileName,  someMatrixName )

in the worksheet where the matrix is generated, and then

anyMatrtixName:= ImportMatrix( somePathName/someFileName )

in the worksheet where it is desired?

if I throw out all the random code in your original worksheet, which isn't actually being used, and put what's leftr in more or less legible form, then I have to agree with you - this so-calle iteration process diverges. See the code snip below.

This gives the same answer as your original worksheet (ie 0.8789295492), for a single iteration (ie m=1), but diverges rapidly as the iteration count increases.

  restart:
  u[0]:= x*(1/2*exp(-t/4) + exp(-1/2) - 1) + 1:
  alpha:= 1;
  printf(  "    Number of\n"):
  printf(  "    Iterations         answer\n"); 
  for m from 1 to 10 do
      for n from 0 to m do
          u[n + 1]:= simplify
                     ( u[n]
                       +
                       alpha*int( x*(1-s)*(diff(u[n], t)-diff(u[n], x$2)-1/4*u[n]),
                                  s = 0..x
                                )
                       +
                       alpha*int( s*(1-x)*(diff(u[n], t)-diff(u[n], x$2)-1/4*u[n]),
                                  s = x..1
                                )
                     ):
       end do:
       printf( "       %2d      %16.8f\n",
               m,
               evalf(eval(u[n], [x = 0.2, t = 0.2]))
             );
  end do:
  
                           alpha := 1

    Number of
    Iterations         answer
        1            0.87892955
        2            0.95670244
        3            0.66722716
        4            2.78900964
        5          -18.82567148
        6          228.99961580
        7        -3004.18166800
        8        47074.66315000
        9      -875743.22920000
       10      18798150.48000000

mann.mw

From this information there are only two possible conclusions

  1. The iteration process has been correctly implemented but does not converge - who claims that it does?
  2. The iteration process should converge, but won't because it has been incorrreclty implemented - and I'm not going to read a few academic papers to work out how to implement it properly!

Pick one of the above - or come up with an alternative!

This image

was generated using the code below (and in attached link). It is 501*501 pixels and (on my machine) took about 90secs to produce. If you have the patience then just change the variable "npts" - a 1001x1001 image ought to be available in ~360secs

  restart:
  with(ImageTools):
  st:=time():
#
# Define the "ranges" for the complex variable
# starting point, and the numbe of points in the
# image. (Actually the number of points in the
# image will be npts+1 by npts+1
#
  bl:=[-2.0, -2.0]:
  ur:=[2.0, 2.0]:
  npts:= 500:
  step:= (ur[1]-bl[1])/npts:
#
# Initialize the "RGB" array
#
  A:=Array(0..npts, 0..npts,1..3):
#
# Specify the three roots. (These are given as
# "floats" just to make convergence testing
# simpler, later on
#
  rts:=[ 1.0+0*I,
        -0.5+evalf(sqrt(3)/2*I),
        -0.5-evalf(sqrt(3)/2*I)
       ]:
#
# Specify the rgb sub-arrays
#
  rgb:= [ Array(1..3, [1,0,0]),
          Array(1..3, [0,0,1]),
          Array(1..3, [0,1,0])
        ]:
#
# Run the calculation: the outer two loops just
# generate the starting value across the range
# -2-2*I to 2+2*I. Yje case 0+0*I has to be trapped!
# 
  for i from 0 by 1 to npts do
      for j from 0 by 1 to npts do
          zold:=bl[1]+j*step+(bl[2]+i*step)*I;
          if   zold=0+0*I
          then break
          else zold:=evalf(zold);
             #
             # Run the iteration
             #
               while true do
                     znew:=zold-(zold^3-1)/3/zold^2;
                     if   abs(znew-zold) < 1e-06
                   #
                   # if true, convergence has been achieved
                   # so check which root is obtained and
                   # populate the output arra with the
                   # appropriate sub-array, and brak from
                   # this inner loop
                   #
                     then A[i,j,...]:= rgb[min[index](abs~(rts-~znew))];
                          break;
                     else zold:=znew;
                     fi;
               od;
          fi;
      od;
   od;
#
# Check time
#
   time()-st;
#
# Display image
#
   pic:=Create( npts+1, npts+1, A):
   Embed(pic);
#
# Write out image - path name will need to be changed
# depending on user
#
  Write("C:/Users/TomLeslie/Desktop/basin.jpg", pic):
                             87.703

basins.mw

There are a few problems with this worksheet

  1. A couple of simple typos whihc I have corrected (and indicated) in the attached.
  2. After fixing these, dsolve the produces a solution module, but it cannot be evaluated, beyond the initial time. My thought was that thsi might be related to the use of piecewise functions which rely on dependent variables, so I replaced these with hs1:=, hs2:=1, hd3:=1 in the attached. This produces a graph - but not for long!
  3. the stop_cond option which you uses in the dsolve command triggers at around t=5e-03, so I have graphe all dependent variables up to thhis point
  4. I think if you want to reinstate the" piecewise behaviour", you are probably going to have to use the 'events' option with discrete variables. See the help at ?Events for dsolve[numeric]. I have added a link to this in the attached. I warn you, that (in my experience of using events, particularly with discrete variables), this is a non-trivial exercise, and possibly best practised on a "simpl"e ODE system

Good luck with the attached, which for some reason won't display inline on this site - becuase somebody has broken something AGAIN!

odeProb.mw

See the attached for the simple plot in your question.

For future reference, please use the big-green up-arrow in the Mapleprimes toolbar to upload worksheets, rather then "cutting-and pasting"

So else what could be wrong?

  1. Wrong Maple version for your OS? Unlikely but possible. The attached shows how to use the interface(version) command to return bith you OS and the Maple version you are running.
  2. Installation problem - try reinstalling, but not before you have read to the bottom of this list for other possibilities
  3. "Legacy" .init file. If your previous Maple installation had some "clever" stuff in an .init file this *may* be an issue, because installing an upgrade will probably copy your legacy .init file
  4. Licensing issues. It has often been reported ithat if some aspect of you licensing is incorrect, then (unlike many other software packages, whihc will just refure to run), Maplee will actuall "run", but just produce completely erroneous output. I think this is just the method used by those at Maplesoft to annoy anyone trying to pirate the software

I would display the output of the attached worksheet inline on this site, but once again some lamebrain at Maple has "improved" thsi site so that simple inline display of trivial worksheets is no longer possible

Download simpleplot.mw

 

 

Probably several ways to do this, but the attached is the "simplest" I can come up with - basically divide the current time (t) by the "period" (T) and if the remaindet is less than "duty", signal is ON, otherwise signal is OFF.

Check the attached


 

  restart;
  with(plots):

  f1 := (x, t ) -> A * (1 + B * sin(a*x-b*t));

proc (x, t) options operator, arrow; A*(1+B*sin(a*x-b*t)) end proc

(1)

#
# Apply a "duty cycle" to f2
#
  f2 := (x, t, period, duty) -> `if`( frem(t, period) < duty,
                                      C * (1 + E * sin(c*x-d*t)),
                                      0
                                     ):

  A:=1: B:=0.5: a:=100: b:=1:
  C:=1: E:=0.5: c:=10: d:=10:
#
# This should reproduce OP's original, ie the
# with f2 always ON
#
  T:=1:
  q:=1:
  densityplot('f1(x, t) * f2(x, t, T, q)',
               x = 0..1,
               t = 0..10,
               style=patchnogrid
             );
#
# This shoudd turn the signal f2 ON 1/3 of the
# time and OFF 2/3 of the time
#
  T:=1:
  q:=1/3;
  densityplot('f1(x, t) * f2(x, t, T, q)',
               x = 0..1,
               t = 0..10,
               style=patchnogrid
             );

 

1/3

 

 

 


 

Download duty.mw

of which the attached (using the geometry() package) is just one

  restart:
  with(geometry):
  _EnvHorizontalName:=x:
  _EnvVerticalName:=y:
  parabola(p1, x^2-6*x+3*y+18=0);
#
# Get the details of the parabola
#
  detail(p1);
  draw(p1);
#
# Extrct the coordinates of the focus and
# the equation of the directrix
#
  coordinates( focus(p1));
  Equation(directrix(p1));
#
# The axis of the parabola is the line
# perpendicular to the directrix which goes
# through the focus, so construct this
# perpendicular and return its equation
#
  Equation
  ( PerpendicularLine
     ( ax,
       point( f,
              coordinates
              ( focus(p1)
              )
            ),
       line( d1,
             Equation
             ( directrix(p1)
             )
           )
     )
  );

p1

 

GeometryDetail(["name of the object", p1], ["form of the object", parabola2d], ["vertex", [3, -3]], ["focus", [3, -15/4]], ["directrix", y+9/4 = 0], ["equation of the parabola", x^2-6*x+3*y+18 = 0])

 

 

[3, -15/4]

 

y+9/4 = 0

 

x-3 = 0

(1)

 

Download pbola.mw

 

 

for log(), which states (emphasis added)

For complex-valued expressions x,  ln(x) = ln(abs(x)) + argument(x)*I, where  -Pi < argument(x) <= Pi. Throughout Maple, this computation is taken to be the definition of the principal branch of the logarithm.

So when you take the imaginary part of  the logarithim of a complex number, the answer will always lie within the range -Pi..Pi

- just because I like alternatives

  restart;
  with(Student[Precalculus]):
  with(plots):
  pt:=[0.504244923, 0.3781836925]:
  f:=7*x^2 - 7*y(x)^2 - 12.0*x + 9.0*y(x) + 2.25 - 2*x*y(x)=0:
  myLine:=Line( pt,
                eval
                ( rhs
                  ( isolate
                    ( diff(f, x), diff(y(x), x))
                  ),
                  [ x=pt[1], y(x)=pt[2]]
                )
              )[1];
  display( [ implicitplot
             ( f,
               x=-1..2,
               y=-2..2,
               color=red
             ),
             plot
             ( rhs(myLine),
               x=-1..2,
               color=blue
             ),
             pointplot
             ( pt,
               color=green,
               symbol=solidcircle,
               symbolsize=16)
           ],
           scaling=constrained,
           view=[-1..2, -2..2]
         );

y = 2.112372436*x-.6869693835

 

 

 

Download hypLine.mw

but I suspect that the existence of the print() statement means that your procedure cannot "autocompile" successfully so when you invoke it, you are actually running "native" Maple code.

When you comment out the print() statement, Maple does its best to compile the procedure (and may even think it has succeeded), but as a general rule, you cannot "compile" commands from a Maple package (such as combinat) - so the "compiled" version will not run.

So your problem is more with the statement

options threadsafe,autocompile:

than anything else.

Note that the attached runs correctly with/without the print statement

  restart;
  fib_even_sum := proc()::integer;
                      uses combinat:
                      local val:= 0, n:= 1, FN:= 0:
                      print(%);
                      while FN <= 4000000 do:
                            if   type(FN,even)
                            then val += FN:
                            fi:
                            n++:
                            FN := fibonacci(n):
                      end do:
                      return val:
                      end proc:
  fib_even_sum();

0

 

4613732

(1)

 restart;
  fib_even_sum := proc()::integer;
                      uses combinat:
                      local val:= 0, n:= 1, FN:= 0:
                     # print(%);
                      while FN <= 4000000 do:
                            if   type(FN,even)
                            then val += FN:
                            fi:
                            n++:
                            FN := fibonacci(n):
                      end do:
                      return val:
                      end proc:
  fib_even_sum();

4613732

(2)

 

Download fib.mw

 

First 31 32 33 34 35 36 37 Last Page 33 of 207