Advanced Search

tomleslie

13876 Reputation

20 Badges

15 years, 163 days
Contact tomleslie

MaplePrimes Activity


These are answers submitted by tomleslie

Hmmmm...

tomleslie 13876
June 16 2022
1 0

looks like Python is producing some sort of listlist format. I can see two possibilities

  1. Persuade Pythin to output data in a 'csv' format, mainly because for a listisit, Maple likes "commas" between list entries and "commas" between lists. As a general rule, almost any software can produce "csv" format, and almost any software can read it
  2. If this is not posssible (or too difficult|) then using the big green up-arrow in the Mapleprimes toolbar, upload the file which Python outputs, which *ought* to some sort of human-readable text, and then we can figure out whattsort of edits might be needed by a Maple read operation

Not exactly sure...

tomleslie 13876
June 15 2022
1 2

what you mean by "isolate"

Maybe as shown in the attached?

> 

  restart;
  expr:= diff(u(t,x,v),t)+v*diff(u(t,x,v),x)+diff(u(t,x,v),x)*int(u(t,x,v)*v*exp(x),v)-v*diff(u(t,x,v),v)*sin(v)*int(v*diff(u(t,x,v),x),v)=0;
  isolate(expr, sin(v));
  isolate(IntegrationTools:-Expand(expr), exp(x));

diff(u(t, x, v), t)+v*(diff(u(t, x, v), x))+(diff(u(t, x, v), x))*(int(u(t, x, v)*v*exp(x), v))-v*(diff(u(t, x, v), v))*sin(v)*(int(v*(diff(u(t, x, v), x)), v)) = 0

  

sin(v) = -(-v*(diff(u(t, x, v), x))-(diff(u(t, x, v), x))*(int(u(t, x, v)*v*exp(x), v))-(diff(u(t, x, v), t)))/(v*(diff(u(t, x, v), v))*(int(v*(diff(u(t, x, v), x)), v)))

  

exp(x) = (-v*(diff(u(t, x, v), x))-(diff(u(t, x, v), t))+v*(diff(u(t, x, v), v))*sin(v)*(int(v*(diff(u(t, x, v), x)), v)))/((diff(u(t, x, v), x))*(int(u(t, x, v)*v, v)))

 (1)
> 

 

Download isol.mw

Hmmmm...

tomleslie 13876
June 15 2022
2 0

Since yoiur function delta(s) appears to define a square wave, there are far better ways of computings its integral, rather than a brute-force numerical integration. I have show one possibility in the attached and compared its result with the "brute force" numerical integration. It is both faster and more accurate.

Since all other quantities are "known", it is not then too difficult to construct the function x(t) and plot it. This is done in the attached for arbitrary values of a, p, q, r, phi - since OP does not provide these

> 

  restart;
  ulim:=124.75:
#
# Define the function delta, suitable for use with
# the int() command
#
  delta:= proc( s )
                if   type(s, numeric)
                then if   type(floor(s), even)
                     then return 1;
                     else return 0;
                     fi;
                else return 'thisproc'(s);
                fi;
          end proc:
#
# Numerically integrate the function delta over the range 0..ulim
#
# Note that this is *BAD* way to compute this quantity
#
  evalf( Int(delta, 0..ulim));

62.70582111

 (1)
> 

#
# The functiondelta() above defines a "square wave". The integral
# of a square wave can be calculated in many other ways, of which
# the following is just one
#
  intDelta:= proc( t )
                   if   type(t, numeric)
                   then return floor((t+1)/2)
                               +
                               `if`( type(floor(t), even),
                                     t-floor(t),
                                     0
                                   ):
                   else return 'thisproc(_passed)';
                   fi:
             end proc:
  intDelta(ulim);

62.75

 (2)
> 

#
# So now OP can define the function x(t) as
#
  x:= proc(tau)
           if  type(tau, numeric)
           then return exp(-p*tau)*( p*q*tau+ r*intDelta(tau)+phi+a*(exp(p)-1)/p/(exp(p)+1));
           else return 'thisproc(_passed)';
           fi;
      end proc:
