tomleslie

13876 Reputation

20 Badges

15 years, 173 days

MaplePrimes Activity


These are answers submitted by tomleslie

?DETools:-convertsys()

much more than I expected!

So I am going to try to deall with many/most of the points raised.

OP stated

I might be wrong, but since turning the inequalities into equalities makes them become the equations of the lines, I guessed they could then be used as constraint equations for LagrangeMultipliers.

Emphatically NO.

The process you suggest will allow you to check for extrema along the boundaries of your quadrilateral - since here the constraints can be expressed as equalities. However it will not allow you to check for extrema within the quadrilateral: because within the quadrilateral the constraints cannot be expressed as equalities

OP's second example

For  defined in the range  the candidate points are , the maximum is  (f = 66) and the minimum is (f = -245), but Maple suggests (4273/2650, 177/212) (f = 15.28) as maximum and (-1, 11/2) (f = -77) as minimum.

I agree that the in-built Maple optimisation routines return incorrect answers for this problem. However, using the (free) DirectSearch add-on available from the Maple Applications Centre, this finds both the maximum at [-7,-6] and the minimum at [6,-7]

I have extended my original worksheet to include your second example. This is attached below.

Note that I have uploaded this worksheet including the output, so that you can see the resuts from the DirectSearch package. Note that whilst you can read this worksheet, I'm guessing that you will not be able to execute it - because you probably(?) don't have the DirectSearch package installed

ineqProb2.mw

You are perfectly justified in asking the question - why does DirectSearch work when the inbuilt Maple routines do not?! Such a question would (probably) start a lengthy and ultimately fruitless discussion. For non-convex functions, there is no 100% guaranteed method of finding a global maximum(minimum), even over a restricted region. It is (probably) best to accept that some optimisers are better than others - inconvenient but true!!

Carl Love's response

He is absolutely correct (he usually is!) that there are two possibilities for the location of extrema of a function over a restricted region. These are

  1. an inflection point within the restricted region: Generally such points can be found as zeros of the gradient. In other words, for any function f(x,y), one needs to find inflection points satisfying both diff( f(x,y), x)=0 and diff( f(x,y),y)=0. This is do-able. Having generated all such points (and depending on the supplied function, there may be none, a finite number, or an infinte numer of these). One then need to filter these to determine which are within the region of interest. Then one need to apply a second differential test to sort these into potential maxima, minima or saddle-points
  2. However the maximum (minimum) of the required function over a defined region need not be an inflection point, if it occurs on the boundary of the region. So one laso need to chck the function values along the region boundaries. Note that the region boundaries can be given by an equality constraint, so this part of the problem can be addressed using Lagrange multipliers
  3. Having obtained  all of the potential maxima (minima) in the interior of the region, and also along the region boundaries, simple function evaluation can then find the 'global' maximum (minimum)
  4. I considered doing this (honest), but it is quite a lot of work, and for your original problem, use of the in-built optimiser was much quicker/simpler. Even for your second example, the use of the DirectSearch optimiser was much quicker/simpler. However I agrre with Carl that optimisation approaches produce 'only' numerical answers, whose correctness cannot be guaranteed - but this 'analytic' approach *ought* to produce a 'symbolic' solution.
  5. I can however envisage situations where thsi symbolic approach might/would fail - probably not applicable to your "relatively simple" objective function. But consider, for example, the case of a non-continuous (hence non-differentiable) objective function - or possibly a periodic objective function where simple differentiation might produce an infinte number of possible extrema.....
  6. So I'm pretty sure that this approach could be made to work for relatively simple objective functions. And I'm also pretty sure that I could come up with a function which would 'break' this solution method

Kitonum's Response

This is essentially an observation on the efficacy of optimisers. No optimisation startegy I am aware of is guaranteed ot find a global optimum for a non-convex function. However if one provides a 'suitable' starting point then one can 'force' a solution - basically by finding the local optimum 'closest' to the initial point. However there are many situations where determining a sutable initial point is not really possible. Consider - if you could guarantee a method of providing a suitable 'initial point' which would guarantee the calculation of a global optimum, then you have solved the general global optimization problem - and I'm pretty sure no-one has achieved this

Not sure that you can use Lagrange Multipliers when you have inequality constraints.

But it is easy to use the optimisation package to find a numerical solution, as in

  constr:={ y>=(1/21)*x-283/42,
                   x<=(3/29)*y+329/58,
                   y<=(1/9)*x+131/18,
                   x>=-(1/9)*y-113/18
               }:
#
# plot region defined by constraints
#
   plots:-inequal( constr,
                             x=-10..10,
                             y=-10..10
                          );
