Preben Alsholm

13471 Reputation

22 Badges

20 years, 261 days

MaplePrimes Activity


These are replies submitted by Preben Alsholm

@Carl Love In the help page dsolve, numeric, DAE we find the statement:
"The following options are also available for some or all of the DAE methods:
....
implicit = boolean
... "
I tried the first example given on the help page. It uses the methods rkf45_dae (the default), rosenbrock_dae, and mebdfi.
For the first two of these the option implicit=true was accepted and gave slightly different results than without thereby indicating that a different computation was made. With method = mebdfi and implicit=true the error message was

Error, (in dsolve/numeric) 'implicit' can only be used with the rkf45, ck45, rosenbrock, dverk78, lsode, classical or gear methods

 

@sherek The imaginary unit in Maple is by default named I (capital). You are using the letter i instead.
If you try
 

interface(imaginaryunit);

and the answer is I, then you have the default setting.
You can change to another name such as j if you wish (as some people in electical engineering do).
Just do
 

interface(imaginaryunit=j);

The response this time is the previous setting I. Don't let that confuse you.

As always an example of this error happening will help you get an answer.

It is easy enough to produce the error if you feed simplify with some nonsense as in this example:
 

simplify(xxx,a=yyy);

where xxx, yyy, and a are unassigned.
But I may assume you are not talking about nonsense like that?

@Christopher2222 As we agree (I think) the ball stops bouncing in a finite time.
When I said that if we don't stop at some finite number of hits, then dsolve/numeric will do it for us I was referring to that error/warning message. Also there is roundoff/truncation or whatever to take into account.
That the integration stops because of Event2 could be considered a more graceful stop than a stop due to maxfun exceeded, but it is probably not worth the trouble.

If we try 1000 hits and raise maxfun nothing much is obtained. Try e.g. maxfun=10^6. Then Event2 is triggered at 
t = 5.19806251391305  (try dsol(5.2)  ) .
Not a whole lot different from the setting of 100 hits, where Event2 is triggered at
t = 5.19790291862644. (again try dsol(5.2) ).
Note: If we try lowering abserr and relerr we easily run into warnings about a singularity.

@Christopher2222 You have to stop the infinitely many events from happening or dsolve/numeric will do it for you.
After that it becomes another ball game, which is that of a rolling ball constrained to move on a given surface.
This modified version of Rouben's worksheet uses a discrete variable b(t) to count the number of hits and stops after 100 hits.
The surface is just a parabola.
MaplePrimes18-04-21_bouncing-ball-on-hilly-terrain.mw

 

@Rouben Rostamian  Does the total energy decrease in your revised model?
Shouldn't it?
OK, I see that it appears so:

odeplot(dsol, [t,(diff(x(t),t)^2+diff(y(t),t)^2)/2+y(t)*9.81],0..10);


Your problem has already been pointed out by acer:  e is just a name with no value, unless it has been assigned to.
In all but very old versions of Maple a warning message comes up if you try dsolve/numeric on an initial value problem with unassigned parameters. The warning from

res:=dsolve({D(y)(x)=(-2*e^x)/(2+e^x), y(0)=1/9}, type=numeric);

is
Warning, The use of global variables in numerical ODE problems is deprecated, and will be removed in a future release. Use the 'parameters' argument instead (see ?dsolve,numeric,parameters)

at least from Maple 12 and up (and including the present release Maple 2018).

I see that you have Maple V (release number irrelevant)  I checked in Maple 8 and there was no warning message either.

I tried the code above in Maple 8 and in Maple 2018, and after that made the assignment

 

e:=evalf(exp(1)); # getting e := 2.718281828

Then tried:
 

res(1.2345);
res(-1.2345);
plots[odeplot](res,[x,y(x)],x=-1..2);

All 3 lines worked in Maple 12 and Maple 2018 (but the method certainly cannot be recommended).
The first line worked in Maple 8 and so presumably also in your version.
The other two lines didn't work.
Don't mess with this: Use exp(x).
##Note: I recall that in old versions (older than Maple 8) the name E (capital) stood for exp(1).
I cannot recall when that stopped.
Try in your Maple:

evalf(E);

if 2.7128 ... comes up you have such a version. If E is the answer, you don't.

@Carl Love I checked the following trivial pdes in Maple 12, 15, 16, 17, and 18:
 

eq1:=diff(u(x,t),x)=0;
eq2:=diff(v(x,t),t)=0;
dsolve({eq1,eq2},{u(x,t),v(x,t)});

Maple 12, 15, and 16 responded with
Error, (in dsolve) not an ODE system, please try pdsolve
whereas the more recent ones acted as pdsolve would.