#
# which can be evaluated for any value of the passed argument
#
  x(ulim);
#
# which in turn can be evaluated for any supplied values of
# a, p, q, r, phi and so can be plotted
#
  eval( x(ulim), [a=1.0, p=2.0, q=3.0, r=4.0, phi=5.0]);
#
# and so can be plotted
#
 a:=1.0: p:=2.0: q:=3.0: r:=4.0: phi:=5.0:
 plot( x(z), z=0..10);

exp(-124.75*p)*(124.75*p*q+62.75*r+phi+a*(exp(p)-1)/(p*(exp(p)+1)))

  

0.4422229855e-105

  

 
> 

 

> 

 

> 

intDelta(0.5);

.5

 (3)
> 

 

Download sqWave.mw

It depends a bit...

tomleslie 13876
June 15 2022
4 0

on exaclt ywhat you want to appear in the figure eg just the intersection lines, or just the plane and the cube, or maybe all three?

See the attached (NB the plot renders a lot better in a Maple worksheet than it does on this site)

> 

  with(plots):
  with(plottools):
  p1:= cuboid( [0,0,0],
               [1,1,1],
               color=blue,
               transparency=0.5
             ):
  p2:= implicitplot3d( x+z=1,
                       x=-1..2,
                       y=-1..2,
                       z=-1..2,
                       style=surface,
                       color=red
                     ):
  l1:=line( [1,0,0], [1,1,0], color=yellow, thickness=10):
  l2:=line( [1,0,0], [0,0,1], color=yellow, thickness=10):
  l3:=line( [1,1,0], [0,1,1], color=yellow, thickness=10):
  l4:=line( [0,0,1], [0,1,1], color=yellow, thickness=10):
  display([p1,p2, l1,l2,l3,l4]);

 
> 

 


 

Download iplot2.mw

Free advice...

tomleslie 13876
June 14 2022
1 2

You state

I want to somehow import two arrays and an equation from maple to said python script, use them in a procedure there, and then import the output of said procedure, a 2d array, backwards to maple.

Don't!!!

Either

  1. Do everything in Python
  2. Do everything in Maple

As a general rule the pain of switching between major software packages just isn't worth it

Rather than...

tomleslie 13876
June 14 2022
1 0

trying to persuade Maple to show you its step-by-step calculation, why not write your own? That way you know exactly what is being calculated, and how. It really is pretty trivial - see the attached

> 

  restart;
#
# The trapezoidal rule
#
# See https://en.wikipedia.org/wiki/Trapezoidal_rule
#
  trapRule:= x-> local j;
                 (x[1]-x[0])/2*(   f(x[0])
                                 + add
                                   ( 2*f(x[j]),
                                     j=1..op([2,2], x)-1
                                   )
                                 + f( x[op([2,2], x)])
                               ):
#
# Simpson's "1/3-rule"
#
# See https://en.wikipedia.org/wiki/Simpson%27s_rule
#
  simpRule:= x-> local j;
                 (x[1]-x[0])/3*add
                               (   f(x[2*j-2])
                                 + 4*f(x[2*j-1])
                                 + f(x[2*j]),
                                 j = 1..op([2,2], x)/2
                               ):
#
# Limits of integration and integral to be approximated
#
  a:= -1:
  b:= 2:
  f:= x-> evalf(x^2):
#
# Number of intervals
#
  numIntervals:= 4:
#
# x-values
#
  xv:= Array
       ( 0..numIntervals,
         j-> a+j*(b-a)/numIntervals
       );
#
# Apply above trapezoidal rule to these x-values
#
  trap:= trapRule( xv );
#
# Apply above Simpson's rule to these x-values
#
  simp:= simpRule( xv );
#
# Exact result
#
  exact:= int(f(x), x=a..b);

`Array(0..4, {(1) = -1, (2) = -\`/\`(1, 4), (3) = \`/\`(1, 2), (4) = \`/\`(5, 4)})`

  

3.281250000

  

3.000000000

  

3

 (1)
> 

 

Download trapSimp.mw

Just -- overkill??...

tomleslie 13876
June 12 2022
1 0

