Preben Alsholm

MaplePrimes Activity


These are replies submitted by Preben Alsholm

Try raising Digits above 10, then fsolve finds nothing.

restart;
Digits:=11:
fsolve([2.32=a+b,pa=.4310344828*a,pb=.4310344828*b,pa+pb=1]);

(See perhaps the addendum to my answer 'No solution indeed').


Try raising Digits above 10, then fsolve finds nothing.

restart;
Digits:=11:
fsolve([2.32=a+b,pa=.4310344828*a,pb=.4310344828*b,pa+pb=1]);

(See perhaps the addendum to my answer 'No solution indeed').


@icegood Making rather obvious changes in the code for IntegrationTools:-Expand we get the following procedure ExpandSum. To see the code for IntegrationTools:-Expand, do showstat(IntegrationTools:-Expand);

ExpandSum := proc(v)
local f, i, c, g, w;
      if type(v,specfunc(anything,{'Sum','sum'})) then
        f := expand(op(1,v));
        i := op([2,1],v);
        if type(f,`*`) then
          g, c := selectremove(type,f,(':-dependent')(i));
          c*op(0,v)(g,op(2 .. -1,v))
         elif type(f,`+`) then
          return map(procname,map(op(0,v),f,op(2 .. -1,v)))
         else
          return v
         end if
       else
        if not hastype(v,specfunc(anything,{'Sum','sum'})) then
          return v
         else
          return evalindets(v,specfunc(anything,{'Sum','sum'}),'procname')
        end if
       end if
end proc;

@icegood Making rather obvious changes in the code for IntegrationTools:-Expand we get the following procedure ExpandSum. To see the code for IntegrationTools:-Expand, do showstat(IntegrationTools:-Expand);

ExpandSum := proc(v)
local f, i, c, g, w;
      if type(v,specfunc(anything,{'Sum','sum'})) then
        f := expand(op(1,v));
        i := op([2,1],v);
        if type(f,`*`) then
          g, c := selectremove(type,f,(':-dependent')(i));
          c*op(0,v)(g,op(2 .. -1,v))
         elif type(f,`+`) then
          return map(procname,map(op(0,v),f,op(2 .. -1,v)))
         else
          return v
         end if
       else
        if not hastype(v,specfunc(anything,{'Sum','sum'})) then
          return v
         else
          return evalindets(v,specfunc(anything,{'Sum','sum'}),'procname')
        end if
       end if
end proc;

If you want to use lambda as a parameter you somehow have to avoid it appearing as an argument to f as it does in one of the boundary points.

So replace x with x(t) satisfying x'(t) = 1 and x(0) = -lambda, and let g(t) = f(x(t)). Then you get the following system

sys:={g(t)^2*diff(g(t),t,t)+x(t)*g(t)*(1-2*beta)+beta*diff(g(t),t)*(x(t)^2-1)=0,diff(x(t),t)=1};
#Now we use 3 parameters p (as in my previous answer), lambda, and beta.
#Otherwise the idea is as before although there are some changes.
init:=dsolve(sys union {g(0)=beta*(lambda^2-1)/p,D(g)(0)=p,x(0)=-lambda},numeric,parameters=[p,lambda,beta],output=listprocedure);
G,X:=op(subs(init,[g(t),x(t)]));
#We don't really need X, but use it as a check.
q:=proc(P,L,B)
   if not type([P,L,B],list(numeric)) then return 'procname(_passed)' end if;
   init(parameters=[P,L,B]);  
   G(L+1);
