## 4862 Reputation

6 years, 362 days

## @litun    Are you thinking t...

Are you thinking to something like that?

Run.mw

RunRun.mw

Run.mw contains the code (please note the line writeto(terminal) before the "return' line in order to redirect the results to your active worksheet)/

An example of a double loop using DocumentTools:-RunWorksheet  is given in RunRun.mw

## @Kitonum   [{2}\$2] could be c...

@Kitonum

[{2}\$2] could be cosidered as a little cheating...
... but an astute one :-)

## @Joe Riel    The double quot...

The double quotes are "ugly" indeed,  and your last use of map is definitively simpler.

BTW: before coming here I didn't know about element-wise operators (funny thing is that I thought at this time that Maple suffered from lack when comparred to Matlab or R on this point ...). As I was reading more and more answers using "~" I gradually became a fan of it and began to set aside "map", at least when possible.
I guess the truth lies in a reasonable use of each of them.

## @acer  @Carl Love Thanks to y...

To acer: "Did you mean   unapply('P'(i,1,1,x), x)  so that a procedure P doen't get its call evaluated there?"
No, I didn't mean 'P' but simply P. So to understant the difference between '¨' and P I did this

restart:
P := proc(a, x)
a+x
end proc:

for k from 1 to 1 do
f[k] := unapply(P(k,x), x);
g[k] := unapply('P'(k,x), x);
end do;
f1 := x -> 1 + x
g1 := x -> P(1,x)

f(x);
1+x
g(x);
1+x

Then, even if the final results are the same, using your 'P' instead my naïve P is definitively better, simply because
having g1(x)= P(1,x)  is probably the thing one implicitely looks for when one writes the inloop assignment.

Good weekend to both of you

## I complete my last reply...

@Carl Love

I send you three worksheets withe very slight modifications between them