The attached shows two "simple" ways

  1. a loop which will keep executing until a preset error tolerance setting is achieved - has the advantage of showing the intermediate approximations
  2. a recursive procedure call which keep going until the erro tolerance is achieved - only the final root approximation is available
> 

  restart:
  f:= x-> x^3+x+1:
  xold:=-1.0:
  tol:= 1e-09:
  while true do
      xnew:=eval(x-f(x)/diff(f(x),x), x=xold):
      if   abs(xold-xnew) < tol
      then break;
      else xold:=xnew:
      fi:
  od;
  xnew;

-.7500000000

  

-.6860465116

  

-.6823395824

  

-.6823278040

  

-.6823278040

  

-.6823278040

 (1)
> 

  restart:
  f:= x-> x^3+x+1:
  tol:= 1e-09:
  rt:=proc(guess)
           local x:
           x:=  eval(x-f(x)/diff(f(x),x), x=guess):
           if   abs(guess-x) < tol
           then return x;
           else rt(x):
           fi:
     end proc:
  rt(-1.0);

-.6823278040

 (2)
> 

 

Download getRoot.mw

Even simpler...

tomleslie 13876
June 12 2022
1 0

is to realize that you need to write procedures for this: just use Maple's sort() command with appropriate optione, as shown in the attached

> 

  restart;
  vals:=[ a = -1/3, b = 0, c = 7/3, d = -2];
#
# sort() vals, biggest to smallest
#
  biggest:=  vals[ sort
                   ( rhs~(vals),
                     `>`,
                     output=permutation
                   )
                 ];
#
# sort() vals, smallest to biggest
#
  smallest:= vals[ sort
                   ( rhs~(vals),
                    `<`,
                     output=permutation
                   )
                 ];

[a = -1/3, b = 0, c = 7/3, d = -2]

  

[c = 7/3, b = 0, a = -1/3, d = -2]

  

[d = -2, a = -1/3, b = 0, c = 7/3]

 (1)
> 

 

Download simpleSort.mw

 

Confirming previous answers...

tomleslie 13876
June 11 2022
0 1

You need to define the coordinate system, and there are a couple of ways to do this.

Pick either of the ways shown in the attached

> 

  restart:
  with(geometry):
  _EnvHorizontalName := x:
  _EnvVerticalName := y:
  circle(c, (x+4)^2+(y-3)^2=81);
  draw(c);

c

  

 
> 

  restart:
  with(geometry):
  circle(c, (x+4)^2+(y-3)^2=81, [x,y]);
  draw(c);

c

  

 
> 

 

Download geomCirc.mw

 

Well...

tomleslie 13876
June 10 2022
4 3
  1. I can't check this in Maple 12, which is now close to 15 years old. So I checked this in the current release (Maple2022)
  2. You define four parametric curves as[ xi(r), yi(r), r=-10..10]. Your code plots these four parametric curves. There are however a couple of oddities
  3. The first and second parametric curves are just a reflection of each other in the x-axis as are the third and fourth. In other words x1(r)=x2(r) and y1(r)=-y2(r).  Similarly x3(r)=x4(r) and y3(r)=-y4(r). 
  4. Each of the four parametric curves has a "gap". This is because in the range (approximately) r=-2..2, all of the expressions for yi(r) generate complex values. You can check this by plotting each of the yi(r) as a simple function of 'r'

See the attached for more detail

