pagan

5142 Reputation

23 Badges

16 years, 258 days

 

 

"A map that tried to pin down a sheep trail was just credible,

 but it was an optimistic map that tried to fix a the path made by the wind,

 or a path made across the grass by the shadow of flying birds."

                                                                 - _A Walk through H_, Peter Greenaway

 

MaplePrimes Activity


These are replies submitted by pagan

I remembered that someone else asked this a while back, but it took me a while to find it.

In terms of the Help system (apart from the manuals) it's a little too difficult to uncover this stuff. The Help page for topic colondash isn't cross-referenced as much as it could be. The content and examples on that page are also mostly about package members' names. It'd be improved with another example explaining the subject of this thread, eg. for option keywords of procedures.

Is this in 16.01? If still problematic for that point-release then sending it to Maplesoft tech support might help to get a fix for some future release.

Good going, for figuring it out. But what meaning can we give to normalizing in just one direction (IFFT but not FFT)?

Good going, for figuring it out. But what meaning can we give to normalizing in just one direction (IFFT but not FFT)?

@PatrickT You've used StringTools:-Has backwards, I think. After redefining WARNING in Preben's fashion, and running your example:

> StringTools:-Has(lastwarning,"event");
                                     true


> StringTools:-Has(lastwarning,"cannot");
                                     true

BTW, you could emit your own call to `error`, instead of ln(0), to force an exception. That'd also let you customize the error string for it too (in case you ever want to catch it specifically).

But I would agree that better would be to use the event halt mechanism, and trying to clear events between tries.

@PatrickT You've used StringTools:-Has backwards, I think. After redefining WARNING in Preben's fashion, and running your example:

> StringTools:-Has(lastwarning,"event");
                                     true


> StringTools:-Has(lastwarning,"cannot");
                                     true

BTW, you could emit your own call to `error`, instead of ln(0), to force an exception. That'd also let you customize the error string for it too (in case you ever want to catch it specifically).

But I would agree that better would be to use the event halt mechanism, and trying to clear events between tries.

@xiao2788 You are not getting a result because the call to fsolve is failing. It is failing because the independent data values (in cx) are not unique.

Your new `cx` and `cy` lists do not describe a function, because the same value in cx can correspond to multiple values in cy. At best, you might decide that this new data represents a discontinuous function (piecewise continuous, say).

But in that case, how will you join up the paired cx and cy points? Your original post said that the x and y data in the lists may not be ordered. Does that mean that you are supposed to re-order cx and cy (as well as ax and ay)? Your post was unclear there.

You have several problems, then. The first is that you will have to decide whether cx and cy can be re-ordered, or must be taken in the order given (for this new example, and presumably any other examples to follow...). The next problem is that you'll have to decide whether to create a piecewise continuous curve. After that, you'll have to decide whether you want linear or some higher order interpolation of the point values. (How do you intend on matching derivatives, if at all, between pieces?)

In summary, there are enough ambiguities in your new problem, and you have not specified well enough what you want. Sure, someone could give you some answer, but as things stand it would be somewhat arbitrary in nature.

@xiao2788 You are not getting a result because the call to fsolve is failing. It is failing because the independent data values (in cx) are not unique.

Your new `cx` and `cy` lists do not describe a function, because the same value in cx can correspond to multiple values in cy. At best, you might decide that this new data represents a discontinuous function (piecewise continuous, say).

But in that case, how will you join up the paired cx and cy points? Your original post said that the x and y data in the lists may not be ordered. Does that mean that you are supposed to re-order cx and cy (as well as ax and ay)? Your post was unclear there.

You have several problems, then. The first is that you will have to decide whether cx and cy can be re-ordered, or must be taken in the order given (for this new example, and presumably any other examples to follow...). The next problem is that you'll have to decide whether to create a piecewise continuous curve. After that, you'll have to decide whether you want linear or some higher order interpolation of the point values. (How do you intend on matching derivatives, if at all, between pieces?)

In summary, there are enough ambiguities in your new problem, and you have not specified well enough what you want. Sure, someone could give you some answer, but as things stand it would be somewhat arbitrary in nature.

@PatrickT The matching attempt of the strings is done... from the left only.

That's why you have sucess when you use  "cannot evaluate the solution" as the catchstring. But "singularity" or "probably a singularity" occur toward the right. And so are not matched. The help is not explicit about this, but I have only ever seen behaviour that can be explained by the rule: it must match, from the left end of the string.

Having said that:

> restart;

> sol:=dsolve({diff(f(x),x)=1/x, f(-1)=1/10},numeric,output=listprocedure):

> sol(1);
Error, (in unknown) cannot evaluate the solution further right of -0.10619384e-304,
     probably a singularity

> lastexception;
  unknown, "cannot evaluate the solution further right of %1, 
                                           -305
     probably a singularity", -1.0619384 10    

So now try that using string shown by lastexception, including the placeholder %1.

> g:=proc()
>   try
>      sol(1);
>   catch "cannot evaluate the solution further right of %1, probably a singularity":
>   end try;
> end proc:

