tomleslie

13876 Reputation

20 Badges

15 years, 175 days

MaplePrimes Activity


These are answers submitted by tomleslie

just because I like alternatives: Get the remainder when dividing the number by 10, as in

c1:=3: c2:=32:c3:=1.234567:c4:=1.274e-23:
irem(op(1,c1),10);
irem(op(1,c2),10);
irem(op(1,c3),10);
irem(op(1,c4),10);
                              
3
                               2
                               7
                               4

I'm pretty(?) sure this will work for all number formats

First of all I should have noticed that the problem of 'r' being unknown had arisen as a result of using the same name for indexed and unindexed variables - sorry about that.

You never use the fact that these names are indexed, so I I have changed all such names, to have 'inert' suffixes: this approach is just 'safer'

The dsolve() command then 'errored' with the message "initial Newton iteration is not converging". I overcame this by multiplying various terms in your ODE system with a variable 'omega', which is used as a continuation parameter by the dsolve() command. This has no effect on the final solution (other than permitting it to be obtained).

The dsolve() command then complained that it was unable to obtain sufficient accuracy, so I increased the 'maxmesh'.

The attached shows the solutions obtained. I can't display these inline here because, I used an ArrayPlot which this site doesn't like - so link only

odeProb.mw

 

your differential equations are first order in s(t), x(t), y(t), lambda[1](t), lambda2(t) and lambda[3](t) - that requires six boundary/initial condtions. However the ode system also contains the unknown variable 'r', which requires another boundary condition (ie some way to evaluate 'r'), making a total of seven.

You only have 6 boundary/initial conditions, with no way of determining the variable 'r'. Hence the system cannot be solved

restart;
expr:=diff(c(x,y),x,x)*diff(a(x,y),x)*c(x,y)*x^6 - 2*diff(b(x,y),x)^3*x^6 + 2*diff(c(x,y),x)^2*c(x,y)*x^5 + a(x,y)^5 = 0;
member(a(x,y), indets(expr));
member(diff(a(x,y),x), indets(expr));
member(diff(a(x,y),y), indets(expr));

which will return

true
true
false

 

 

 

As supplied the definition of your epression contains the "interesting" variable named 'RootOf8'. I have no idea how it got there, since it is a 'name', not RootOf(8). I have taken the liberty of replacing the name RootOf8 with sqrt(8), then used allvalues() to compute all branches of the resulting function. Two of these can be plotted directly, and two seem to be irredeemably complex, so will need special treatment for plotting: if you want to plot the complex curves then please state how - eg absolute value, real value, 3D plot with real/complex components whatever, and I'm sure I can figure something out

The two real branches of you expression are plotted in the attached