> 

  restart;
  expr:=[ [ (1.428571429*(r^2+.49-4*r*(r^2+.49-2*exp(-.1/r))/(2*r-.2*exp(-.1/r)/r^2)))*csc((1/3)*Pi),
          sqrt(16*r^2*(r^2+.49-2*exp(-.1/r))/(2*r-.2*exp(-.1/r)/r^2)^2+.49*cos((1/3)*Pi)^2-2.040816327*(r^2+.49-4*r*(r^2+.49-2*exp(-.1/r))/(2*r-.2*exp(-.1/r)/r^2))^2*cot((1/3)*Pi)^2)
        ],
        [(1.428571429*(r^2+.49-4*r*(r^2+.49-2*exp(-.1/r))/(2*r-.2*exp(-.1/r)/r^2)))*csc((1/3)*Pi),
          -sqrt(16*r^2*(r^2+.49-2*exp(-.1/r))/(2*r-.2*exp(-.1/r)/r^2)^2+.49*cos((1/3)*Pi)^2-2.040816327*(r^2+.49-4*r*(r^2+.49-2*exp(-.1/r))/(2*r-.2*exp(-.1/r)/r^2))^2*cot((1/3)*Pi)^2)
        ],
        [(5.*(r^2+0.4e-1-4*r*(r^2+0.4e-1-2*exp(-.2/r))/(2*r-.2*exp(-.2/r)/r^2)))*csc((1/3)*Pi),
          sqrt(16*r^2*(r^2+0.4e-1-2*exp(-.2/r))/(2*r-.2*exp(-.2/r)/r^2)^2+0.1e-1*cos((1/3)*Pi)^2-25.*(r^2+0.4e-1-4*r*(r^2+0.4e-1-2*exp(-.2/r))/(2*r-.2*exp(-.2/r)/r^2))^2*cot((1/3)*Pi)^2)
        ],
        [(5.*(r^2+0.4e-1-4*r*(r^2+0.4e-1-2*exp(-.2/r))/(2*r-.2*exp(-.2/r)/r^2)))*csc((1/3)*Pi),
         -sqrt(16*r^2*(r^2+0.4e-1-2*exp(-.2/r))/(2*r-.2*exp(-.2/r)/r^2)^2+0.1e-1*cos((1/3)*Pi)^2-25.*(r^2+0.4e-1-4*r*(r^2+0.4e-1-2*exp(-.2/r))/(2*r-.2*exp(-.2/r)/r^2))^2*cot((1/3)*Pi)^2)
        ]
       ]:
> 

#
# Do some checking on these expressions
#
# Check "x-functions'
#
  is(expr[1,1]= expr[2,1]);
  is(expr[3,1]= expr[4,1]);
#
# Check 'y-functions'
#
  is(expr[1,2]= -expr[2,2]);
  is(expr[3,2]= -expr[4,2]);

true

  

true

  

true

  

true

 (1)
> 

#
# From the above, the first and second parametric
# curves will the have the same "x-value" and
#  "y-values" of opposite sign, so the second
# parametric curve is identical to the first
# parametric curve exceptfor a reflection in the x-axis.
#
# The same applies to the third and fourth parametric
# curves - identical except for reflection in x-axis
#
plot( [seq( [expr[j][], r=-10..10], j=1..4)],
        -13 .. 5,
        -5 .. 5,
        color=[red, blue, green, black]);

 
> 

#
# Why is there a "gap" in each of these four curves?
#
# Plot the 'y-expressions' on their own
#
# Check the value of these expressions at say r=-1 and r=1.
# These are all complex
#
  seq( seq( eval( expr[j,2], r=-i),i in[-1,1]), j=1..4);

  plot( [seq( expr[j,2], j=1..4)],
        r=-10..10,
        view=[-13 .. 5, -5 .. 5],
        color=[red, blue, green, black]
      );

2.166736088*I, 2.738430349*I, -2.166736088*I, -2.738430349*I, 6.965649860*I, 10.43539041*I, -6.965649860*I, -10.43539041*I

  

 
> 

 

Download parPlot2.mw

When you say "solutions"...

tomleslie 13876
June 09 2022
3 3

do you mean solutions to (x-y)^2+(1-z)^2 =0?  Since both terms on the right hand side are quadratic (ie >=0), then the only solutions occurs when each of these terms is zero - in other words y=x, z=1. You can plot this using spacecurve(), as shown in the attached

> 

  restart;
  with(plots):
  expr:=(x-y)^2+(1-z)^2;
  eval(expr, [y=x, z=1]);
  spacecurve([x,x,1], x=-10..10, color=red, thickness=5);

(x-y)^2+(1-z)^2

  

0

  

 
> 

 

Download spCurve.mw

Are you doing this the hard way...

tomleslie 13876
June 08 2022
2 1

for a reason?. Why not just use dsolve() with the series option, as in the attached. Or am I missing something?

