Preben Alsholm

13728 Reputation

22 Badges

20 years, 254 days

MaplePrimes Activity


These are replies submitted by Preben Alsholm

@yemisi I know nothing about Block BDF.
But certainly it should be easy to modify Rouben Rostamian's Euler code to this new case. The major difference is that y(t) is not constant for t < 0, but equal to sin(t).
Something like this modification should work:

restart;
doit := proc(h, T)
  local n, y, i, t, y_arg, k;
  n := ceil(T/h);   # number of steps
  #Defining y for t <=0:
  for i from 0 to -ceil(evalf(Pi/2)/h) by -1 do y[i] := evalf(sin(i*h)) end do;
  for i from 1 to n do
    t := i*h;
    y_arg := t-evalf(Pi/2);
    # ... we convert that to a grid index  
    k := round(y_arg/h);
    # Euler's method:
    y[i] := y[i-1] - h*y[k]
  end do;
  return [seq([i*h, y[i]], i=0..n)];
end proc:
#Trying with stepsize h=0.01:
res := doit(0.01, 10):
plot(res, style=point, color=blue, symbol=point, labels=[t,y]); p1:=%:
plots:-display(p1,plot(sin(t),t=0..10));
############
Rouben Rostamian's code can be found here:
http://mapleprimes.com/questions/204395-DELAY-DIFFERENTIAL-EQUATIONS

You got plenty of response to your first problem. This problem has constant delay so is much simpler. In fact it can be solved symbolically.
Incidentally (but you probably know already) the solution is y(x) = sin(x) for all x.

######
When is your homework due?

@Markiyan Hirnyk When  I said my ddesolve I meant that it is my private software. The name ddesolve is surely not owned by anyone.

@Rouben Rostamian  I tried eq on my own ddesolve, which I have had for a while. I get a graph that is the same as Maple's dsolve up to 0.84 after which my code won't work (reasonably enough).
So what is going on in Maple for x>0.84 when the delay becomes negative?
One guess is that there is extrapolation going on, as there could be for x < 0. Try plotting for x < 0 where certainly  dsolve's algorithm for computing values for x > 0 assumes that y(x) = 1.
The blue curve is done with my ddesolve, the red with Maple's dsolve. Visual agreement for x > 0 (x<0.84), but a striking difference for x < 0.

@Carl Love While examining the beta-version of dsolve/numeric/delay which was in fact available in Maple 18, I found a problem with output=listprocedure, which makes me suspicious when that (or output=operator) is involved.

restart;
dde:=diff(y(t),t)=-y(t-(1+cos(t)));
res:=dsolve({dde,y(0)=1},numeric,output=listprocedure,delaymax=2);
Y:=subs(res,y(t));
plots:-odeplot(res,[t,y(t)],0..10); #OK
plot(Y,0..10); #OK
X:=Vector([seq(0..10,.1)]);
Y~(X);
#This is still broken in Maple 2015.1:
plot(X,Y~(X));
plot(Y,0..10);
## But this is fine:
plots:-odeplot(res,[t,y(t)],0..10);
##############################
###########
The image you show is basically the same in Maple 2015.1.
One wonders why the difference in the beginning is not roughly 0.
Clearly the understanding in dsolve/delay is that y(0) = 1 means that y(x) = 1 for x<=0. Comparing results from Maple's dsolve with other methods clearly shows that to be true. Plotting for x < 0 should not confuse you, but is likely to:
plots:-odeplot(res,[x,y(x)],-1..1);
#Maybe extrapolation is used?


@Rouben Rostamian  I cannot tell you,but surely it is a good question.

@Markiyan Hirnyk You have a totally different equation, but OK.
Now the "delay" is given by t-tau = y(t)+1, thus tau = t-y(t)-1. Therefore tau < 0 at least for t=0, where tau = 0-0.5-1 = -1.5 < 0, i.e. tau is NOT a delay.
The facility for doing delay differential equations was written with the understanding that the delay must be just that, a delay.

@Markiyan Hirnyk Since the delay tau is given by y(x)-2=x-tau, we have (as also stated in the question posed) tau = x-y(x)+2, which certainly will be considered a variable delay.

@Carl Love Considering that the presumably much simpler system odesys3 obtained from odesys by replacing exp by 1, the sigmas by 0, and all floats by their signs is difficult (if not impossible?) I would think it is a waste of time to let this run for days:

odesys2:=eval({odesys},{:-exp=1,sigma[1,1]=0,sigma[1,2]=0});
odesys3:=evalindets(odesys2,float,signum);
:-dsolve(odesys3);

The systems are quadratic in the y's. Even for a quadratic system in 2 dependent variables solution is not trivial.



@acer It may boil down to automatic simplification in one case and not another:
restart;
0.1*(0.4*eta+1); #Automatically simplified
0.1*(0.4*eta+1)*(2-x); #Not so
0.1*(0.4*eta+1)*(2-x)/(0.1*(0.4*eta+1)); #Consequently no cancellation

#Mimicking to some degree what is going on in isolate which uses `isolate/isolate1`:
restart;
eq:=.4*eta-0.8e-1+(.1*(.4*eta+1))*(2-x)=0;
EQ:=subs(x=_XX,eq);
EQ1:=EQ:
to 10 do EQ1:=`isolate/isolate1`(EQ1) end do;
isolate(EQ,_XX,10); #10 iterations instead of 100000


@acer A strange coincidence. I sent an SCR giving the very same example yesterday.

@emmantop The way I usually do it is to use debug (aka trace)as in
debug(`dsolve/numeric/bvp/convertsys`);
#To stop debugging use undebug:
undebug(`dsolve/numeric/bvp/convertsys`);
#Also taking a look at the procedure is a good idea. Here one can use showstat:
showstat(`dsolve/numeric/bvp/convertsys`);

These facilities (debug and showstat) have been around forever.

##############
Having found the problem: that it lies with the use of isolate (instead of solve), it is tempting to try to replace `dsn/solve` with a version where isolate is replaced with solve:

`dsn/solve`:=subs(isolate=(proc(eq,var) solve(eq,{var})[1] end proc),eval(`dsn/solve`));

After that `dsolve/numeric/bvp/convertsys` will work fine on the problem we have, and so dsolve(....) should work too.
However, I think it is wiser to stick with one of the workarounds.

It is interesting that the help page for isolate says:
"For direct solutions of equations, solve may be more efficient than isolate, which is intended, primarily, for use in an interactive Maple session".
But isolate has the advantage of being a much less complicated procedure than solve.


@Alejandro Jakubi If we assume that x>0 it suffices to use Parts:

restart;
with(IntegrationTools):
J:=Int(ln(1+sqrt((x+1)/x)),x):
J1:=Parts(J,ln(1+sqrt((x+1)/x)));
simplify(J1) assuming x>0;
value(%);


@jjcp You are not doing anything wrong, but if you only carry 4 digits in a complex computation it is not surprising  ( or rather it should be expected) that the result obtained is not the same as the one obtained by rounding to 4 digits a result computed carrying 10 digits.

Could you explain why you are asking?

After all what does it mean to set Digits to 4 or 10 (the default)? Surely you know of rounding?

In case the point is that not even the first 4 digits in the first result are all correct then remember that setting Digits to 4 does not guarantee 4 correct significant digits in the result.

Trivial example:
Digits:=10:
evalf(Pi)-3.14159;
evalf(Pi-3.14159); #Same result, but the point is that there are only 4 significant digits, not 10


First 115 116 117 118 119 120 121 Last Page 117 of 230