sigma2 := RootOf(-90526382422649463214540800000000*Pi^8-18587959930253464168320000000000*Pi^6-5863377073505044924800000000000*Pi^4+43980465111040000000000000000*sqrt(3)*Pi^25*sqrt(32*Pi^2+2)*sigma+21990232555520000000000000000*sqrt(3)*Pi^23*sqrt(32*Pi^2+2)*sigma-98268851732480000000000000000*sqrt(3)*Pi^21*sqrt(32*Pi^2+2)*sigma-44495861186560000000000000000*sqrt(3)*Pi^19*sqrt(32*Pi^2+2)*sigma+82188225740800000000000000000*sqrt(3)*Pi^17*sqrt(32*Pi^2+2)*sigma+33095407370240000000000000000*sqrt(3)*Pi^15*sqrt(32*Pi^2+2)*sigma-30136000839680000000000000000*sqrt(3)*Pi^13*sqrt(32*Pi^2+2)*sigma-10618895073280000000000000000*sqrt(3)*Pi^11*sqrt(32*Pi^2+2)*sigma+3822293002240000000000000000*sqrt(3)*Pi^9*sqrt(32*Pi^2+2)*sigma+1210118016000000000000000000*sqrt(3)*Pi^7*sqrt(32*Pi^2+2)*sigma+118805400000000000000000000*sqrt(3)*Pi^5*sqrt(32*Pi^2+2)*sigma+5028750000000000000000000*sqrt(3)*Pi^3*sqrt(32*Pi^2+2)*sigma+79101562500000000000000*sqrt(3)*sigma*Pi*sqrt(32*Pi^2+2)+111484894360500000*Pi^2*20^RootOf8+1765920726670320000*Pi^4*20^RootOf8-569534208772147200*Pi^6*20^RootOf8-4505569481375428608*Pi^8*20^RootOf8+972005049637797888*Pi^10*20^RootOf8+5143616921914048512*Pi^12*20^RootOf8-554194415829123072*Pi^14*20^RootOf8-2216777663316492288*Pi^16*20^RootOf8+(-9231519020818020433920000000000*Pi^22+195541371952408496701440000000000*Pi^20+89300299589267320995840000000000*Pi^18-333503605675043554590720000000000*Pi^16-115500365322956203622400000000000*Pi^14+204706142659640339988480000000000*Pi^12+55783620627641021399040000000000*Pi^10-43454880575740151285760000000000*Pi^8-9286786763553830541120000000000*Pi^6-635208422610519981000000000000*Pi^4-16054449064166199375000000000*Pi^2-85686765999732421875000000)*_Z+(1683627180032000000000000000000*Pi^28+947040288768000000000000000000*Pi^26-243897798836910985052160000000000*Pi^24-105849518880314282213376000000000*Pi^22+543806205557386676011008000000000*Pi^20+206745517628405562998784000000000*Pi^18-493535946568048375234560000000000*Pi^16-161556685841710476165120000000000*Pi^14+209521703041307302907904000000000*Pi^12+57932333046211895115008000000000*Pi^10-32606166808014116503296000000000*Pi^8-7574931806403147431400000000000*Pi^6-916854325001083153125000000000*Pi^4-60848666758777034179687500000*Pi^2-1531121744500488281250000000)*_Z^2+(14538675656595603456000000000000*Pi^20+6360670599760576512000000000000*Pi^18-24363640065154351104000000000000*Pi^16-9459367828326973440000000000000*Pi^14+10040437028153917440000000000000*Pi^12+3693930616897744896000000000000*Pi^10+1609933205706216192000000000000*Pi^8+58674582771546096000000000000*Pi^6-1202653471578517170000000000000*Pi^4-149668239567146343750000000000*Pi^2-4663745768352832031250000000)*_Z^3+(-8723205391669003498291200000000*Pi^24-4361602695834501749145600000000*Pi^22+19490912047010429691494400000000*Pi^20+8825430454852624633036800000000*Pi^18-16301436833461042151424000000000*Pi^16-6564233354119088386867200000000*Pi^14+5977256592087501137510400000000*Pi^12+2106180608207148770918400000000*Pi^10-758124018049754123827200000000*Pi^8-240018107472837924480000000000*Pi^6-23564187036740637000000000000*Pi^4-997415989180706250000000000*Pi^2-15689219628471679687500000)*_Z^4-13647882752248245117187500000-261292721157421875*20^RootOf8-535230827832343213125000000000*Pi^2+305811336261213249011712000000000*Pi^12+79115470702645314657484800000000*Pi^10-239241111641945951698944000000000*Pi^16-79895480796476508576153600000000*Pi^14+7986315188014109687808000000000*Pi^22-60346149989113268482867200000000*Pi^20-18258684357505568263372800000000*Pi^18+14855623787650488886886400000000*Pi^24)

indets(sigma2, 'name');
fred:=allvalues(subs(RootOf8=sqrt(8), sigma2)):

{Pi, RootOf8, sigma}

(1)

F := plot([fred], sigma = -10 .. 10, color = [red, green, blue, orange]);

Warning, unable to evaluate 2 of the 4 functions to numeric values in the region; complex values were detected

 

 

NULL


 