> 

ode := xi^2*(diff(psi[m](xi), xi, xi))+2*xi*(diff(psi[m](xi), xi))+(xi^2-m*(m+1))*psi[m](xi) = 0; Order := 8; dsolve(ode, psi[m](xi), series); convert(%, polynom)

xi^2*(diff(diff(psi[m](xi), xi), xi))+2*xi*(diff(psi[m](xi), xi))+(xi^2-m*(m+1))*psi[m](xi) = 0

  

psi[m](xi) = _C1*xi^m*(series(1+(1/(-4*m-6))*xi^2+(1/((-8*m-20)*(-4*m-6)))*xi^4+(1/((-12*m-42)*(-8*m-20)*(-4*m-6)))*xi^6+O(xi^8),xi,8))+_C2*xi^(-m-1)*(series(1+(1/(m^2-(-m-1)^2+6*m-1))*xi^2+(1/((m^2-(-m-1)^2+10*m-11)*(m^2-(-m-1)^2+6*m-1)))*xi^4+(1/((m^2-(-m-1)^2+14*m-29)*(m^2-(-m-1)^2+10*m-11)*(m^2-(-m-1)^2+6*m-1)))*xi^6+O(xi^8),xi,8))

  

psi[m](xi) = _C1*xi^m*(1+xi^2/(-4*m-6)+xi^4/((-8*m-20)*(-4*m-6))+xi^6/((-12*m-42)*(-8*m-20)*(-4*m-6)))+_C2*xi^(-m-1)*(1+xi^2/(m^2-(-m-1)^2+6*m-1)+xi^4/((m^2-(-m-1)^2+10*m-11)*(m^2-(-m-1)^2+6*m-1))+xi^6/((m^2-(-m-1)^2+14*m-29)*(m^2-(-m-1)^2+10*m-11)*(m^2-(-m-1)^2+6*m-1)))

 (1)
> 

``

Download odeSer.mw

My $0.02...

tomleslie 13876
June 08 2022
4 4

I think that the only thing you are really missing is the fact that when you write omega=2*Pi*f1, you have to realize that the combined term (2*Pi) should have the units of 'radians', and you have to insert this