> g(); # caught!

@PatrickT The matching attempt of the strings is done... from the left only.

That's why you have sucess when you use  "cannot evaluate the solution" as the catchstring. But "singularity" or "probably a singularity" occur toward the right. And so are not matched. The help is not explicit about this, but I have only ever seen behaviour that can be explained by the rule: it must match, from the left end of the string.

Having said that:

> restart;

> sol:=dsolve({diff(f(x),x)=1/x, f(-1)=1/10},numeric,output=listprocedure):

> sol(1);
Error, (in unknown) cannot evaluate the solution further right of -0.10619384e-304,
     probably a singularity

> lastexception;
  unknown, "cannot evaluate the solution further right of %1, 
                                           -305
     probably a singularity", -1.0619384 10    

So now try that using string shown by lastexception, including the placeholder %1.

> g:=proc()
>   try
>      sol(1);
>   catch "cannot evaluate the solution further right of %1, probably a singularity":
>   end try;
> end proc:

> g(); # caught!

@xiao2788 

1) "axycurve := x -> CurveFitting:-ArrayInterpolation(axy, x, method = spline, extrapolate) " produces a procedure which accepts an x-point returns the corresponding y-point on the interpolated curve made from the ax and ay data. So, axycurve(1.41) will return the y-value on that curve which corresponds to x=1.41

 

2) The curves axycurve and cxycurve appear to intersect in at least two places. One of those places is clearly seen to lie within the range of the x-data. The second place of intersection might be slightly off to the right, not lying in the x-range of the data, and only seen by extrapolating the curves out a tiny bit more to the right.

You are right, in general the intersection of two curves produced from data in this way will not necessarily fall within the subranges min(ax)..max(ax) and min(cx)..max(cx) in particular. But for your particular data,,, they do. I observed that from the plot. You could also do it this way:


fsolve(axycurve-cxycurve,min(min(ax),min(cx))..max(max(ax),max(cx)));

                          17944.44743

fsolve('axycurve'(t)-'cxycurve'(t),
       t=min(min(ax),min(cx))..infinity,avoid={t=%});

                          25199.98161

There are other root-finders which likely could accomplish what fsolve is doing here, eg. NextZero, DirectSearch. The root-finding isn't so tricky. The key is to get the interpolated curves (be they linear, or spline, etc) whose difference can be solved for points of intersection. I expect that BSplineCurve could be another way to generate the interpolated curves (as piecewise expressions, rather than as procedures).

 

3)

axycurve(17944.44743);

                        0.21786664019881208

cxycurve(17944.44743);

                        0.21786664122100774

@xiao2788 

1) "axycurve := x -> CurveFitting:-ArrayInterpolation(axy, x, method = spline, extrapolate) " produces a procedure which accepts an x-point returns the corresponding y-point on the interpolated curve made from the ax and ay data. So, axycurve(1.41) will return the y-value on that curve which corresponds to x=1.41

 

2) The curves axycurve and cxycurve appear to intersect in at least two places. One of those places is clearly seen to lie within the range of the x-data. The second place of intersection might be slightly off to the right, not lying in the x-range of the data, and only seen by extrapolating the curves out a tiny bit more to the right.

You are right, in general the intersection of two curves produced from data in this way will not necessarily fall within the subranges min(ax)..max(ax) and min(cx)..max(cx) in particular. But for your particular data,,, they do. I observed that from the plot. You could also do it this way:


fsolve(axycurve-cxycurve,min(min(ax),min(cx))..max(max(ax),max(cx)));

                          17944.44743

fsolve('axycurve'(t)-'cxycurve'(t),
       t=min(min(ax),min(cx))..infinity,avoid={t=%});

                          25199.98161

There are other root-finders which likely could accomplish what fsolve is doing here, eg. NextZero, DirectSearch. The root-finding isn't so tricky. The key is to get the interpolated curves (be they linear, or spline, etc) whose difference can be solved for points of intersection. I expect that BSplineCurve could be another way to generate the interpolated curves (as piecewise expressions, rather than as procedures).

 

3)

axycurve(17944.44743);

                        0.21786664019881208

cxycurve(17944.44743);

                        0.21786664122100774

A nice bug. It has some characteristics of the issues that can occur when Maple generates an expression sequence with a single member. (Not really supposed to happen, I think, but sometimes gets through the cracks...)

This is 64bit 15.01 on Windows.

restart:

Test:=Array(1..3,[1,2,3]):

member(3,Test,'l');
                              true

l;
                               3

lprint(l);
3

dismantle(l);

INTPOS(2): 3


addressof(l), addressof(3);
                              7, 7

l-3;
                               0

L:=3:

addressof(l), addressof(L);
                              7, 7

if L < 100 then hey else whoa end if;
                              hey

if l < 100 then hey else whoa end if;
Error, cannot determine if this expression is true or false: (3) < 100

disassemble(l);
                              2, 1

disassemble(L);
                              2, 1

l[1]; # ahah
                               3

L[1];
                              3[1]