Download getPlots.mw

  1. Maple has a probe tool for plots which appears on the Plot2D and Plot3D Toolbars (as well as on the plot context menus). However this only works for 2D plots :-(
  2. If you are plotting an algebraic exprssion then you can obtain the expression value at any point by using the eval() command, whihc will be something like evalf(eval(expression, [x=xvalue,y=yvalue])). One can also compute maximum, minimum and their locations symbolically using evalf(maximize(expression, location=true))[2][];
  3. If you do not have an algebraic expression, eg you are producing a surface plot from data (or similar), then things get a bit trickier. It is possible to retrieve the data from the plot structure and use this to generate an interpolating function, using the Interpolation() package. This seems a bit long-winded becuase it would probably be better to run the interpolation on the data which was used to produce the plot in the first place. One can then supply any coordinate values to this interpolating function to generate the associated (interpolated) function value

Thes alternatives are illustrated in the attached, where I have used a 'test' 2-D function, since the actual function was not supplied


 

  restart;
#
# Define/plot a test function
#
  f:=x*exp(-x^2-y^2);
  p1:=plot3d(f, x=-2..2, y=-2..2, grid=[100,100]);

x*exp(-x^2-y^2)

 

 

#
# Get max/min from the function definition
#
  evalf(maximize(f, location=true))[2][];
  evalf(minimize(f, location=true))[2][];
#
# Define a utility which will return the
# function value from its definition
# for any specified cooordinate values
#
  getVal:=(p,q)->evalf(eval(f, [x=p,y=q])):
  getVal(1.5,1.5);

[{x = .7071067810, y = 0.}, .4288819424]

 

[{x = -.7071067810, y = 0.}, -.4288819424]

 

0.1666349481e-1

(1)

#
# Define a procedure which  numerically interpolates
# the plot data to approximate the original function
#
  fInterp:= proc(p::function)
                 uses plottools, Interpolation:
                 local p1d:= getdata(p),
                       xos:= op(1,p1d[2][1]),
                       xstep:= (op(2,p1d[2][1])-op(1,p1d[2][1]))
                               /
                               (op([2,1,2], p1d[3][])-1),
                       yos:= op(1,p1d[2][2]),
                       ystep:= (op(2,p1d[2][2])-op(1,p1d[2][2]))
                               /
                               (op([2,2,2], p1d[3][])-1),
                       coords:= Array( 1..op([2,1,2],p1d[3][])*op([2,2,2],p1d[3][]),
                                       1..2,
                                       [ seq
                                         ( seq
                                           ( [ xos+xstep*(j-1), yos+ystep*(k-1)],
                                             j=op([2,1],p1d[3][])
                                           ),
                                           k=op([2,2],p1d[3][])
                                         )
                                       ]
                                     ),
                       data:= Array( 1..op([2,1,2],p1d[3][])*op([2,2,2],p1d[3][]),
                                     [ seq
                                       ( seq
                                         ( p1d[3][][i,j],
                                           i=op([2,1],p1d[3][])
                                         ),
                                         j=op([2,2],p1d[3][])
                                       )
                                     ]
                                   );
                return Interpolate(coords, data):
            end proc:
  fint:= fInterp(p1):

#
# generate a few value from the interpolation and
# the original expression for comparison purposes
#
  fint(1.5,1.5);
  fint(0,0);
  eval(f, [x=1.5,y=1.5]);
  eval(f, [x=0,y=0]);

HFloat(0.01669542421616841)

 

HFloat(0.0)

 

0.1666349481e-1

 

0

(2)

#
# Use DirectSearch to find global/local maxima/minima
# from the interpolation function
#
  with(DirectSearch):
  GlobalSearch( fint(x,y),
                [x=-2..2, y=-2..2]
              );
  GlobalSearch( fint(x,y),
                [x=-2..2, y=-2..2],
                maximize=true
              );

Matrix(5, 3, {(1, 1) = -.42870694049523717, (1, 2) = [x = -.707070717258705, y = 0.2020198891840375e-1], (1, 3) = 117, (2, 1) = -.42870693884298994, (2, 2) = [x = -.7070706480684299, y = -0.20201993982392433e-1], (2, 3) = 117, (3, 1) = -.41826581753161496, (3, 2) = [x = -.5923523377305716, y = -0.5058121049210635e-1], (3, 3) = 92, (4, 1) = 0.6709252837965911e-3, (4, 2) = [x = 1.9999999924168164, y = 1.9999999983516001], (4, 3) = 62, (5, 1) = 0.6709253119671453e-3, (5, 2) = [x = 1.999999978704228, y = -1.9999999997442524], (5, 3) = 46})

 

Array(%id = 18446744075418877214)

(3)

 


 

Download plotInterp.mw

help(DirectSearch);

or

?DirectSearch;

both work for me. You can then navigate to any subtopic like GlobalSearch(). Interestingly(?), if I use

?DirectSearch:-GlobalSearch;

this takes me to the 'wrong' DirectSearch() help page. I get the first  (alphabetically) of the DirectSearch sub-topics, which happens to be DirectSearch[BoundedObjective].

If you are not getting any DirectSearch help pages at all, an obvious thing to check is whether you have actually installed the 'help' for this package. Also ensure that the help database is in tha appropriate format for your Maple version. Prior to Maple 18, help files were in .hdb format: from Maple 18 onwards these files should be in .help format.

The first paragraph of the DirectSearch:-GlobalSearch() help page is as shown below - my highlighting

The GlobalSearch command attempts to find numerically all local and global minimums (maximums) of a real-valued nonlinear multivariate function f with (without) constraints constr. The GlobalSearch is multimodal derivative-free multistart direct searching optimization method, i.e. it does not require the objective function f(x1,x2,…,xn) and constraints to be differentiable and continuous.

The word "attempts" occurs because there is no guaranteed way to find a global optimium for a non-convex problem.

The word "multistart" occurs because (like almost all other "global" optimisers), the basic strategy is to run mutiple optimizations starting from a different intial points

the VectorCalculus() commands SurfaceInt() and Flux() to make this calculation a lot simpler, as in the attached

restart;
with(VectorCalculus):
SetCoordinates(cartesian[x, y, z]):
u1:=x: u2:=1:u3:=z: s:=x^2+y^2+z^2-2:
u:= VectorField([u1, u2, u3]):
n:= Normalize(Gradient(s)):
SurfaceInt( u.n, [x,y,z] = Sphere( <0,0,0>, sqrt(2) ) );
Flux(u, Sphere(<0,0,0>, sqrt(2), outward));

(16/3)*2^(1/2)*Pi

 

(16/3)*2^(1/2)*Pi

(1)

 


 

Download gaussTheo.mw

so you can't integrate that term. You can integrate the others by expanding the integral as in
 

restart;
u[1](r,z):=(C[o]^2*exp(2*lambda*z)-D[o]^2*exp(-2*lambda*z))+(lambda/2)*(C[o]*exp(lambda*z)-D[o]*exp(-lambda*z))*(1-r^4)+A[1](z)*r^2;
M1:=expand(int(u[1](r,z), z));

C[o]^2*exp(2*lambda*z)-D[o]^2*exp(-2*lambda*z)+(1/2)*lambda*(C[o]*exp(lambda*z)-D[o]*exp(-lambda*z))*(-r^4+1)+A[1](z)*r^2

 

(1/2)*C[o]^2*(exp(lambda*z))^2/lambda+(1/2)*D[o]^2/(lambda*(exp(lambda*z))^2)-(1/2)*r^4*C[o]*exp(lambda*z)-(1/2)*r^4*D[o]/exp(lambda*z)+(1/2)*C[o]*exp(lambda*z)+(1/2)*D[o]/exp(lambda*z)+r^2*(int(A[1](z), z))

(1)

 


 

Download intProb.mw

bu this will work,

restart;
GST := AE+AI-ANB-AR-CP;
subs([seq( j=cat(j, [i]), j in indets(GST))], GST);

 

of the code I previously posted will handle variables in steps of 0.1. The brute forceapproach no loner looks quite so good becuase ther are 100*50*10 evaluations to perform, so time-wise it is starting to hurt. But it still works! But since you now have 7599 possible solutions, I doubt if there any way to find them (all) faster!

See the attached - which for some reason still will not display here

Download getsols2.mw

 

like this, where the function is trivial to evaluate and there are only ~100 combinations, to check, a brute force approach isn't necessarily(?) bad. See the attached (whose contents for some reason won't display here!)
 