I prefer keeping things separate. It helps us think what we are actually doing. It is part of clarifying the problem before starting to solve it.

Thus I am not pleased with the command PDEtools:-Solve either. It is described this way:

"The Solve command computes the value of solving_variables that solves a system of equations sys. The system being solved can involve algebraic or differential equations, or both. You can request an exact (default), numeric, or series solution (respectively use the option numeric or series). In this sense, Solve is a unified command that understands when to call solve, fsolve, dsolve, or pdsolve according to your input, also facilitating the analysis of different types of solutions by just adding the keywords series or numeric." [my emphasis]

This makes me think about my response to students, who faced with some problem would automatically try solve.
My response: solve( "What is the meaning of life?" );

 

Your ball stops at a finite time after an infinite amount of jumps, doesn't it?
Some quick calculation suggests that it stops at t = 3*sqrt(2), which is approximately 4.242640686 thus it may not be so surprising that you got problems. But you were probably aware of that.
## So dsolve/numeric will be expected to handle an infinity of events!  :-)

This brings up the often lamented fact that while the help page for dsolve,numeric,events mentions quite a lot of possibilities, there are only two examples given on that page.
Those two examples don't illustrate (and cannot be expected to illustrate) all the many features apparently available according to the help page. Thus this presumably great work is basically out of reach to the average user.
Documentation is badly needed!

@tomleslie An exact solution (as in y(x) = sin(x)*exp(x) ) is nice to have if the functions appearing there have well documented properties (as do sin and exp). It is an additional advantage if those functions are also well known to the actual user.
An exact solution as in our case involving AiryAi, AiryBi, and unevaluated integrals involving those functions surely could find some use in some situation. In the present case, however, where the only important thing is to grind out numbers for a comparison with an approximate solution (found by using a wavelet approach) I would strongly favor using dsolve/numeric as I have already stated.
After all, the uniquely given maximal solution to the initial value problem

ode:=diff(Y(t),t,t)-diff(Y(t),t)+t*Y(t)=t^2;
ics:= {Y(0)=1,D(Y)(0)=-1});

is defined on all of the real axis and may be given a fancy name, why not Fellini, and the recipe for computing that exact solution Y(t) = Fellini(t) is simple: Use dsolve/numeric.

Even the help page for AiryAi and AiryBi starts by referring to the basic property of those functions:

"The Airy wave functions AiryAi and AiryBi are linearly independent solutions for w in the equation w'' - z*w = 0."

 

@tomleslie You say " Maple can solve this ODE exactly, so I have used this as the "exact" solution ".
OK, in some sense yes. But the solution is given in terms of integrals that Maple cannot find exactly:
 

ode:=diff(Y(t),t,t)-diff(Y(t),t)+t*Y(t)=t^2;
sol:=dsolve({ode,Y(0)=1,D(Y)(0)=-1});
value(sol);

Y(t) = -(1/2)*exp((1/2)*t)*AiryAi(1/4-t)*(3*AiryBi(1/4)-2*AiryBi(1, 1/4))/(AiryAi(1/4)*AiryBi(1, 1/4)-AiryAi(1, 1/4)*AiryBi(1/4))+(1/2)*exp((1/2)*t)*AiryBi(1/4-t)*(3*AiryAi(1/4)-2*AiryAi(1, 1/4))/(AiryAi(1/4)*AiryBi(1, 1/4)-AiryAi(1, 1/4)*AiryBi(1/4))+exp((1/2)*t)*Pi*(AiryAi(1/4-t)*(int(AiryBi(1/4-_z1)*_z1^2*exp(-(1/2)*_z1), _z1 = 0 .. t))-AiryBi(1/4-t)*(int(AiryAi(1/4-_z1)*_z1^2*exp(-(1/2)*_z1), _z1 = 0 .. t)))

Thus these integrals will have to be computed numerically. The advantage of the "exact solution" is therefore hard to see.

This is not an answer, but just a comment about the production of the data from the worksheet.
Your worksheet doesn't give the data, but it can be generated like this after your worksheet is run:
 

interface(rtablesize=16); #To see some numbers (default is 10)
resM:=dsolve({ode,Y(0)=1,D(Y)(0)=-1},numeric,abserr=1e-15,relerr=1e-13,output=Array([seq((2*i-1)/32,i=1..16)]));
A:=resM[2,1]:
yap:=y~(A[..,1]);
err:=A[..,2]-yap;
M:=<yap|A[..,2]|err>;

Notice that the author of the paper from which your image is taken uses a fifth order polynomial as representing the exact solution. I'm using the result from dsolve/numeric which represents the exact solution much better. I added a comment about that in your other question:
https://www.mapleprimes.com/questions/224390-The-Code-Doesn-T-Work-Why#comment247856