end proc;
q(.0236,.3,.01);
#I use the plotting data to find the value of p where G(lambda+1)=0 (i.e. f(1)=0).
plot(q(p,.3,.01),p=.023..0.024,-1..1);
op(indets(%,listlist)):
L:=remove(has,%,HFloat(undefined)):
Ls:=sort(L,(x,y)->evalb(abs(x[2])Ls[1];
p1,q1:=op(Ls[1]);
# fsolve has problems with fsolve(q(p,.3,.01),p=p1), which is not so surprising since f'(1) is infinite when f(1) = 0. So the following is just a check.
fsolve(q(p,.3,.01)-q1,p=p1);
X(0);
G(0);
G(.3+1);
#Another value for lambda, same beta.
plot(q(p,.35,.01),p=0..0.2,-1..1);
op(indets(%,listlist)):
L:=remove(has,%,HFloat(undefined)):
Ls:=sort(L,(x,y)->evalb(abs(x[2])Ls[1];
p1,q1:=op(Ls[1]);
fsolve(q(p,.35,.01)-q1,p=p1);
X(0);
G(0);
G(.35+1);

This method could be automated, but maybe somebody has a better idea.

Added: Instead of the plotting idea which makes use of adaptive plotting you could use the following procedure R, which takes lambda, beta, and a p-interval as its first 3 arguments and outputs a p-value with abs(f(1))< a tolerance, by default 10^(-Digits+2).

R:=proc(L::numeric,B::numeric,interval::range(numeric),{initstep::positive:=1e-2,tolerance::positive:=10^(-Digits+2),iterationlimit::posint:=10})
  local i,j,pp,M,a,b,d,A;
  a,b:=op(interval);
  d:=initstep;
  A:=tolerance+1;
  for j from 1 to iterationlimit while A>tolerance do
    M:=Matrix(0);
    i:=0;
    for pp from a to b by d do
      try
        i:=i+1;
        M(1,i):=pp;
        M(2,i):=q(pp,L,B);
      catch:
        break
      end try
    end do;
    if i=1 then error "Try another interval" end if;
    A:=abs(M(-3));
    a:=M(-4);
    d:=d/10;
  end do;
  if A>tolerance then WARNING("Result not within tolerance") end if;
  M(-4),M(-3);
end proc;
R(.35,.01,0.01..1);
R(.36,.01,0.01..1);
# A loop through lambda-values, initially crude since f(1) may never become zero.
for L from .3 to .5 by .02 do L,R(L,.01,.01..1,iterationlimit=2) end do;
for L from .3 to .37 by .01 do R(L,.01,.01..1) end do;


R(.35,.2,.01..1,iterationlimit=7,tolerance=1e-3,initstep=.01);
R(.35,.5,.01..1,iterationlimit=7,tolerance=1e-3,initstep=.07);

If you want to use lambda as a parameter you somehow have to avoid it appearing as an argument to f as it does in one of the boundary points.

So replace x with x(t) satisfying x'(t) = 1 and x(0) = -lambda, and let g(t) = f(x(t)). Then you get the following system

sys:={g(t)^2*diff(g(t),t,t)+x(t)*g(t)*(1-2*beta)+beta*diff(g(t),t)*(x(t)^2-1)=0,diff(x(t),t)=1};
#Now we use 3 parameters p (as in my previous answer), lambda, and beta.
#Otherwise the idea is as before although there are some changes.
init:=dsolve(sys union {g(0)=beta*(lambda^2-1)/p,D(g)(0)=p,x(0)=-lambda},numeric,parameters=[p,lambda,beta],output=listprocedure);
G,X:=op(subs(init,[g(t),x(t)]));
#We don't really need X, but use it as a check.
q:=proc(P,L,B)
   if not type([P,L,B],list(numeric)) then return 'procname(_passed)' end if;
   init(parameters=[P,L,B]);  
   G(L+1);
end proc;
q(.0236,.3,.01);
#I use the plotting data to find the value of p where G(lambda+1)=0 (i.e. f(1)=0).
plot(q(p,.3,.01),p=.023..0.024,-1..1);
op(indets(%,listlist)):
L:=remove(has,%,HFloat(undefined)):
Ls:=sort(L,(x,y)->evalb(abs(x[2])Ls[1];
p1,q1:=op(Ls[1]);
# fsolve has problems with fsolve(q(p,.3,.01),p=p1), which is not so surprising since f'(1) is infinite when f(1) = 0. So the following is just a check.
fsolve(q(p,.3,.01)-q1,p=p1);
X(0);
G(0);
G(.3+1);
#Another value for lambda, same beta.
plot(q(p,.35,.01),p=0..0.2,-1..1);
op(indets(%,listlist)):
L:=remove(has,%,HFloat(undefined)):
Ls:=sort(L,(x,y)->evalb(abs(x[2])Ls[1];
p1,q1:=op(Ls[1]);
fsolve(q(p,.35,.01)-q1,p=p1);
X(0);
G(0);
G(.35+1);

This method could be automated, but maybe somebody has a better idea.

Added: Instead of the plotting idea which makes use of adaptive plotting you could use the following procedure R, which takes lambda, beta, and a p-interval as its first 3 arguments and outputs a p-value with abs(f(1))< a tolerance, by default 10^(-Digits+2).

R:=proc(L::numeric,B::numeric,interval::range(numeric),{initstep::positive:=1e-2,tolerance::positive:=10^(-Digits+2),iterationlimit::posint:=10})
  local i,j,pp,M,a,b,d,A;
  a,b:=op(interval);
  d:=initstep;
  A:=tolerance+1;
  for j from 1 to iterationlimit while A>tolerance do
    M:=Matrix(0);
    i:=0;
    for pp from a to b by d do
      try
        i:=i+1;
        M(1,i):=pp;
        M(2,i):=q(pp,L,B);
      catch:
        break
      end try
    end do;
    if i=1 then error "Try another interval" end if;
    A:=abs(M(-3));
    a:=M(-4);
    d:=d/10;
  end do;
  if A>tolerance then WARNING("Result not within tolerance") end if;
  M(-4),M(-3);
end proc;
R(.35,.01,0.01..1);
R(.36,.01,0.01..1);
# A loop through lambda-values, initially crude since f(1) may never become zero.
for L from .3 to .5 by .02 do L,R(L,.01,.01..1,iterationlimit=2) end do;
for L from .3 to .37 by .01 do R(L,.01,.01..1) end do;


R(.35,.2,.01..1,iterationlimit=7,tolerance=1e-3,initstep=.01);
R(.35,.5,.01..1,iterationlimit=7,tolerance=1e-3,initstep=.07);

No problem with your last attempt except of course if i doesn't have a numerical value, or if k does have a numerical value.

Example.

restart;
i:=5:
for j to i do c[j]:=plot(x^j,x=0..1) end do:
plots:-display(c[k]$k=1..i);

#Since j has a numerical value (in this case 6) you need unevaluation quotes as shown below, therefore the use of seq in place of $ (as suggested by Clare So) is in general better.
plots:-display('c[j]'$'j'=1..i);

The elementwise operation ~ used in Robert Israel's answer was introduced in Maple 13. You can use map instead.

plots[display](map(plottools[scale],{p1,p2,p3},1/(2*Pi*sqrt(2)),1),
            labels=[x/(2*Pi*'sqrt(2)'),h]);

The elementwise operation ~ used in Robert Israel's answer was introduced in Maple 13. You can use map instead.

plots[display](map(plottools[scale],{p1,p2,p3},1/(2*Pi*sqrt(2)),1),
            labels=[x/(2*Pi*'sqrt(2)'),h]);

Why do you want a perfect echo of the input?

In some cases I see the point. As an example it may be a good idea to assign a differential equation to a variable and ending the assignment with a semicolon so the output is seen. This can prevent attempts to solve the wrong equation.

The purpose of seeing an echo of (5*x+4)/(3*x)=7 could be to catch an error like writing (5*x+4)/3*x=7. If that is the problem you could use 2D input, which I personally never do, however.

Why do you want a perfect echo of the input?

In some cases I see the point. As an example it may be a good idea to assign a differential equation to a variable and ending the assignment with a semicolon so the output is seen. This can prevent attempts to solve the wrong equation.

The purpose of seeing an echo of (5*x+4)/(3*x)=7 could be to catch an error like writing (5*x+4)/3*x=7. If that is the problem you could use 2D input, which I personally never do, however.

@Markiyan Hirnyk I said that the discriminant is 8*l^2 and that the roots are correct, and indeed they are, since (when l is real) {l,-l} and {abs(l),-abs(l)} are the same sets.

Maple's answer to the roots has the advantage that it is correct even if l is not real.

 

@Markiyan Hirnyk I said that the discriminant is 8*l^2 and that the roots are correct, and indeed they are, since (when l is real) {l,-l} and {abs(l),-abs(l)} are the same sets.

Maple's answer to the roots has the advantage that it is correct even if l is not real.

 

The discriminant of the polynomial:

pol:=(2+l^2-4*l)*n^2-(4-4*l)*n+2;
discrim(pol,n);
                                 2
                              8 l
and the roots of pol are correct.


The discriminant of the polynomial:

pol:=(2+l^2-4*l)*n^2-(4-4*l)*n+2;
discrim(pol,n);
                                 2
                              8 l
and the roots of pol are correct.


First 206 207 208 209 210 211 212 Last Page 208 of 231