Preben Alsholm

MaplePrimes Activity


These are replies submitted by Preben Alsholm

Is this actually in Maple?
I tried the following and it works as expected:
 

A:=LinearAlgebra:-RandomMatrix(3,generator=-99..99);
B:=LinearAlgebra:-RandomMatrix(3,generator=-99..99);
A.B;
A[2,3]:=777;
A;
A.B;

 

I tried Generate PDF under More ... . I asked that the file be opened in Foxit PhantomPDF. I received the error:
Format Error: Not a PDF or corrupted.

It appears to have 0 bytes.

@Rouben Rostamian  It may be worth while posting a proof that the present bvp problem doesn't have a solution y(x) that is positive for all 0 < x < 2945.
##
Notes:
The inequality ineq4 below will give L__crit when solved for L in terms of b=D(y)(L).
The inequality becomes an equality when D(y)(0) = y1 = 0.
##

## Proof that the given problem ode with the given bcs has no solution y(x) > 0 for all x > 0, x<=2945.

restart;
ode:=diff(y(x),x,x)=a*y(x)^q;
ics:={y(0)=0, D(y)(0)=y1}; # y1 >= 0 is necessary to ensure y(x) > 0 for x > 0.
## We assume q > 0, q < 1, a > 0, and that y(x) > 0 for x > 0.
## Thus y''(x) > 0, and y1 >= 0 is necessary to ensure y(x) > 0 for x > 0.
##
diff(y(x),x)*ode;
eq:=map(int,%,x)+(0=C/2);
eval[recurse](convert(eq,D),ics union {x=0});
eq1:=eval(eq,C=y1^2);
solve(eq1,{diff(y(x),x)})
## We need y'(x) >=0 all x, so pick the first:
ode1:=simplify(op(%[1])) assuming q>0;
## Now estimate
rhs(ode1) > eval(rhs(ode1), y1=0);
## So we have
lhs(ode1)> eval(rhs(ode1), y1=0);
ineq1:=simplify(%) assuming a>0;
## Thus by separation:
diff(y(x),x)*y(x)^(-(q+1)/2)>sqrt(2*a/(q+1));
map(int,%,x=0..X,continuous);
ineq2:=eval(%,{X=x,y(0)=0});
ineq3:=y(x)>(lhs(ineq2)*(1-q)/2)^(2/(1-q));
## So y(L) (L>0) is at least
ymin:=eval(lhs(ineq3),x=L);
eval(convert(ineq1,D),x=L);
## So the following inequality holds:
ineq4:=simplify(eval(%,y(L)=ymin)) assuming q>0,q<1,a>0,L>0;
## Now in the present case we have:
params:={a=0.00003019,q=0.337,x=2945,L=2945};
evalf(eval(ineq4,params)); #Result 0.2984244716 < D(y)(2945)
## But we also have
D(y)(2945)=0.0116;
## A contradiction. So we cannot have y(x) > 0 for 0 < x < 2945.

Now the inequality ineq4 in our case was 0.2984244716 < D(y)(2945), so it could be fun to see if dsolve/numeric/bvp can handle a value for D(y)(2945)  like 0.3 which is not much higher than 0.2984244716:
 

## Trying D(y)(2945)=0.3 (slightly higher than the lower limit found above.
## Using continuation:
ode1:=diff(y(x),x,x)=0.00003019*abs(y(x))^(c*0.337);
res:=dsolve({ode1,y(0)=0,D(y)(2945)=0.3}, numeric,method=bvp[midrich],initmesh=8192,continuation=c,abserr=0.3e-5,mincont=1/300);
plots:-odeplot(res,[x,y(x)],caption=typeset(D(y)(2945)=0.3)); #Looks good!

@Rouben Rostamian  Very nice, and the solution is nonnegative for all x, whereas my numerical solution is negative for all x > 0.
I was just about to give a proof that no solution exists with y(x) > 0 for 0 < x <2945 (in the present case).
Making use of the nonuniqueness of solutions starting with y(x0) = y'(x0) = 0 is a great idea.

 

@Mehdi Khalifeh Reading the description you provide I see that H(t) = (H1(t), 0, 0) and that M'(t) = (0, M2'(t), M3'(t)).
Thus H(t) . M'(0) = 0 for all t and so E = 0.
???

@dineshchawde You gave us an image. Thus to experiment with it we would have to type the whole thing ourselves. Needless to say, we don't like that. So please use text or upload a worksheet using the fat green arrow in the MaplePrimes editor.

@Axel Vogt Yes, following your rewriting we find that the integral over the fourth quadrant is given by 26.3.19 in A&S:
 

AS:=1/4-arcsin(B/sqrt(G))/2/Pi; #rho = -R/sqrt(G) in A&S 26.3.19
## Value from the guess:
v4:=arctan(sqrt(-B^2+G),B)/2/Pi;
##
df:=simplify(v4-AS) assuming G>B^2;
convert(df,arctan);
simplify(%) assuming G>B^2,B>0; # 0
simplify(%%) assuming G>B^2,B<0; # 0

 