#
# Find maximum of function subject to constraints
#
   Optimization:-Maximize( 2*x^2+4*x*y+2*x-2*y^2+y+4,
                                             constr
                                          );

 

  1. Make sure that the file is not open in another application (ie Excel?)
  2. Make sure that Maple knows exactly where the file is located by supplying a full path name: ie something like

ExcelTools:-Import("C:/Users/TomLeslie/myMaple/test2.xlsx", "Sheet1", "A1:C3");

This is one of these cases where it is actually simpler to write completely new code from the description of your problem, rather than debug the code you supplied. My attempt is attached

oddProc.mw

Only issue  I have with this code is the definition of the list you want as the second element of the returned Array() when repetitions occur. This is not clear from your description - but if my guess is incorrect then I could modify the attached code to return any list you want (in more or less any order)

Tricky one

The error in yur worksheet are a combination of several factors

The first problem, with the assign() command, is almost certainly because you re-executed this command without a clean restart. The assign() command evaluates both of its arguments. In your specific case, this will work properly provided that the first argument (ie 'X[k]') has not already been assigned to a numerical value- but if it has, then the command will fail. Consider the command sequence

restart;
assign(a,1);
assign(a,2);

The first assign() command will work, becuase 'a' is unassigned. The second assign() command will fail, becuase in this command, 'a' will be evaluate to 1 (from the first command) and then you try to make the assignment 1:=2. Sorry but this is never going to happen! Apart from anything else, Maple names (ie things you can assign to) have to start with a non-numeric

