Preben Alsholm

13743 Reputation

22 Badges

21 years, 118 days

MaplePrimes Activity


These are replies submitted by Preben Alsholm

I fail to see the connection between the original ode and the rest of the code.

Your first order system would be
y'(t) = u(t)
u'(t) = -3*u(t)-2*y(t) + cos(t)

Just one indication of the missing connection: Where is cos(t) in the code?

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


First 122 123 124 125 126 127 128 Last Page 124 of 232