Download getsols.mw

 

you have to remember to distinguish between an equation and an assignment, As in


 

  restart;
  expr:=com2(e,e,f)=eh+he;
  subs({eh=<1,0>, he=<0,1>}, expr);

expr := com2(e, e, f) = eh+he

 

com2(e, e, f) = Vector[column](%id = 18446744074395821470)+Vector[column](%id = 18446744074395821590)

(1)

#
# or
#
  restart;
  com2(e,e,f)=eh+he;
  subs({eh=<1,0>, he=<0,1>}, %);

com2(e, e, f) = eh+he

 

com2(e, e, f) = Vector[column](%id = 18446744074395821230)+Vector[column](%id = 18446744074395821350)

(2)

 


 

Download dosubs.mw

but the attached fulfils your request

restart:
with(plots):
data := <.150, .175, .150, .175, .200, .300, .250, .250, .350, .400, .450, .550, .550>:
fit:= minimize
      ( add
        ( (data[n+1]-a*data[n]-b)^2,
           n=1..numelems(data)-1
        ),
        location
      );
p1:= pointplot( data[1..-2],
                data[2..-1],
                symbol=solidcircle,
                symbolsize=20,
                color=red
              ):
p2:= pointplot( data[1..-2],
                eval( a*data+~b,
                      fit[2,1,1]
                    ),
                connect=true
              ):
display( [p1, p2],
         labels=[x[n], x[n+1]],
         labelfont=[times,bold,20]
       );

0.2642696629e-1, {[{a = 1.035955056, b = 0.2314606742e-1}, 0.2642696629e-1]}

 

 

 

NULL

Download oddFit.mw

Basically the process is to expand your 2-D integral into two 2-D integrals, then solve one of them by parts.

See the attached.

  restart;
  I1:=int(diff(f(x, y), x)*g(x, y)+f(x, y)*diff(g(x, y), x), x);
#
# Expand to two integrals. Solve first integral
# by parts. Note that user has to 'eyeball' which
# 'part' is u and which is v', although this could
# be automated
#
# Then add the second integral
#
  IntegrationTools:-Parts( op(1,expand(I1)),
                           g(x,y),
                           f(x,y)
                         )
                         + op(2,expand(I1));

int((diff(f(x, y), x))*g(x, y)+f(x, y)*(diff(g(x, y), x)), x)

 

f(x, y)*g(x, y)

(1)

 

Download intParts.mw

First 120 121 122 123 124 125 126 Last Page 122 of 207