Preben Alsholm

13728 Reputation

22 Badges

20 years, 254 days

MaplePrimes Activity


These are replies submitted by Preben Alsholm

@Carl Love Yes, indeed it is faster. By a factor 4 on my machine.

@Carl Love Using the code I gave in my reply "Limiting behavior" I found the following modification of FindSingularity useful.
The purpose is to distinguish between maxfun being exceeded and a singularity.

q:=proc(t) sol(parameters=[t]); sol(-2) end proc;

FindSingularity:= proc(q, t::numeric)
     try q(t)
     catch "cannot evaluate the solution further left of %1, probably a singularity":
          return lastexception[-1];
     catch "cannot evaluate the solution":
     end try;
     FAIL
end proc:

FindSingularity(q,-0.6504); #Testing
# Now trying to narrow down the t-value at which a singularity begins to occur
seq([-0.65-0.00005*i,FindSingularity(q,-0.65-0.00005*i)],i=0..8);
seq([-0.65025-0.00001*i,FindSingularity(q,-0.65025-0.00001*i)],i=0..8);
# Then one could continue with
p2:=proc(t,itvl::range) sol(parameters=[t]); plots:-odeplot(sol,[[z,u(z)],[z,v(z)]],itvl,_rest) end proc;
plots:-animate(p2,[t,-2..0,view=0..2,caption=typeset("t = ",t)],t=-0.65030..-0.65029,paraminfo=false);





@patient Regardless of the relevance of the following, it is interesting (I think) to observe the changes as t passes from -0.6504 to -0.65 and then later the convergence of the curves to the solution to the system evaluated at t=0 (sys0 below).

restart;
mu := 0.1e-2; nu := .1; sigma := 1; k := 1;
ode1:= mu*(diff(u(z),z,z))/t+(z-2*sigma*(diff(1/v(z)^(1/2), z)))*(diff(u(z), z))+(3-2*sigma*(diff(1/v(z)^(1/2),z,z)))*u(z) = 0;
ode2:= nu*diff(v(z),z,z)-z*t*diff(v(z), z)-2*t*v(z)-k*v(z)^(1/2)*u(z)=0;
ics:= u(0)=1, v(0)=1, D(u)(0)=0, D(v)(0)=0:
sol:= dsolve({ode1,ode2,ics}, numeric, parameters=[t]);
p:=proc(t) sol(parameters=[t]); plots:-odeplot(sol,[[z,u(z)],[z,v(z)]],-1..0,_rest) end proc;
p(-1,color=[red,blue]); #Just testing p.
# Various animations in t.
plots:-animate(p,[t,view=0..2],t=-1..-0.3);
plots:-animate(p,[t,view=0..2],t=-0.6504..-0.65);
plots:-animate(p,[t,view=0..2],t=-1..-0.01);
# The system for t=0:
sys0:=diff(u(z),z,z)=0,eval(ode2,t=0);
sol0:=dsolve({sys0,ics},numeric);
plots:-odeplot(sol0,[[z,u(z)],[z,v(z)]],-1..0,thickness=2,linestyle=3,view=0..2); p0:=%:
plots:-animate(p,[t,view=0..2],t=-.001..0,background=p0);


Actually the following is simpler and (I think) better.
The procedure Frem returns unevaluated if it doesn't receive numerical input.
Notice the use of the 'known' option in dsolve.

restart;
Frem:=proc(t,u)
   if not [t,u]::list(numeric) then return 'procname(_passed)' end if;
   frem(t,u)
end proc;

sol:=dsolve({diff(x(t), t)=eval(piecewise(0 <= t and t < 10/3, 60, 10/3 <= t and t < 13/3, 0),t=Frem(t,13/3)),diff(y(t),t)=eval(piecewise(0 <= t and t < 4, 50, 4 <= t and t < 5, 0),t=Frem(t,5)),x(0)=0,y(0)=0},numeric,output=listprocedure,known=[Frem]);
plots:-odeplot(sol,[[t,x(t)],[t,y(t)]],0..50);

@maple fan 
restart;
p:=proc(N,t,Y,YP) local tt;
   YP[1]:=eval(piecewise(0 <= tt and tt < 10/3, 60, 10/3 <= tt and tt < 13/3, 0),tt=frem(t,13/3));
   YP[2]:=eval(piecewise(0 <= tt and tt < 4, 50, 4 <= tt and tt < 5, 0),tt=frem(t,5))
end proc;
sol:=dsolve(numeric,procedure=p,number=2,start=0,procvars=[x(t),y(t)],initial=Array([0,0]));
plots:-odeplot(sol,[[t,x(t)],[t,y(t)]],0..50);

###  It is striking how slow the above is. The following modification speeds it up quite a bit:

restart;
p:=proc(N,t,Y,YP) local fr1,fr2;
   fr1:=frem(t,13/3);
   fr2:=frem(t,5);
   YP[1]:=piecewise(0 <= fr1 and fr1 < 10/3, 60, 10/3 <= fr1 and fr1 < 13/3, 0);
   YP[2]:=piecewise(0 <= fr2 and fr2 < 4, 50, 4 <= fr2 and fr2 < 5, 0)