if l[1] < 100 then ok else nok end if;
                              ok

if op(1,l) < 100 then ok else nok end if;
                              ok

if op(1,L) < 100 then ok else nok end if; # form that works also for usual posints, fwiw
                              ok

A nice bug. It has some characteristics of the issues that can occur when Maple generates an expression sequence with a single member. (Not really supposed to happen, I think, but sometimes gets through the cracks...)

This is 64bit 15.01 on Windows.

restart:

Test:=Array(1..3,[1,2,3]):

member(3,Test,'l');
                              true

l;
                               3

lprint(l);
3

dismantle(l);

INTPOS(2): 3


addressof(l), addressof(3);
                              7, 7

l-3;
                               0

L:=3:

addressof(l), addressof(L);
                              7, 7

if L < 100 then hey else whoa end if;
                              hey

if l < 100 then hey else whoa end if;
Error, cannot determine if this expression is true or false: (3) < 100

disassemble(l);
                              2, 1

disassemble(L);
                              2, 1

l[1]; # ahah
                               3

L[1];
                              3[1]

if l[1] < 100 then ok else nok end if;
                              ok

if op(1,l) < 100 then ok else nok end if;
                              ok

if op(1,L) < 100 then ok else nok end if; # form that works also for usual posints, fwiw
                              ok

@PatrickT I suspect some misunderstanding.

Firstly, in your earlier versions you made assignments such as sol:=dsol along with the comment that this produced a copy of dsol. But that is not so. That assignment does not make a copy of dsol, or of the three procs within dsol. (More on copying, below...)

Even if you assign the value of (read: the procedure assigned to) dsol to various other names, you still have to reset the initial conditions (over again) of the earlier names each time that you wish to go back and recalculate with them.

I don't see any cheap way around having to reset the initial conditions each time you wish to recalculate. You don't have to reinvoke dsolve, I think, but you do have to reinstate the various initial conditions. For example,

restart:

dsol :=
  dsolve(
      [ diff(u(t),t)=u(t)-v(t), diff(v(t),t)=u(t)-w(t), diff(w(t),t)=v(t)-w(t)
      , u(0)=2, v(0)=1, w(0)=2 ]
      , [u(t),v(t),w(t)]
    , 'type' = numeric
    , 'output' = listprocedure
  ):

solA:=dsol:
solA(initial=[2,1,2]):
solA(-135.0);

[t(-135.0) = -135.0, u(t)(-135.0) = HFloat(3.9077543506678216), 

  v(t)(-135.0) = HFloat(4.992270555376236), 

  w(t)(-135.0) = HFloat(4.0845162047084)]

solB:=dsol:
solB(initial=[6,7,8]):
solB(-135.0);

 [t(-135.0) = -135.0, u(t)(-135.0) = HFloat(8.08460488495403), 

   v(t)(-135.0) = HFloat(7.17680411000912), 

   w(t)(-135.0) = HFloat(6.092199225055099)]

solA(initial=[2,1,2]):
solA(-135.0);

[t(-135.0) = -135.0, u(t)(-135.0) = HFloat(3.9077543506678216), 

  v(t)(-135.0) = HFloat(4.992270555376236), 

  w(t)(-135.0) = HFloat(4.0845162047084)]

Are you seeing behaviour where this kind of approach is not good enough, and still has "contamination"?

In light of the above, there's not much point to having different names solA, solB, etc, with this particular scheme above. But... see at end for an alternative which does merit the various names.

I don't think that you can get around this calling the `copy` command on each of the three procs within dsol either. I mean, by doing solAu:=copy(eval(u(t),dsol)), etc, or with various judiciously placed calls like forget(solB), etc. And that's not only because of the way Maple uniquifies procs (ie. "simpls them"). You can subsop(1=..,eval(u(t),dsol)) and get various instances with actual different addresses, and still have same behaviour.

If you want to have various named procedure, for convenience of reexecution, then you can write them as simple procedures which wrap around preliminary calls to reset the initial conditions. Eg.

solA:=proc(tt) dsol(initial=[2,1,2]); dsol(tt); end proc:

solB:=proc(tt) dsol(initial=[6,7,8]); dsol(tt); end proc:

solA(-135.0);

[t(-135.0) = -135.0, u(t)(-135.0) = HFloat(3.9077543506678216), 

  v(t)(-135.0) = HFloat(4.992270555376236), 

  w(t)(-135.0) = HFloat(4.0845162047084)]

solB(-135.0);

 [t(-135.0) = -135.0, u(t)(-135.0) = HFloat(8.08460488495403), 

   v(t)(-135.0) = HFloat(7.17680411000912), 

   w(t)(-135.0) = HFloat(6.092199225055099)]

solA(-135.0);

[t(-135.0) = -135.0, u(t)(-135.0) = HFloat(3.9077543506678216), 

  v(t)(-135.0) = HFloat(4.992270555376236), 

  w(t)(-135.0) = HFloat(4.0845162047084)]

Hopefully this may be of some use to you.

2 3 4 5 6 7 8 Last Page 4 of 81