Preben Alsholm

13728 Reputation

22 Badges

20 years, 250 days

MaplePrimes Activity


These are replies submitted by Preben Alsholm

Yes, your experience with increasing Digits is one I have had often in bvp-problems. It is a wild goose chase. Generally if I have Digits at 15, I give up raising it because the goose is always going to be ahead of me.
So the solution is generally to make use of other options.

In the present case, I quite agree with I_Mariusz that you should not have right endpoint at 20, but considerably lower. He chose 5, I had success also with 7.
The problems appear because U and V become roughly constant so that U'(eta) = V'(eta) = 0 (almost). When solving for the highest derivatives the quantity U'(eta)^2 +V'(eta)^2 appears in the denominators. That will cause numerical problems although the numerators also will be approaching zero (and I think at a higher rate).

@vv I must have been working on the expanded version of my answer while you posted your reply.
Yes, even for that concrete value of y we don't get a definite answer.

I copied your code into Maple 2015.2 and it ran without any error.
After that I did:
seq(P[j][1,1],j=1..n-1);

and I received as expected 1,1,1.

Did you run the code you showed us just after a restart?
If not try that.

That should be possible in Statistics:-Fit.

?Statistics:-Fit

@roya2001 I don't see any way to get an answer that looks like your approximate solution. That obviously doesn't mean that such a solution doesn't exist, however.

######## An observation:
Let S and G be the expressions for your approxsoln for s(x) and g(x):

S:=cosh(upsilon*x)-cos(upsilon*x)-(cosh(upsilon)+cos(upsilon))*(sinh(upsilon*x)-sin(upsilon*x))/(sinh(upsilon)+sin(upsilon));
G:=sin(((2*n+1)*(1/2))*Pi*x);
##Now one of your boundary conditions is D(g)(1)+1/2*(D(s)(1))^2=0.
#We find
eval(diff(G,x),x=1); # Zero
eval(diff(S,x)^2/2,x=1); #46.0533267839276

Thus that boundary condition is very far from being satisfied by your proposed approximate solution.


Besides what J4James mentions, you also need to enclose the sequence of pdes in {}:
{PDE1, PDE2, PDE3}

@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.

First 95 96 97 98 99 100 101 Last Page 97 of 230