@student_md I wouldn't be the right one to answer that one.
But my guess is that the best way is to export the numbers as a matrix to a text file using either Export or ExportMatrix.
Then do the latex work outside of Maple.
But feel free to ask a new question in this forum about that general problem.
######################################
Note: It appears that the author of the paper uses a truncated series solution using highest degree 5.
That is not particularly accurate and can hardly be claimed to represent the exact solution.
You are much better off using the numerical solution as representing the exact solution.
You could of course increase the order as in this example:
 

resS:=dsolve({ode,Y(0)=1,D(Y)(0)=-1},Y(t),series,order=13);
yS:=convert(rhs(resS),polynom);

But yS is still an approximation and so is the numerical solution, but I suggest using the latter because of the error control available. With the default setting of Digits (10) you can even use abserr=1e-15, relerr=1e-13 as in this:

res:=dsolve({ode,Y(0)=1,D(Y)(0)=-1},numeric,abserr=1e-15,relerr=1e-13);

 

@student_md I have cleaned the code up some. Removed t as formal parameter from hi and pn since t is being used globally in other places. One thing I haven't done is to replace the indefinite integrals int(hd[i],t) and int(pn(i,n-1),t) by definite ones. You should be aware of that. Are you sure that the integral returned is the one you want?
 

Reference:
http://www.m-hikari.com/ams/ams-2012/ams-125-128-2012/sunmonuAMS125-128-2012.pdf

restart;
h1 := t-> piecewise(0 <= t and t < 1, 1):
###### We will be using t as the global time parameter throughout. #####
hi:=proc(j,k) local a,b,c,m; 
  m:=2^(j);
  a := k/m; 
  b := (k+1/2)/m; 
  c := (k+1)/m;
  return piecewise(a <= t and t < b, 1, b <= t and t < c, -1)
end proc:
##
J:=3: 
N :=2^J:
hd := Vector(N): 
H := Matrix(N, N): 
T := Vector(N):
hd[1] := h1(t):
for i from 1 to N do T[i] := (i-1/2)/N end do:
##
for j from 0 to J-1 do
  m := 2^j;
  for k from 0 to m-1 do
    i := m+k+1;
    hd[i] := hi(j, k)
  end do
end do:
#Now Compute H at the collocation points
for i from 1 to N do
  for j from 1 to N do
    H[i, j] := eval(hd[i], t = T[j])
  end do
end do:
##
pn:=proc(i,n) 
  if n=1 then 
    return int(hd[i],t) ## Notice the use of an indefinite integral!
  else
    return int(pn(i,n-1),t) ## Notice the use of an indefinite integral!
  end if
end proc:
##
##EXAMPLE 2 from the reference given:
f:=t->t^2;
alpha1:=t->-1:
alpha2:=t->t:
y0:=1:
y1:=-1:
##
RHS:= t->f(t)-alpha1(t)*y1-alpha2(t)*(t*y1+y0);
R:=Vector(N):
TMP:=Matrix(N,N):
A:=Matrix(N,N):
##
for i from 1 to N do
  R[i] := evalf(RHS(T[i]));
  tmp := alpha1(t)*pn(i,1)+alpha2(t)*pn(i,2);
  for j from 1 to N do
    TMP[i,j]:=eval(tmp, t = T[j])
  end do
end do:

##Compute the wavelet coefficients:
use LinearAlgebra in
  A := Transpose(LinearSolve(Transpose(H+TMP), R))
end use:
#Now compute the approximate solution
sol := y1*t+y0+add(A[m0]*pn(m0,2), m0 = 1 .. N):
#Convert to a function of t for easy comparison with exact solution
y:=unapply(sol,t):
#############################################
## Compare with the solution found by dsolve/numeric:
## EXAMPLE 2 ode:
ode:=diff(Y(t),t,t)-diff(Y(t),t)+t*Y(t)=t^2; 
res:=dsolve({ode,Y(0)=1,D(Y)(0)=-1},numeric,abserr=1e-10,relerr=1e-9);
plots:-odeplot(res,[t,Y(t)-y(t)],0..1,caption="The error on the interval 0..1");

########## In the text in the reference I see that integrals from 0 to t are used.
In recent versions of Maple (certainly from 12 and up, but not in Maple 8) you can use the physicists' abuse of notation t=0..t: Example:

u:=exp(t):
int(u,t); # not zero at t=0.
int(u,t=0..t); # t is used in two roles.
## Alternatively, without abuse of notation and works in all versions:
subs(tt=t,int(u,t=0..tt));


 

First 39 40 41 42 43 44 45 Last Page 41 of 225