For reasons which are slightly(?) mysterious (to me at least) Maple will return the resut with the unit m/(s*m('radius')) - wihc I admit *looks* pretty weird, but is actually correct! Consider the simple-minded proble of computing arclength in terms of raidua and angle. One would normally write s=r*theta where s is arclengt in (say) metres, r is angle in (say) metres and theta is angle in radians. So if you rearrange this as theta=s/r, then one has angle(radians)=arclength(metres)/(radius(metres). Now the problem with this representation is that the right-hand-side is "unitless", because it is expressed in metres/metres - yet the left-hand side has units of radians. So Maple circumvents this logical contradiction by experessing the right-hand-side as metres/metres(radius). It is a method of overcoming the unit cancellation which would otherwise occur.

Although this workaround looks a little clumsy, it is not difficult to convert your result for angular speed to either rad/sec or rpm - see the attached.

NB on this site, anything with units still seems to appear in double bracket notation (ie [[unit]]). This does not happen in the actual worksheet- honest!

> 

restart
> 

with(Units[Standard])
> 

 

> 

acc1 := 9.0*Unit('gn')

9.0*Units:-Unit(gn)

 (1)
> 

arm1 := 10.0*Unit('ft')

10.0*Units:-Unit(ft)

 (2)
> 

vel1 := (acc1*arm1)^(1/2)

16.40170792*Units:-Unit(m/s)

 (3)
> 

circ1 := 2*Pi*arm1

62.83185308*Units:-Unit(ft)

 (4)
> 

T1 := circ1/vel1

1.167631378*Units:-Unit(s)

 (5)
> 

f1 := 1/T1

.8564346752*Units:-Unit(1/s)

 (6)
> 

omega1 := 2*Pi*Unit('rad')*f1;
convert(omega1, units, rpm)
convert(omega1, units, radians/sec)
 

5.381137767*Units:-Unit(m/(s*m(radius)))

  

51.38608049*Units:-Unit(rpm)

  

5.381137767*Units:-Unit(rad/s)

 (7)
> 

 

Download unitProb.mw

It's simpler than you think...

tomleslie 13876
June 06 2022
2 0

see the attached. (The plot renders better in the worksheet than it does on this site -honest!)

> 

  restart:
  with(plots):
  theta := (2*Pi)/17:
  pointplot
  ( [ seq
      ( [cos(k*theta), sin(k*theta)],
        k=1..16
      )
    ],
    symbol = solidcircle,
    color = red,
    symbolsize = 10
  );

 
> 

 

Download ptplt.mw

Your basic problem...

tomleslie 13876
June 06 2022
2 0

is the use of the entry alpha*gamma=0.004 in the list parameters, and in particular when you use this to evaluate terms in the ODEs. Although your ODEs only ever contain the product alpha*gamma (and never either of these parameters separately), you cannot use the eval() command to "substitute for a product term. The simplest way around this is to write alpha=0.004/gamma.

I tidied up a fair bit, because there was a lot of "stuff" in this worksheet which wasn't being used, and a lot of things which could be done in simpler(?) ways

> 

restart
> 

with(plots)
> 

U := W(t), X(t), Y(t), Z(t); funcs := proc (W, X, Y, Z) options operator, arrow; A-mu*W-delta1*Y-delta2*Z end proc, proc (W, X, Y, Z) options operator, arrow; beta*(W-X-Y-Z)*X-mu*X-X end proc, proc (W, X, Y, Z) options operator, arrow; alpha*gamma*X-(mu+delta1)*Y end proc, proc (W, X, Y, Z) options operator, arrow; (-alpha*gamma+1)*X-(mu+delta2)*Z end proc; sys := seq(diff(U[j], t) = funcs(U)[j], j = 1 .. 4)

proc (W, X, Y, Z) options operator, arrow; A-mu*W-delta1*Y-delta2*Z end proc, proc (W, X, Y, Z) options operator, arrow; beta*(W-X-Y-Z)*X-mu*X-X end proc, proc (W, X, Y, Z) options operator, arrow; alpha*gamma*X-(mu+delta1)*Y end proc, proc (W, X, Y, Z) options operator, arrow; (-alpha*gamma+1)*X-(mu+delta2)*Z end proc

  

diff(W(t), t) = A-mu*W(t)-delta1*Y(t)-delta2*Z(t), diff(X(t), t) = beta*(W(t)-X(t)-Y(t)-Z(t))*X(t)-mu*X(t)-X(t), diff(Y(t), t) = alpha*gamma*X(t)-(mu+delta1)*Y(t), diff(Z(t), t) = (-alpha*gamma+1)*X(t)-(mu+delta2)*Z(t)

 (1)
> 

params := {A = 1, alpha = 0.4e-2/gamma, beta = 0.9e-3, delta1 = 0.139e-2, delta2 = 0.134e-2, mu = 0.132e-2}; ics := {W(0) = 750, X(0) = .1, Y(0) = .2, Z(0) = .6}, {W(0) = 755, X(0) = .2, Y(0) = .4, Z(0) = .9}, {W(0) = 760, X(0) = .3, Y(0) = .6, Z(0) = 1}, {W(0) = 765, X(0) = .4, Y(0) = .7, Z(0) = 1.4}

{W(0) = 750, X(0) = .1, Y(0) = .2, Z(0) = .6}, {W(0) = 755, X(0) = .2, Y(0) = .4, Z(0) = .9}, {W(0) = 760, X(0) = .3, Y(0) = .6, Z(0) = 1}, {W(0) = 765, X(0) = .4, Y(0) = .7, Z(0) = 1.4}

 (2)
> 

 

> 

for j to numelems([ics]) do sol[j] := dsolve(eval([sys, ics[j][]], params), type = numeric); p[j] := odeplot(sol[j], [t, Y(t)], 0 .. 1800) end do; display([seq(p[j], j = 1 .. numelems([ics]))], color = [red, green, blue, black])

 
> 

``

Download badODE.mw

First 23 24 25 26 27 28 29 Last Page 25 of 207