Preben Alsholm

13743 Reputation

22 Badges

20 years, 340 days

MaplePrimes Activity


These are replies submitted by Preben Alsholm

@surnizai The range was chosen to include the result of fsolve(x*y1,x) which I actually did before I plotted.

From the plot I could see that there was (at least) one more zero and that one was near x=1+I*epsilon. So giving x=1 as a start to fsolve might get us the second roots. It did!

@patient I noticed that you changed both of the equations so you have to think again. I shall be busy with travelling the next few days so won't have time to look at it.

@matmxhu I tried the following. When doing it several times the error was successfully caught some of the time and sometimes not:

restart;
g:=proc() solve(exp(x)=1,x) end proc:
try
 timelimit(0.02, g());
catch "time expired":
 123456789
end try;

I don't know why the error is not caught all the time.
This was done in Maple 2015.
I tried in Maple 15 several times. The error was caught every time I tried.
In Maple 16, 17, and 18 the error was caught sometimes and sometimes not.

With the following code the error was caught all the times I tried in Maple 15, 16, 17, 18, and in 2015.

restart;
f:=proc() local i; for i to 10000 do evalf(exp(i)) end do end proc:
try
 timelimit(0.02, f());
catch "time expired":
 123456789
end try;

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

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