• Observation 1 : for high values of elements of the matrix Maple can fail to verify if the product of the eigenvalues is equal to the determinant (generator = 0..10^(10^3),  eig1.mw

• Observation 2 : trying to simplify the formal expressions of the eigenvalues and of the determinant may lead to dramatic failures of the test of equality  eig2.mw

• Observation 3 : In  eig3.mwi the simplifications are made on the valued expressions of the product of the eigenvalues and of the determinant. The test of equality seems to perform well even for an upper bound of 10^(10^3),

• Observation 4 : An upper bound of 10^(10^4),​​​​​​ is beyond the capacities of my laptop.
But !!! The time spent is mainly due to the equality test

## @Carl Love  Maybe. It depends on w...

Maybe.
It depends on what kristavaldes uses as a generator.
There will be no problem if generator=a..b where a abd are two floats
If a and b are integers (let's say a=0 and b "large") problems may arise, but they are more problems from manipulating large numbers than really problems that could be solved by using expand :

restart:
with(LinearAlgebra):
N := 100:
A := RandomMatrix(3, 3, generator=0..N):
CodeTools:-Usage(Eigenvalues(A)):
memory used=6.37MiB, alloc change=32.00MiB, cpu time=134.00ms, real time=2.89s, gc time=7.14ms

N := 10^100:
A := RandomMatrix(3, 3, generator=0..N):
CodeTools:-Usage(Eigenvalues(A)):
memory used=11.06MiB, alloc change=32.00MiB, cpu time=215.00ms, real time=447.00ms, gc time=20.99ms

N := 10^1000:
A := RandomMatrix(3, 3, generator=0..N):
CodeTools:-Usage(Eigenvalues(A)):
memory used=47.43MiB, alloc change=0 bytes, cpu time=6.18s, real time=6.51s, gc time=0ns

On my laptop (Mac OS X, Maple 2015.2, 8 Gb of memory) N = 10^(10^4) generates a "connection to server lost"

## @tomleslie  Thanks for the advices...

@tomleslie

I know very few about the plottools package beyond some functions. I didn't even know that plot displacement functions had been developped in Maple.
I see you use it to modify the tickmarks too, something I never find obvious.

I'm going to see right now how I can still improve your last plot.

Thanks again

## @Preben Alsholm  You're ri...

You're right, a well posed BC problem invoking a third order ODE cannot contain a third order derivative.
So  D[1,1,1](f)(6)=2 is probabibly a tipo error

Such an approximate solution that verifies the boundary conditions is guess := eta -> eta+(1/72)*eta^4
Probably a poor guess for the error message now becomes (I don't understand it)

dsolve( {ode,  D(f)(0)=1, D[1,1](f)(0), D[1, 1, 1](f)(6)=2}, numeric, range=0..6, approxsoln=[f(eta)=guess(eta)]);
Error, (in dsolve/numeric/BVPSolve) unable to store '((D@@3)(HFloat(HFloat(undefined))))(6)-2.' when datatype=float

For the bc (D@@2)(f)(6) = rhs(bc3) an initial guess could be  guess := eta -> eta+(1/18)*eta^3
I tried to modify the problem by setting f(eta) = g(eta)+guess(eta).
In this case g is the solution of an ODE with homogenous bcs  D(g)(0) =0, (D@@2)(g)(0)=0, (D@@2)(g)(6) = 0

I had hoped it could help but I obtain a

Error, (in dsolve/numeric/bvp) initial Newton iteration is not converging

error I don't kow how to circumvent

Last but not least(but you already noticed that), the same ODE is solved in  Numerical Solution of Difficult ODE Boundary Value Problems help page (equation 12) but with other bcs (equation 13)

## @acer    I'm really not ...

I'm really not familiar with subsindets.
I see this command is often used here and that it seems very powerfull. I guess I have to learn more about it and begin to understand how to use it.
Your last file will help me a lot.

Than you Acer

## @acer Thank you Acer for all those ...

@acer

Thank you Acer for all those details.

You write < Your question mentioned the possibility to "modify the ellipses", but without "ad hoc modifications" >. : I wanted to say that it is always possible to "extract" the content of a plot through plottools:-getdata, to remove the undesired part from the returned list (or to modify parts of it), and eventually using PLOT to obtain something that meets out expectations.
But it can be lengthy and sometimes written in a hurry, to answer a specific issue (so the "ad hoc" term).
More of this,    P := DrawGraph(...):   MyPlot:=plottools:-getdata(P):  doesn't contain all the informations about P: for instance there is colouring information.
That is why modifying MyPlot is rather lengthy.

<  For example, do you mean that you would only want the ellipse re-jigged to the longest label if it could be done fully automatically and robustly? > ... the truth is that I'm interested in any idea that could help me to "adjust" the result of DrawGraph, to a reasonable extent.
Among these adjustements are the aspect ratios of the ellipses (maybe rectangles instead), the font (thank you for the tip), the linestyle of the edges/arcs (I'm not yet familiar with Maple 2018  but I remember that, even this last point, was not obvious with Maple 2015).

Thanks again

## @Carl Love I had thought to that, &...

@Carl Love

I had thought to that,  of course.
But inuitively I was certain that it would not be very efficient because of this loop, you know sort of "flee the loop like the pleague".

After a few tests your solution appears to be faster than mine's.
More of this, a slight modification of your's enables computing the first hitting times for increasing calues of c in a very efficient way (obviously c2 > c1 ==> k2 > k1)

FirstHit:= proc(L::{Array,list}(float), C::realcons, n::posint)
local c:= evalf(C), k;
for k from n to numelems(L) do
if L[k] > c then return k fi
od;
FAIL
end proc:

h := FirstHit(XC, 300, 1);
FirstHit(XC, 305, h);

Thank you Carl

## @Joe Riel  Thank you for the confi...

@Joe Riel

Thank you for the confirmation.
Unfortunately I don't have the possibility to test this at home (Mac OSX) but I have some notes in hand from which I can give you some more of the details.
In fact the problem with ExcelTools:-Export occured when we used cmaple with Windows 7 (we used Windows XP before that). With Windows XP Excel files are ".xls" files while they are ".xlsx" since Windows 7.
ExcelTools:-Export(MyMatrix, "MyFile.xls") runs correcty, but ExcelTools:-Export(MyMatrix, "MyFile.xlsx")  doesn't.

Then you're right, ExcelTools:-Export runs correctly when launched from cmaple, the "xlsx file problem" is another issue.

Thanks again

## @Joe Riel  With all due respect, a...

With all due respect, are you sure this is the case for every operational system ?

(If I'm not mistaken I faced problems with Maple 2015 and 2016 when using ExcelTools with cmaple ... but I could be wrong)

## Insufficient...

@vv

SolveTools:-SemiAlgebraic({x^2 + y^2 + a*x*y<z},[x,y], 'parameters'=[a,z]) i

is not sufficient: the good command is

SolveTools:-SemiAlgebraic({x^2 + y^2 + a*x*y<z, x > 0, x < 1, y > 0, y < 1},[x,y], 'parameters'=[a,z])

Now try and run this command without no more information, that is for any real value a....
The memory used to compute the solution (CodeTools:-Usage) is larger than 15 GiB, giving a small idea of how complex the solution is.

Let's try a simple thing with a=4/3, as set in my previous file

SolveTools:-SemiAlgebraic({x^2 + y^2 + 4/3*x*y<z, x < 1, y > 0, y < 1},[x,y], 'parameters'=[a,z]);   ?

and give me please the formal expression of "the (double) integral (not necessarily a closed form, of course, even if such a closed form exists)." given the output of SemiAlgebraic.
If you can and Maple can't we can say Maple has some weaknesses. But if you can't how could Maple ?

I think you deeply underestimate the difficulty of the problem

SemiAlgebraic.mw

BTW: I probably made a mistake in my previous post. I think I have confused z and sqrt(z) in the reduced form of the ellipses. I can correct the file ToyProblem.mw if needed

## @vv Not really a severe limita...

@vv

Not really a severe limitation of PDF.
Or you have to consider in this case that "int has severe limitations" too ...

The problem here reduces to the integration of some function z=f(x,y) over a domain D(x,y | z) where "|" stands for "conditional to the value z of Z".
Generally difficult, without even considering the fact that the CDF of a sum of independent random variable already invokes a convolution product.

Let Z = X^2 + Y^2 + a*X*Y.

First point: it is easy to show that the level curves Z = constant are ellipses as long as  |a| < 2 (hyperbolas for |a| > 2 and a degenerate conocs for |a|=2).

Second point: there are a few obvious cases that makes obvious the computation of  Prob(Z < z)

• a=2 (Z = (X+Y)^2) here the level curves of Z are parallel lines of slope -1
• a=-2 (Z = (X-Y)^2)  here the level curves of Z are parallel lines of slope -1
• a=0 (Z = X^2 + Y^2)  here the level curves of Z are concentric circles of center (0,0)

Apart those 3 elementary situations (Maple 2015 doesn't even solve without a little bit of help), let us consider a  simpler problem than your original one in order to introduce a geometrical approach.
Let X any be two Uniform random variable over [-H, H] where H is some strictly positive number.
Then, provided z is "sufficiently small" regarding to H, Prob(Z < z) is just the surface of some ellipsa.
Under the change of variables {X=U*cos(Pi/4)+V*sin(Pi/4), Y=-U*sin(Pi/4)+V*cos(Pi/4)} Z = X^2 + Y^2 + a*X*Y
writes Z = (1-(1/2)*c)*U^2+(1+(1/2)*c)*V^2.
With  A = 1/sqrt((1-(1/2)*c)) and B = 1/ sqrt((1+(1/2)*c)) this relations becomes Z = U^2/A^2+V^2/B^2
Then Prob(Z < z) = Pi*(A*z)*(B*z).
This result holds only if ellipsa U^2/A^2+V^2/B^2 = z^2 is entirely contained in the unit square [-H,  H]^2.

Now a more complex case.
First observe that U^2/A^2+V^2/B^2 = z^2 describe ellipses centered at point (0,0): then there exist no complete ellipsa within the square (X,Y)=[0, H]^2 and Prob(Z < z) is the measure of some angular sector of the ellipsa defined by angles alpha and beta.
A general formula for this measure, with a proper definition of alpha and beta, is
A*B/2*(arctan(A/B*tan(alpha)) - arctan(A/B*tan(beta))) (... = Prob(Z < z) )
Even if  it is not very difficult to compute "by hand", you already see that some auxiliary difficulties occur:

• express Z = X^2 + Y^2 + a*X*Y in the more suitable form Z = (1-(1/2)*c)*U^2+(1+(1/2)*c)*V^2
• identify the ellipsa U^2/A^2+V^2/B^2 = z^2
• compute alpha and beta

It is not very surprising that PDF (or int) is not capable to do tje job without any human help.

To give you a better idea of the details of the computations you can look to the attached file.
It doesn't contain a complete solution for any value of a in [-2, 2], but just an example for the case a=4/3 that PDF doesn't handle (generalisation of the given example is straightforward).

If you want a more general solution for any value of "a", you have to analyse the different situations when the iso-probability curves are hyperbolas.
At the end a complete solution of your toy-problem becomes a very complex issue.

﻿