If you really want keep the assign(0 statement as written, then ensure that any time you perform a re-execution of the worksheet, do so from the start - including the restart command, to ensure that all names are unassigned.

Your second problem arises from the command

Convert([seq(Distance.....`)

'Convert' should be 'convert' and 'Distance' will be unevaluated becaase you haven't loaded the required Student[Precalculus] package, of which 'Distance' is a method.

In the attached I have fixed all of the above errors (and possibly a few more which I have forgotten about!)

geodFix.mw

However you have a bigger and more fundamental problem.

The method in RJLopez' original worksheet (which I have reproduced), works very well when the surface over which geodesics have to be calculated can be expressed in the form

unique z-coordinate=f(x,y)=some function of x-coordinate and y-coordinate

Notice, however, that for a sphere there is no unique z-value for each (x,y) pair. There are two such z-values. Under this circumstance, the algorithm will fail. It *ought* to work for a cone (axis parallel to z), but not for a cylinder (with any axis, because this will have somewhere betwee 2 and infinity z-values, for each (x-y) pair).

The problem of calculating geodesics on a sphere (or cylinder) is not difficult - in fact it is much simpler than the case considered in RJLopez' original worksheet. But the algorithm in the latter cannot be applied to surfaces which have multiple z-values for given (x-y) pairs

 

Download/install the DirectSearch package from

http://www.maplesoft.com/applications/view.aspx?SID=101333

Replace your fsolve() command with

DirectSearch[SolveEquations]([entries(EQQ, nolist)], evaluationlimit=1000000)

After about five minutes (on my aging machine) this will return the answer

[0.1921937446143878480084068181265523577652e-2,
 Vector[column](4, [` 1 .. 36 `*Vector[column],
`Data Type: `*anything,
`Storage: `*rectangular,
`Order: `*Fortran_order]),
[AA[1] = .5168920149275052365972754691683655791422,
AA[2] = 0.5100492565996477252826303801823045271666e-1,
AA[3] = -1.158385090558286218804863895556307296586,
AA[4] = .2837803429654482579543478405366628528333,
AA[5] = 0.6397053215209891666517565145617507182090e-1,
AA[6] = 0.1673392802413355420838191730032496945114e-2,
AA[7] = 1.172248947532987136781629960218909444346,
AA[8] = 0.1230694399937193179788176515101367092991e-1,
AA[9] = -3.104692044771448609087019829812828207207,
AA[10] = .1197607490553816240419646122284835815162,
AA[11] = 0.5147642793817588352199585647008031691750e-1,
AA[12] = 0.6536847665431726689386446404460938833486e-2,
AA[13] = .9701114867162291146972436326178433865351,
AA[14] = 0.1230808678759375090184337192592145944061e-1,
AA[15] = -2.562299401678170502388271917932774627788,
AA[16] = .1197616195743942198905331008246897752656,
AA[17] = 0.5147602896928951020492183251662888176406e-1,
AA[18] = 0.6536840335570423709881333883288829634921e-2,
AA[19] = 1.172147873904211810098287793264506643493,
AA[20] = 0.1222379533592305040268583423079962871502e-1,
AA[21] = -3.104746264844391836304232194300071559982,
AA[22] = .1191871284123785751064343877594609923404,
AA[23] = 0.5126197690004876021034942475180200887251e-1,
AA[24] = 0.6367317266868154020034897132110644127898e-2,
AA[25] = .9700104194646210892588490090699183639792,
AA[26] = 0.1222494446215037663371613117651060499250e-1,
AA[27] = -2.562353668835540388315280742862321917827,
AA[28] = .1191880315089809310831574377384966398293,
AA[29] = 0.5126154541158110772951899195632797975374e-1,
AA[30] = 0.6367342728951599377764301128147395737756e-2,
AA[31] = -.3477888526305781039106651445931594305721,
AA[32] = 0.5101234853088738537914792274773168084661e-1,
AA[33] = 1.161784920127796330121514485801180941873,
AA[34] = .2837816172198358297004446938437821760649,
AA[35] = 0.6396973710694829468788724620135892531260e-1,
AA[36] = 0.1673529439981824938314270905104991298476e-2],
121537]

Use one of the midpoint methods for dsolve(). The following worked for me

numsol1 := dsolve( { BCSforNum1,
                                    BCSforNum2,
                                    ODEforNum1,
                                    ODEforNum2
                                 },
                                numeric,
                                output = listprocedure,
                                method = bvp[midrich]
                             );

  1. Your ' academic books' are rubbish
  2. You can't type accurately

Because

  1. There is no command 'polygons' in Maple
  2. plot3d() would not accept polygons() as an argument (see (1) above) - in fact plot3d() would not accept polygon() as an argument eithe
  3. A single polygon is described by a list of vertices. Each vertex is a list of coordinates. Hence each polygon is a list of lists. And if you have multiple polygons the you will have a list of lists of lists - in other words three levels of nested square brackets - your command only has two

The following shows a couple of ways to plot a list of polygons - there are probably many others


 

restart;
with(plottools):
with(plots):

display
( polygon
  ( [ [[0, 0, 0], [1, 0, 0], [1, 0, 0], [0, 1, 0]],
      [[0, 0, 0], [0, 1, 0], [0, 1, 1], [0, 0, 1]],
      [[1, 0, 0], [1, 1, 0], [1, 1, 1], [1, 0, 1]],
      [[0, 0, 0], [1, 0, 0], [1, 0, 1], [0, 0, 1]],
      [[0, 1, 0], [1, 1, 0], [1, 1, 1], [0, 1, 1]],
      [[0, 0, 1], [1, 0, 1], [1, 1, 1], [0, 1, 1]]
    ]
  ),
  color=[red, green, blue, orange, yellow, black]
);

polygonplot3d
( [ [[0, 0, 0], [1, 0, 0], [1, 0, 0], [0, 1, 0]],
    [[0, 0, 0], [0, 1, 0], [0, 1, 1], [0, 0, 1]],
    [[1, 0, 0], [1, 1, 0], [1, 1, 1], [1, 0, 1]],
    [[0, 0, 0], [1, 0, 0], [1, 0, 1], [0, 0, 1]],
    [[0, 1, 0], [1, 1, 0], [1, 1, 1], [0, 1, 1]],
    [[0, 0, 1], [1, 0, 1], [1, 1, 1], [0, 1, 1]]
  ],
   color=[red, green, blue, orange, yellow, black]
);

 


 

Download polyplot.mw

Making a few guesses about what you latest worksheet intended, try the following.

At least it provides all function values for all specified time values


 

restart;

#
# Define a whole bunch of parameters
#
  c:= 0.4e-1:    f:= .4:        beta:= .2:  epsilon:= .8:
  theta[1]:= .1: theta[2]:= .3: alpha:= .9: rho:= .7:
  eta := .99:    delta := .3:   a:= 0.4e-1: phi:= 1:
  lambda:= 0.96e-1:

#
# Define OP's ODEs
#
  odesys:= diff(e(t), t)
           =
           lambda*s(t)-(delta+a+epsilon)*e(t),

           diff(i(t), t)
           =
           delta*e(t)-(eta+a+epsilon)*i(t),

           diff(r(t), t)
           =
           eta*i(t)+f*theta[2]*v(t)-(a+epsilon)*r(t),

           diff(s(t), t)
           =
           (1-phi)*epsilon+(1-rho)*a+(1-f)*alpha*v(t)-(lambda+theta[1]+a+epsilon)*s(t),

           diff(v(t), t)
           =
           phi*epsilon+rho*a+theta[1]*s(t)-((1-f)*alpha+f*theta[2]+a+epsilon)*v(t):
#
# Specify the initial conditions
#
  ICS :=  e(0) = .24, i(0) = .17, r(0) = .13,
          s(0) = 0.6e-1, v(0) = .4:

NULL

#
# Solve the ODE system with given initial conditions
# 'sols' contains an 'exact' solution in rationals,
# but since the OP wants solutions at (float) values
# of the independent variable 't', I'm guessing that
# the 'floating point' versions in fsols will probably
# be adequate
#
  sols := dsolve({odesys, ICS}):
  fsols := evalf~(sols);

#
# Generate numeric solutions at time values specified
# by the OP.
#
# NB This data could be reorganised/presented
# in so many different ways!!
#
  seq( eval~( fsols, t=tVal), tVal=0.1..1.5,0.2);

 

NULL


 

Download odeProb.mw

You have a system of four first order ODEs. These ODEs contain 11 unknown parameters , ie

a,alpha,delta,eta,f,lambda,varphi,rho,varepsilon,theta[1],theta[2],

Solcing these four ODEs, in the absence of boundary/initial conditions,will introduce a further four undetermined constants. Hence the solution which Maple finds is a bit on the *long* side!

The existence of the completely unknown and impossible-to-determine function 'i(t)' probably isn't helping


 

NULL

restart; odesys := {diff(e(t), t) = lambda*s(t)-(delta+a+epsilon)*e(t), diff(r(t), t) = eta*i(t)+v(t)*f*theta[2]-(a+epsilon)*r(t), diff(s(t), t) = (1-phi)*epsilon+(1-rho)*a+(1-f)*alpha*v(t)-(lambda+theta[1]+a+epsilon)*s(t), diff(v(t), t) = phi*epsilon+rho*a+theta[1]*s(t)-((1-f)*alpha+f*theta[2]+a+epsilon)*v(t)}; `union`(`~`[indets](odesys)[]); dsolve(odesys)

 

``


 

Download odeProb.mw

  1. variable 'xf' is used but never defined, I have asssigned xf:=xv
  2. 'Delta' and 'delta' are two differnet variables. Maple is case-sensitive. I have used 'delta' throughout
  3. I tidied up the plot commands

The attached now works

.myParab.mw

As in

foo:=proc(x)
local y,z:
x^2;
end proc;
maplemint(foo);

 

Wrap the  ODE probelm and subsequent integral evaluation within a procedure, which accepts initial values as arguments. It is then possible to use this procedure to return specific values of the integral for specific values of the initial conditions. It is also possible to supply this procedure as the first argument of a plot3d() and hence obtain a plot of the required integral as R(0) and mu(0) vary


 

  restart:

#
# Define procedure which accepts IC values as arguments
# and returns integral value
#
  getInteg:=proc(x,y)
                 local A0:= 10^(-3),
                       a:= 10^5,
                       cond:= [ R(0) = x, mu(0) = y, (D(mu))(0) = 0],
                       sys := [ diff(R(theta), theta) = A0*exp(2*mu(theta))*sin(theta)/(2*a),
                                R(theta) = 2*exp(-2*mu(theta))*(1-(diff(mu(theta), theta$2))-cot(theta)*(diff(mu(theta), theta)))],
                       F := dsolve({cond[], sys[]}, numeric, output = listprocedure):
               return  evalf(Int(sin(theta)*exp(2*rhs(F[3])(theta))*(rhs(F[2])(theta)+rhs(F[2])(theta)^2+1), theta = 0 .. Pi));
         end proc:

#
# Check procude for the OP's original case
# ie R(0)=1.0e-05 and mu(0)=100
#
  getInteg(1.0e-05, ln(100));

#
# Plot for ranges of R[0] and mu[0] values
# Reduced the number of grid points just to
# just to speed things up - Doesn't look
# like the integral value is doing anything
# particularly exciting anyway. NB the (default)
# 50*50 grid takes ~20mins to produce
#
  CodeTools:-Usage( plot3d
                    ( getInteg,
                      0.5e-05..1.5e-05,
                      ln(50)..ln(100),
                      style=surface,
                      grid=[10,10]
                    )
                  );

 

 


 

Download doInt.mw


 

  restart:
#
# Generate closed form expression for recursion
#
  r:= unapply
      ( rsolve
        ( { F(m)=2*F(m-1)+4*F(m-2),
            F(0)=2,
            F(1)=4
          },
          F
        ),
        m
      );
  seq( round(r(j)),
       j=0..20
      );
#
# Or if you want "pretty" output
#
  seq( printf( "\t\t%4d %d\n",
               j, round(r(j))
             ),
       j=0..20
     );

#######################################################
# Or using recursive procedure with remember table
#
  restart:
  F := proc(n)
            F(n) := 2*F(n-1)+4*F(n-2);
       end proc:
  F(0):=2:
  F(1):=4:
  seq( F(j), j=0..20);
#
# Or if you want "pretty" output
#
  seq( printf( "\t\t%4d %d\n",
               j, F(j)
             ),
       j=0..20
     );

 


 

Download recur.mw

First 165 166 167 168 169 170 171 Last Page 167 of 207