end proc;
sol:=dsolve(numeric,procedure=p,number=2,start=0,procvars=[x(t),y(t)],initial=Array([0,0]));
plots:-odeplot(sol,[[t,x(t)],[t,y(t)]],0..50);




@patient I'm still confused about z and t. Your original problem, which Carl Love answered so nicely, was to find the value of z for which there was a singularity. That at least was how I understood it. Now you are introducing a parameter t (referred to as time) and are asking for the time (t) of the appearance of a singularity. But for each given t (or each in some interval) isn't there a singularity happening at some z (dependent on t)?

Rereading your original question you talk of blow-up. But the singularity encountered by dsolve/numeric is not a blow-up.
It is a value of z for which v(z) = 0, so integration has to stop.

In the present context (your uploaded worksheet) with t involved and set to -1 the plot

plots:-odeplot(res,[[z,u(z)],[z,v(z)]],-1..1);

shows what happens:

Integration stops at z = 0.42086521 going right and at z= -0.42086521 going left because v(z) (the blue graph) is approaching zero., which clearly is a problem for your system.

To take a very simple example exhibiting the same kind of singularity:
restart;
ode:=diff(y(x),x)=-y(x)^(-1/2);
res:=dsolve({ode,y(0)=1}); #Symbolic solution: y(x) = (1/4)*(-12*x+8)^(2/3)
sol:=dsolve({ode,y(0)=1},numeric);
plots:-odeplot(sol,[x,y(x)],0..1); #Singularity at x = 2/3



@patient If this is to be considered an ODE system with a parameter t that would be OK, but initial conditions would have to be given at a particular value of z.

Actually the solution seems to be found, but rejected in favor of a more general "solution":

restart;
debug(solve);
solve({abs(a-b)=0, sqrt(2*b+c)=0, c^2-c+1/4=0});

This seems to be used by solve (and works):

SolveTools:-Identity({c^2-c+1/4 = 0, abs(a-b) = 0, sqrt(2*b+c) = 0},{},{a,b,c});



@patient You have a variable t in ode1, should that be z? You also refer to t (and T) later.
Could you clarify?

There are some missing parentheses. I suppose what you intended was:

ode1:= -0.1*diff(u(z),z,z)+(z-2*diff(v(z)^(-1/2),z))*diff(u(z),z)+(3-2*diff(v(z)^(-1/2),z,z))*u(z)=0;
ode2:= 0.1*diff(v(z),z,z)+0.01*z*diff(v(z),z)+0.02*v(z)-u(z)*v(z)^(1/2)=0;
ics:= u(0)=0, v(0)=0.1, D(u)(0)=0, D(v)(0)=0;

With the zero initial conditions you clearly have a problem just getting started.

@khoirun When I run your code in Maple 2015 I don't get any complaints about singularity.
Be aware that fieldplot3d and display belong to the plots package.

However, if you extend the t-interval to the left e.g. to -10 then you get warnings about singularities.

You could try the first initial condition in IC and use dsolve/numeric like the following.
The construction {SYS[],IC[1][]} just gives you a set containing the odes and the initial condition.

res:=dsolve({SYS[],IC[1][]},numeric);
plots:-odeplot(res,[t,y(t)],-6..5);
res(-5.1383855);

#There is most likely indeed a singularity at about t=-5.1383855. What this means here is that the orbit "reaches" infinity in a finite time. There is nothing weird about that.
The simple example dsolve({diff(y(t),t) = -y(t)^2, y(0)=1}) has a singularity at t=-1.


@yellowcanary If you want to switch to using 1D input either permanently (i.e. until you actively change back) or for just one session, you can do that under the menu item Tools/Options/Display/Input Display.
You should know that Maple Notation is the same as 1D.
Then go to Interface in the same dialog box and under Default format for new worksheets choose Worksheet.
After that click the button Apply Globally (or Apply to Session).

@In-Jee Jeong Since it is not clear what you are actually doing I think you ought to present us with the details, perhaps as an uploaded worksheet.

It is not clear what you are actually doing because in the question itself the word 'numeric' does not appear in the dsolve command as tomleslie points out, and yet you are using odeplot there.
Secondly, in your reply to bogo you write plot(b[0](t),t=0..1). For a numerical solution that will only work if you are using output=listprocedure and b[0] is assigned one of the solution procedures.

@tomleslie I agree. The problem is (again) the 2D input, where spacing is interpreted as multiplication, which apparently is indeed intended in this case. But then there are (it appears) two cases of missing a space after PB.

Since everybody knows that if an expession F is an antiderivative of the expression f on some interval I then so is F+C, where C is any constant, there is no need for Maple to add a constant. What would be the purpose? Knowing one antiderivative (on some interval) means knowing all.
You just have to know that Maple doesn't add such a constant.

I noticed that in the help page for int in Maple 2015 this is pointed out explicitly.

First 121 122 123 124 125 126 127 Last Page 123 of 230