Preben Alsholm

13743 Reputation

22 Badges

21 years, 117 days

MaplePrimes Activity


These are replies submitted by Preben Alsholm

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.


Did you mean

ode1:=(1/r^2)*diff(r^2*diff(u(r,t),r),r) =v(r,t);

ode2:= diff(v(r,t),t)=u(r,t);

IBC:={u(0.2,t)=1,u(0,t)=0.5,v(r,0)=0.01};

?



 

@shinelookat You can use odeplot like this,

plots:-odeplot(sol1,[x(t),y(t),z(t)],0..10,axes=boxed,thickness=3);

Animating the orbit can be done by asking for a number of frames as in

plots:-odeplot(sol1,[x(t),y(t),z(t)],0..10,axes=boxed,thickness=3,frames=50);

@shinelookat You can use odeplot like this,

plots:-odeplot(sol1,[x(t),y(t),z(t)],0..10,axes=boxed,thickness=3);

Animating the orbit can be done by asking for a number of frames as in

plots:-odeplot(sol1,[x(t),y(t),z(t)],0..10,axes=boxed,thickness=3,frames=50);

@Vaclav 

If you try

g:=sum((-p)^n*n/(factorial(k-n)*factorial(n)), n = 1 .. k);

you will see that p=1 is an exceptional case. It could be handled by taking a limit:

limit(g,p=1) assuming k::posint;

but notice that you get the generic result (not valid for k = 1) as mentioned earlier.

@Vaclav 

If you try

g:=sum((-p)^n*n/(factorial(k-n)*factorial(n)), n = 1 .. k);

you will see that p=1 is an exceptional case. It could be handled by taking a limit:

limit(g,p=1) assuming k::posint;

but notice that you get the generic result (not valid for k = 1) as mentioned earlier.

First 207 208 209 210 211 212 213 Last Page 209 of 232