Preben Alsholm

MaplePrimes Activity


These are replies submitted by Preben Alsholm

@roya2001 The new error is clearly due to the fact that no successes were encountered, thus res is still just a name, not tyhe name of a table.
What did you change since the last version? I cannot keep track of it.

@fbackelj If you are the teacher then it seems to me to be your job to explain the situation to your students after having read the help for LinearAlgebra:-ConditionNumber. It uses the word "estimates" and that surely indicates that you are getting an estimate of the condition number. What are your students using the condition number for? Do they need the exact value for anything?

@Thomas Richard Clearly there could be a bug somewhere, but if I understand the statement in the link given by acer correctly, then there need not be.

The statement in https://www.nag.co.uk/numeric/fl/nagdoc_fl24/html/F07/f07agf.html

says:

"The computed estimate RCOND is never less than the true value ρ, and in practice is nearly always less than 10ρ, although examples can be constructed where RCOND is much larger."

Since the value returned (RCOND) is an estimate of 1/c where c is the (true) condition number, I understand this statement as saying that est_c <= c (always), est_c > c/10 (almost always), but that there are examples where est_c << c/10.

(By est_c I mean the estimated condition number).

####### You may try the following and notice (as Markiyan did) that most often there is agreement, but fairly often est_c is less than c (and not because of rounding).

restart;
eps:=0.1;
N:=0;
R:=Vector():
to 1000 do
use LinearAlgebra in
  A:=RandomMatrix(4,generator=-9.0..9.0);
  Ec:=ConditionNumber(A);
  c:=Norm(A)*Norm(A^(-1))
end use:
  if (c-Ec)/c>eps then
    N:=N+1;
    R(N):=Ec/c;
    print(Ec,c,'ratio'=R(N))
  elif (c-Ec)/c<-10^(-Digits+3) then
    print("Shouldn't happen",Ec,c)
  end if;
end do:


In a trial run I found that in 12 percent of the cases (c-Ec)/c > 0.1 and no "Shouldn't happen".
The ratio min and max
(min,max)(R);
were roughly 0.4 and 0.9.
This is in full compliance with the cited accuracy statement.

After answering and commenting your question at

http://mapleprimes.com/questions/207601-Remove-Assumptions-

I got think that one reason for the statement "This functionality is not intended for end users" may be that complication with assumed variables when saving as an .m file instead as an .mpl file.

Below I was trying to answer your question as formulated.
If you use an .mpl extension instead all your problems are gone. There is nothing to remove.

@Ramakrishnan Certainly this a bug. Kitonum didn't enter contradicting assumptions. The inequalities

x>-4*Pi, x<-5*Pi/2

are clearly NOT contradictory, and as Kitonum is pointing out sin(2*x)/cos(x+3*Pi/2)=1 has the unique root x=-11*Pi/3 satisfying those inequalities.

That solve can handle Pi is seen here:

solve([cos(x)=1/2,x>0,x<7*Pi], x, allsolutions,explicit); #Take a look at the correct result

You say

" Experts have made rule for maple to follow and maple follows meticulously!."

I won't disagree with that, but just point out that even experts (being human) make mistakes.

@Markiyan Hirnyk Thanks for your confidence in my abilities. But the problem you refer to was a simple logical error anybody would do on occasion. Too bad it wasn't caught in testing.

Getting solve to return correct solutions for equations with assumptions is clearly a hard task. So no wonder that incomplete solutions are returned as in my example with cos(x)=1/2. But see below for more on that.

The same solutions are returned in Kitonum's example if sin(2*x)/cos(x+3*Pi/2)=1 is replaced by simplify(sin(2*x)/cos(x+3*Pi/2))=1. No error.

I tried the following:
#####
restart;
res:=solve([simplify(sin(2*x)/sin(x))=1,x>-4*Pi,x<-5*Pi/2], x, allsolutions);
about(_Z1);
about(_B1);
##So since  _B1 and _Z1 are integers we have
_B1-3*_Z1<=6 , _B1-3*_Z1>4;
isolve({_B1-3*_Z1<=6 , _B1-3*_Z1>4,_B1>=0,_B1<=1}); #A warning issued
eval(res,%); ##The correct unique solution
#####################################################
##Now without simplify; cos(x+3*Pi/2) is simplified just by evaluation to sin(x), so I use that.
restart;
res:=solve([sin(2*x)/sin(x)=1,x>-4*Pi,x<-5*Pi/2], x, allsolutions); ##The error reported by Kitonum
about(_Z1);
        Originally _Z1, renamed _Z1~:
  is assumed to be: -2
### Remarkable: that was what we found that the _Z1 in the preveious session should be

about(_B1);
_B1:
  nothing known about this object


@law47 

L:='L': #If L is previously assigned

Q(L)(x,y,y1);
      (D[1](L))(x, y, y1)+y1*(D[2](L))(x, y, y1)+y*y1*(D[3](L))(x, y, y1)/x

No problem.

As the one who supplied the code for Q, I would like you to clarify what you want. After all even if L is not defined to be anything, this works:

Q(L)(x,y,y1);
    (D[1](L))(x, y, y1)+y1*(D[2](L))(x, y, y1)+y*y1*(D[3](L))(x, y, y1)/x

In the question in

http://mapleprimes.com/questions/207524-Writing-An-Operator

you asked

"I am trying to write 

D= d/dx + y1*d/dy + y*y1/x*d/dy1 

So that I can apply D to some function L(x,y)."

I think I answered that.
But did you really mean something else?
Or did you want to make a more general (higher level) procedure taking in a particular expression a(x,y,y1)*d/dx +b(x,y,y1)*d/dy+c(x,y,y1)*d/dy1 and as output having a procedure like Q  working on functions?
In our particular case a(x,y,y1)=1, b(x,y,y1)=y1, c(x,y,y1)=y*y1/x.
If so that should be pretty easy.

Yes, this is a mess:

solve([cos(x)=1/2, x>0, x<2*Pi], x, allsolutions, explicit); #OK
      {x = (1/3)*Pi}, {x = (5/3)*Pi}
solve([cos(x)=1/2,  x>-4*Pi, x<-5*Pi/2], x, allsolutions, explicit); #Ignoring the inequalities
      {x = -(2/3)*Pi*_B2+(1/3)*Pi+2*Pi*_Z2}

@Dmitry I'm not about to speculate on this. I'd like to see an example of what you mean by really different?

@Dmitry The only difference between the 2 regimes is the appearance/disappearance of upsilon in the K-equation. That is handled by piecewise. Thus there is only one system.

Notice now that dsolve/bvp works with the change from psi[S] to psi[s] as mentioned in my answer above.

You want to apply this D to a function L of 2 variables (x and y). So what is d/dy1 supposed to mean?

@Carl Love Thanks for calling my attention to frem.

I suppose something like
mp:=(x,p)->frem(x+p/2,p)+p/2;

might work OK in the present context (if 2*Pi is replaced by evalf(2*Pi)).

That mp(8,2) returns 2. and not 0. seems to be irrelevant here.

First 96 97 98 99 100 101 102 Last Page 98 of 231