acer

32348 Reputation

29 Badges

19 years, 330 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are replies submitted by acer

Is the pattern that `max` and `sort` are still using evalb for relation testing? Maybe it's time to consider promoting `is` to a new job.

For `sort` you can supply your own comparator (as I know you're aware) but there's still an argument to made for having it be more automatic. For `max` it's tougher to do it efficiently.

The generated 3D globe plot consists of POLYGONS with many (approx. 30-40k) elements, CURVES for the borders, and a MESH for the sphere/ocean. Those main substructires can all be manipulated separately, even obtaining a few customizations not easy (or possible) directly from the original command.

globe_fun.zip

@Preben Alsholm 

I'm not sure I understand why evalf/Int isn't succeeding (and quickly) by splitting the integral at x=1 the point of discontinuity in the given range x=0..5.

restart;
g:=sqrt(1/(1-x^2)):
CodeTools:-Usage(evalf(int(g, x = 0 .. 5)));
memory used=52.37MiB, alloc change=34.00MiB, cpu time=502.00ms, real time=502.00ms, gc time=39.53ms
                  1.570796327 + 2.292431670 I
restart;
g:=sqrt(1/(1-x^2)):
CodeTools:-Usage(evalf(IntegrationTools:-Split(Int(g,x=0..5),
                          select(x->x>=0 and x<=5,[op(discont(g,x))]))));
memory used=20.81MiB, alloc change=67.01MiB, cpu time=216.00ms, real time=216.00ms, gc time=15.67ms
                  1.570796327 + 2.292431670 I
restart;
g:=sqrt(1/(1-x^2)):
CodeTools:-Usage(evalf(Int(Re(g), x = 0 .. 5))
                 +evalf(Int(Im(g), x = 0 .. 5))*I);
memory used=6.15MiB, alloc change=33.00MiB, cpu time=67.00ms, real time=66.00ms, gc time=5.67ms
                  1.570796327 + 2.292431670 I

Everything has the same sign, approaching the discontinuity at x=1 from both sides.

plot([Re,Im](g),x=0..5,view=0..5,thickness=4,filled,color=[green,red]);

@Samir Khan Yes. I was supposing that the country data or its processing might confused about that disputed region. (I wrote lake in a dry attempt at humour)

@Samir Khan Is that a large lake along the border of Western Sahara and Mauritania?!

@assma My first suggestion (if diff can handle what you call the "known functions" but which you resist showing us in full) is,

w2:=unapply(diff(v2(x,y),x)-diff(u2(x,y),y), [x,y]):

For some strange reason you changed that to the following, which I did not recommend at all,

w2:=proc(x, y) options operator, arrow; unapply(diff(v2(x, y), x)-(diff(u2(x, y), y))) end proc;

@assma No, I did not suggest putting a call to unapply within the body of w2. Yes, what you've done does solve your original runtime error, but it is as bad for performance as calling diff within w2. That's exactly what I suggested not to do.

@itsme The over-conservativeness lies in LinearAlgebra:-Eigenvalues, but I recall that it's been mosty like that since Maple 6 (yr. 2000) when it was introduced. On the one hand mixed floats and symbolics has always been dodgy (and not just in Maple). On the other hand there's always been a struggle between pragmatism ("just give me an answer") and restricting functionality to clearly-defined domains for which correct behavior is understood (in theory).

@nidojan I think that it's not a very good idea to decrease the accuracy of the intermediate results (using evalf, say) duing the looped calculation. So I've left such rounding-down until the Array/Tabulated pretty-printing stage.

Below, I've moved the call to Jacobian outside the loop, on general efficiency grounds. I've gotten rid of the variable pointt. I changed sol[n] to just sol since the values are already stored in x and y. Also, there will always be an x[n],y[n] pair that is computed after the last instance of the Jacobian. (There's not much point computing the last point's Jacobian and not then definitely not using it.)

I've changed from J^(-1) to a call to LinearSolve, on general grounds of numeric robustness.

Ideally all the code could be put into a easily re-usable procedure. I've left that for you, and the same for fancy stopping criteria such as forward-error or x-increment tolerances.

The colored Table at bottom looks nicer in the actual Maple GUI.

restart

f[1] := proc (x, y) options operator, arrow; 3*x^2-y^2 end proc; f[2] := proc (x, y) options operator, arrow; 3*x*y^2-x^3-1 end proc

proc (x, y) options operator, arrow; 3*x^2-y^2 end proc

proc (x, y) options operator, arrow; 3*x*y^2-x^3-1 end proc

var := x, y

x, y

Jgeneral := Student:-MultivariateCalculus:-Jacobian([f[1](x, y), f[2](x, y)], [var], 'output' = 'matrix')

Matrix(%id = 18446883951302186758)

Jfun := unapply(Jgeneral, [var])

x[0], y[0] := 1.0, 1.0

1.0, 1.0

iter := 5

5

for n from 0 to iter-1 do print(cat('x__', n) = x[n], cat('y__', n) = y[n]); f1val, f2val := f[1](x[n], y[n]), f[2](x[n], y[n]); print(('f[1]')(cat('x__', n), cat('y__', n)) = f1val, ('f[2]')(cat('x__', n), cat('y__', n)) = f2val); J[n] := Jfun(x[n], y[n]); print(cat('J__', n) = J[n]); print(); sol := (Vector(2, {(1) = x[n], (2) = y[n]}))-LinearAlgebra:-LinearSolve(J[n], Vector(2, {(1) = f1val, (2) = f2val})); x[n+1], y[n+1] := sol[1], sol[2] end do

x__0 = 1.0, y__0 = 1.0

f[1](x__0, y__0) = 2.00, f[2](x__0, y__0) = 1.000

`#msub(mi("J"),mi("0"))` = _rtable[18446883951302174950]

NULL

x__1 = .611111111111111, y__1 = .833333333333333

f[1](x__1, y__1) = .425925925925926, f[2](x__1, y__1) = 0.449245541838137e-1

`#msub(mi("J"),mi("1"))` = _rtable[18446883951302172294]

NULL

x__2 = .503659080767515, y__2 = .852494422132976

f[1](x__2, y__2) = 0.342706691508959e-1, f[2](x__2, y__2) = -0.296666581679239e-1

`#msub(mi("J"),mi("2"))` = _rtable[18446883951302160982]

NULL

x__3 = .499964121073170, y__3 = .866045636388624

f[1](x__3, y__3) = -0.142677226374066e-3, f[2](x__3, y__3) = -1.25763153524527*10^(-6)

`#msub(mi("J"),mi("3"))` = _rtable[18446883951302157846]

NULL

x__4 = .500000000014920, y__4 = .866025401817003

f[1](x__4, y__4) = 3.45245754207468*10^(-9), f[2](x__4, y__4) = -5.08916730979081*10^(-9)

J__4 = Matrix(%id = 18446883951302154710)

print(cat('x__', n) = x[n], cat('y__', n) = y[n]); print(('f1')(cat('x__', n), cat('y__', n)) = f[1](x[n], y[n]), ('f[2]')(cat('x__', n), cat('y__', n)) = f[2](x[n], y[n]))

x__5 = HFloat(0.5), y__5 = HFloat(0.8660254037844387)

f1(x__5, y__5) = HFloat(-1.1102230246251565e-16), f[2](x__5, y__5) = HFloat(2.220446049250313e-16)

n := 'n'; Tab := Array([[n, 'x__n', 'y__n', 'J__n', ('f__1')('x__n', 'y__n'), ('f__2')('x__n', 'y__n')], seq([i, x[i], y[i], `if`(i = iter, "", J[i]), f[1](x[i], y[i]), f[2](x[i], y[i])], i = 0 .. iter)]); Tab[1 .. -1, 4] := evalf[6](Tab[1 .. -1, 4]); Tab[1 .. -1, 5 .. 6] := evalf[3](Tab[1 .. -1, 5 .. 6]); Tab

Array(%id = 18446883951302147358)

DocumentTools:-Tabulate(Tab, weights = [1, 6, 6, 5, 3, 3], widthmode = pixels, width = 700, fillcolor = (proc (T, i, jj) options operator, arrow; `if`(i = 1, "#E0E0E0", "#F0FFFF") end proc))

 

n

`#msub(mi("x"),mi("n"))`

`#msub(mi("y"),mi("n"))`

`#msub(mi("J"),mi("n"))`

f__1(x__n, y__n)

f__2(x__n, y__n)

0

1.0

1.0

Matrix(2, 2, {(1, 1) = 6.0, (1, 2) = -2.0, (2, 1) = 0., (2, 2) = 6.00})

2.00

1.00

1

.611111111111111160

.833333333333333370

Matrix(2, 2, {(1, 1) = 3.66667, (1, 2) = -1.66667, (2, 1) = .962963, (2, 2) = 3.05556})

.426

0.449e-1

2

.503659080767514533

.852494422132976326

Matrix(2, 2, {(1, 1) = 3.02195, (1, 2) = -1.70499, (2, 1) = 1.41922, (2, 2) = 2.57620})

0.343e-1

-0.297e-1

3

.499964121073170287

.866045636388624196

Matrix(2, 2, {(1, 1) = 2.99978, (1, 2) = -1.73209, (2, 1) = 1.50021, (2, 2) = 2.59795})

-0.143e-3

-1.26*10^(-6)

4

.500000000014919843

.866025401817003271

Matrix(2, 2, {(1, 1) = 3.00000, (1, 2) = -1.73205, (2, 1) = 1.50000, (2, 2) = 2.59808})

3.45*10^(-9)

-5.09*10^(-9)

5

.500000000000000000

.866025403784438708

 

-1.11*10^(-16)

2.22*10^(-16)

 

 

 


q1nwtnnonlinearsys_2.mw

@Markiyan Hirnyk Yes, but that is advanced. People new to Maple should start with handling and manipulating expressions. Operators and procedures should come after they've understood the basics of expressions. I saw enough in this member's first three Questions to make the judgment call that an introduction to integrating with expressions would be a gentle start.

@Markiyan Hirnyk Yes, I deliberately omitted that one because I didn't want to lead him closer to difficulties with 2D Input and the typeset integral (ie. from palette, or via command completion), or toward the poor habit of defining an operator just for sake of only ever referencing it as a function call like f(x).

@tomleslie Since it is not the job of evalb to test `whether expressions are "equivalent"' then it's not the best tool for your example.

An advantage of `is` over evalb (when it's not hitting a bug) is that `is` can return FAIL in a broader set of (nonnumeric) cases where it cannot ascertain whether the result is true or false. In contrast, when evalb returns false you are left wondering whether some stronger intermediary command would produce the opposite result.

Also, within the context of a conditional test, it can be more awkward in my opinion to deal with evalb returning a relational argument unevaluated than it is to deal with `is` returning FAIL.

restart;

ee := cos(a+b) = cos(a)*cos(b) - sin(a)*sin(b):

evalb( expand(ee) );

true

evalb( combine(ee) );

true

is( ee );

true

 

Download evalb_is.mw

[edit] It might also be of interest to some people that, while simplify neither combines or expands the rhs or lhs of that expression assigned to ee (since neither form is canonically simpler), it will simplify their difference to zero.

simplify( ee );

cos(a+b) = cos(a)*cos(b)-sin(a)*sin(b)

evalb( simplify( ee ) );

false

simplify( (rhs-lhs)( ee ) );

0

evalb( simplify( (rhs-lhs)(ee) ) = 0  );

true

 

@Joe Riel Right. This is not supposed to be the job of evalb.

It is supposed to be the job of is, but unfortunately that can return FAIL for the example in question under the given assumptions (and false without them). I am submitting a bug report against is with this example as details.

note. I get exact zero from each of combine, expand, and simplify.

[edited]

What I would have expected to work is using is like so,

ff := c^10*c^n - c^(n+10):

is( ff  = 0 );
                      true

Possibly related, I notice this difference in behaviour,

restart;
ff := c^10*c^n - c^(n+10):
ee := 2^10*2^n - 2^(n+10):

combine(ee,power), combine(ff,power);                    
                       0, 0

simplify(ee,power), simplify(ff,power);           
                     n    (n + 10)
               1024 2  - 2        , 0

I included such details in my bug report for the original example. I also submitted a report that the following allocated large resources quickly (and required a kill from the OS):
   testeq( -4*2^(a+I*b)+2^(a+I*b+2) )

@ThU See this thread for an earlier discussion. You may scroll through the Comments on my Answer there, as there were refinements to the original for handling units, and addition of context-menu items. Here's a file. Phasors11.zip

I keep meaning to turn it into a special maplecloud/maple2017 package thing. (I started it... help page, move some stuff into ModuleLoad, etc, but then got busy.)  Perhaps it ought to be renamed "DegreePolar". Opinions?

What makes your correction to your nlc2 (to use the appropriate and documented form of procedure for the constraint) work is not that it returns W[1], but rather that it updates the entry of its second argument (the Vector W) inplace. It could return NULL and work as well.

Here is an example like what you first suggested (with strict equality constraints), working in both Maple 2017.3 and Maple 14.01 on 64bit Linux. I used both Operator form and Matrix form. For the Matrix form I also passed in procedures for the objective-gradient and constraint-jacobian. (I'm not sure how much general interest might be drummed up for the failings of Maple 14, seven years and seven major releases later. But I'm willing to help with specific examples.) NLPSolve_operator.mw

Can I ask, is it your goal to try and use NLPSolve as a root-finder?

Since this does not pertain specifically to LPSolve or the rest of this Question thread then it would be more useful all round if you started new discussion in some new thread, in future. You could use the Branch tab at the bottom of a Comment to link the two threads, if you'd like. Thanks.

First 259 260 261 262 263 264 265 Last Page 261 of 592