@NorwegianStudent If a given function is already periodic you don't have to do anything: just plot it.
If, on the other hand, you intend to extend a function given on some interval a..b to a function periodic with period b-a then you could just use the following procedure PeriodicExtension:
 

PeriodicExtension:=proc (fip, itvl::{range, name = range}) local x, f, a, b, p, pr; 
   description "Input a procedure or an algebraic expression and an interval on the form a..b or x = a..b on which the function is given by the procedure/expression.\nOutput is a procedure/expression that is the periodic extension with period b-a."; 
  if type(itvl, range) then 
    f := eval(fip); 
    a, b := op(itvl) 
  else 
    x := lhs(itvl); 
    f := unapply(fip, x); 
    a, b := op(rhs(itvl)) 
  end if; 
  p := b-a; 
  pr := proc (x) f(x-floor((x-a)/p)*p) end proc; 
  if type(itvl, range) then eval(pr) else pr(x) end if 
end proc;
###
### Two examples:
f:=PeriodicExtension(x^2,x=0..1);
plot(f,x=-3..3,discont=true);
##
f:=PeriodicExtension(arcsin,-1..1);
plot(f,-3..3,discont=true);

To get the description do
Describe(PeriodicExtension);

Could you give a simple example of what you mean by "segregate different terms"?

Two comments:
(1) The following shows that the number of builtin procedures as found by anames depends on what happens in the session.
Thus it is not a reliable measure of the total number of builtin procedures, but it gives a lower bound.
 

restart;
AB:={anames(builtin)};
nops(AB); # 234 in Maple 2017.2
member(zip,AB); #false
op(3,eval(zip)); # builtin = 589 in Maple 2017.2
## Now continue with
U:=Vector([1,2,3]); #Causes one more builtin to pop up:
AB2:={anames(builtin)};
nops(AB2); # 235 in Maple 2017.2
AB2 minus AB;

AB2 minus AB returns (the set of) the procedure `print/rtable/prepro`.
Added: The list B generated by Carl Love's code contains that procedure:
select(has,B,586); #In Maple 2017.2

(2) In zip the truly builtin procedure rtable_zip is used when rtables are arguments for zip:

member(rtable_zip,AB); #true
op(3,eval(rtable_zip)); # builtin = 401 in Maple 2017.2
showstat(rtable_zip); # The expected error

I tried the code above in Maple releases 12, 15, 16, 17, 18, 2015, 2016, and 2017.
For the sake of clarity the code was:

AB:={anames(builtin)};
nops(AB);
member(zip,AB);
op(3,eval(zip));
showstat(zip);

The results were that
1. nops(AB) returned 195, 216, 222, 226, 228, 233, 232, 234.
2. member(zip,AB) returned false in all those releases.
3. op(3,eval(zip)) returned Copyright 1999 in Maple 8, Copyright 2007 in Maple 12, 15, 16, 17, and 18, in Maple 2015: builtin=578, Copyright 2007, in Maple 2016: builtin=585, Copyright 2007, in Maple 2017: builtin=589, Copyright 2007.
4. All but Maple 8 showed 103 lines of code in showstat(zip). In Maple 8 showstat(zip) showed that the real code was put in ZipImplementation:-zip, which has 99 lines of code (as shown by showstat).

In the code for zip in all but Maple 8 (of the ones I looked at) there is a procedure keep_sparse. I wonder what that is doing.

@AndersDD When clicking I get a "Page Not Found" error.
Try again and this time maybe change the name of the file to e.g.  Afleveringssaet-1.maple or whatever doesn't have æ, ø, or å in it.

@tomleslie The integral over the whole plane Maple finds easily, and it is found to be 1.
 

P := (x,y)->(1/2)*exp(-(1/2)*(x^2+G*y^2-2*B*x*y)/(-B^2+G))/(Pi*sqrt(-B^2+G));
A:=Int(Int(P(x,y),x=-infinity..infinity),y=-infinity..infinity);
value(A) assuming G>B^2;

Furthermore, if we know the value of the integral over e.g. the integral over the fourth quadrant (value = v4, say) then the value over the other quadrants (v1,v2,v3) are simply v2 = v4, v3 = eval(v4,B=-B), and v1 = v3.
It is easily verified that if (as I guessed and _Maxim_ confirms above)
v4 = arctan(sqrt(-B^2+G), B)/Pi/2 then
v1+v2+v3+v4 = 1 all assuming that G > B^2.

@tomleslie Yes that specification of the ranges must be interpreted.
But if the whole plane was meant, why write it that complicated when ranges x= -infinity..infinity, y=-infinity..infinity would do the job?

I have had no such problem with Maple 15 nor have I had with Maple 2015 (just in case that is really the one you are talking about?).

First 58 59 60 61 62 63 64 Last Page 60 of 231