pagan

5057 Reputation

23 Badges

14 years, 60 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 answers submitted by pagan

You might create it as TypeMK, which is a MathML analogue that Maple's Standard GUI uses to represent 2D display Math.

This can usually be done by building up the expression as 2D Math input, and then before executing it (to get 2D Math output) using the right-click context menu action 2D-Math -> Convert To -> Atomic Identifier. Then execute. Then lprint(%) to get it as 1D Maple notation input code.

For example, this input (as one long line, unbroken)

`#mrow(munder(mo("lim"),mrow(mi("x"),mo("→"),msup(mn("3"),mo("+")))),
mfrac(mrow(mo("&uminus0;"),mn("4")),mrow(mi("x"),mo("−"),mn("3")),linethickness = "1"))`

displays as output like

The above is not a programmed approach, and in general requires interaction. So it might not serve for code in that Button component, unless you just want variations on this particular simple example in which case you can build it up using the `cat` command.

Alternative methods can be messy. It can be hard to convince Maple not to pull the negative out front when the numerator is just an integer. You can fake it by converting the integer to a name (in which case the font for it changes from upright to italic). And sometimes you can glue bits together using operators like ``() which just happen to get hidden during 2D output printing. These can get awkward, though opinions might vary. Also, if you don't mind changing the denominator as well, you can sometimes redistribute the negative sign below

 

with(RandomTools):

a1 := Generate(integer(range = -10 .. -1)):

T:=a1/(x-3):

Limit(T, x=3, right);

T:=(numer(T)/signum(numer(T)))/expand(denom(T)/signum(numer(T))):

Limit(T, x=3, right);

 

sol:=eliminate({x1=X1, x2=X2, y1=Y1, y2=Y2},{b,a,M,phi}):

ans:=eval(M,sol[1]):

The local `c` is a table.

You can convert it to a list, and return that. Either of these should work, at the end of that procedure.

return [entries(c,'nolist')],

return convert(c,list);

An alternative is to replace the inner `j` loop with `mul`, the inner `i` loop with `add`, and replace the outer `k` loop with `seq` to build the list directly.

I'm not sure how you want the two columns of your Matrix to come into play here. It sounds like you are trying to find these sums by adding rows, it that right?

> restart:

> M:=Matrix(4,2,fill=2);

                                      [2    2]
                                      [      ]
                                      [2    2]
                                 M := [      ]
                                      [2    2]
                                      [      ]
                                      [2    2]

> N:=copy(M):

> for i from 1 to (4-1) do
>   LinearAlgebra:-RowOperation(N,[i+1,i],1,inplace);
> end do:

> N;

                                   [2    2]
                                   [      ]
                                   [4    4]
                                   [      ]
                                   [6    6]
                                   [      ]
                                   [8    8]

Did you want the two columms treated separately, like that? Or did you want them treated as if they were one big single column (the botton of the first column being succeeded by the top of the second column), in which case the N[4,1] entry of 8 would have to be added to each entry in the second column. Or you could Alias N to an 8x1 Matrix before the operation.

If you call int(..., x=0..1) or evalf(int(...,x=0..1)) then Maple will try to do an exact symbolic analysis, and the answer to your question will relate to the methods it uses for that. I'm pretty sure others will explain this exact symbolic side of things, so I'll restrict my answer and suppose that you are asking about numeric (floating-point approximate) integration.

If you call evalf(int(...,x=0.0..1.0)) then `int` will invoke `evalf/int` to do numeric quadrature, which is also what will get used if you call evalf(Int(...,x=0..1)). How this gets done relates to the various numeric quadrature methods (eg. _NCrule, _Dexp, _CCquad, etc) and the control code which governs how they get used internally. You can set infolevel[`evalf/int`]:=6 and see quite a lot of detail. But not everything, as some of the methods can call external compiled code, or don't report everything. In general, evalf(Int(...)) tries to determine whether the numeric integral converges to within a tolerance, while it also knows that it should handle singularities specially. To make it more complicated, sometimes evalf(Int(...)) may try symbolic analysis of the integral, so getting only numeric determination of convergence may necessitate passing a procedure or operator, eg. evalf(Int(x->1/((exp(x)-1),0..1)) to prevent this.

If you are looking for a quite separate means for numeric convergence determination, you can try and convert it to a discrete numeric summation problem, and use some convergent tests for that. For example,

> restart:

> #showstat(`evalf/Sum/LevinU_Symb/1_inf`);

> #trace(`evalf/Sum/LevinU_Symb/1_inf`):

> G:=t->evalf(Int(x->1/(exp(x)-1),1/(t+1)..1/t)):

> evalf(Sum('G(T)',T=1..infinity));              

                                Float(infinity)

That Monte Carlo block is built from Embedded Components such as PlotComponents, Buttons, etc. Such components can have code underlying them, which you can view (and even change).

Right-click on the Button component labeled "Compute Option Price". In the pop-up menu that appears, select the choice "Component Properties". This in turn should pop open another pane, which includes a button labeled "Edit", which lies beside the description, "Action When Clicked". Select that Edit button. This should open a very simple code editing window, which already contains plaintext 1D Maple notation code. That code makes it work.

See here for a few ways to go about it.

For an indexed RootOf of a univariate polynomial (only indeterminate being _Z say) evalf/RootOf will call fsolve with the complex option so as to get all the roots. And then it will sort the results, and then pick the one which corresponds to the index.

So if you do this for each RootOf index=i, i=1..d for a degree d univariate polynomial, then you will make Maple compute all the d distinct roots d many times. And fsolve doesn't use remember and re-use previous results. So that is d times too much computational effort (excluding evalf's remember table's effect).

But if you extract the univariate polynomial 1st operand of the RootOf, and pass that to fsolve, then you can compute all d roots just once.

The 1..100 loop in this comparison is to help make the timings clear. Otherwise this particular very small example runs so fast that the difference is hard to detect. (Yes, I know that garbage collection is also being measured in the looped test. But garbage collection is what often happens anyway, and it's artificial to try and do timings which only ever exclude it.)

You can also test with other examples.

> restart:

> r := RootOf(189*_Z^6-756*_Z^5+567*_Z^4+819*_Z^3-1386*_Z^2+672*_Z-104):

> st:=time():
> for i from 1 to 100 do
>    fsolve(subs(_Z=x,op(r)), x);
> end do:

> time()-st;
                                    0.234

> restart:

> r := RootOf(189*_Z^6-756*_Z^5+567*_Z^4+819*_Z^3-1386*_Z^2+672*_Z-104):

> st:=time():
> for i from 1 to 100 do
>   evalf(allvalues(r));
> end do:

> time()-st;
                                    2.184

> restart:
> r := RootOf(189*_Z^6-756*_Z^5+567*_Z^4+819*_Z^3-1386*_Z^2+672*_Z-104):
> sort([fsolve(subs(_Z=x,op(r)), x)]);

   [-1.280160435, 0.3260605103, 0.6407586587, 0.8042832092, 1.033180831, 

     2.475877226]

> restart:
> r := RootOf(189*_Z^6-756*_Z^5+567*_Z^4+819*_Z^3-1386*_Z^2+672*_Z-104):
> sort([evalf(allvalues(r))]);

   [-1.280160435, 0.3260605103, 0.6407586587, 0.8042832092, 1.033180831, 

     2.475877226]

Or you can confirm the claim, by tracing fsolve to shows each different instance that it gets called.

restart:
r := RootOf(189*_Z^6-756*_Z^5+567*_Z^4+819*_Z^3-1386*_Z^2+672*_Z-104):
trace(fsolve):
Digits:=15:
for i from 1 to 1 do
   fsolve(subs(_Z=x,op(r)), x, complex);
end do:
sort([%]);

restart:
r := RootOf(189*_Z^6-756*_Z^5+567*_Z^4+819*_Z^3-1386*_Z^2+672*_Z-104):
trace(fsolve):
for i from 1 to 1 do
   evalf(allvalues(r));
end do:
sort([%]);

It sounds as if you are running the animation in the same Document in which you are coding. Is that right? I'll suppose so for now, unless you say otherwise.

You could edit the code in a separate Document from that in which you run the animation. There are several ways to do this (and maybe you're already doing one of them).

- You could have multiple Documents open, but sharing the same kernel engine. There's a GUI preference that controls such behaviour, when you open multiple sheets in the same session. I wouldn't be surprised if this also were sluggish, maybe like the scenario in which you edit and run all in a single document.

- You could have multiple Documents open, but using distinct kernel engines. In one of them you could edit the code for the procedures that do the work. (Hopefully, your code is structured and assembled into procedures, since they can be re-used and redistributed quite separate from worksheets in which they run. It's good programming practice, in most respects.) In another document you could run the animation. You could share the procedures between the documents by saving them to maple a "library". You could even do this with two completely difference instances of the whole GUI, one opened for programming and one opened for running.

- Or, as mentioned by others, you could simply use a text editor to program plaintext .mpl source files, and then `read` them into the Document that runs the animation. Or you could use a text editor, and use its facilities to savelib the .mpl source file in a library, instead using `read`.

Maybe it would be good to look for the cause of the slowness, though. Does the entire Document get slow, when you run the animation? Or is it just slow to scroll the plot on and off the visible portion? How are you animating the plot rotation? With what commands, specifically? Are you using a Plot Component in the Document? Are you animating it using the GUI or Plot Component's own animation facility, or are you mimicing that using a sequence of plots controlled by a Button Component?

Have you tried using the Standard GUI with 1D Maple notation input, inside a Worksheet instead of a Document?

How many frames do you want to be shown? This can relate to whether you use something like the plots[animate] command (which has to build and store all frames at once) or whether you mimic animation using Components (which can form only one frame at a time, and keep total allocation down).

I don't know how it is on the Mac, but on Linux and Windows the GUI and the kernel engine consume memory separately. Getting an accurate view on total memory allocation by Maple can involve finding out the separate memory use of each.

I usually see a smoother transition of the color on a surface, with less artefacts, with plotdevice=bmp.

Eg.

plotsetup(bmp): # or whichever of the ways you choose to do this

plot3d((x+y)^2,x=-12..12,y=-12..12,axes=box,lightmodel=light4,glossiness=0.2);

With plotdevice=gif or cps I see rough sections of the surface image, as the color changes splotchily.

Are you looking for less jagged lines (axes, etc) or smoother color transition, or something else?

Unless using new features like 2D Math captions, etc, I think you're likely more often better off by setting the plotdevice than you are with right-click export.

Instead of having to raise Digits, you can use evalf[d](...) instead of evalf(...)

The `d` wil be the working precision used by evalf for that call.

For example

evalf[20](subs(n = 4.24568325257496060, rhocut = 1.66125593854788400,
               m = 2.00000021862965971, rho = .956914761577446150, A));

It's important to realize that Digits and `d` above will control the working precision (sometimes commands use a few more guard digits) but they do not specify a requested accuracy for floating-point evaluation of compound expressions such as your A. There is no such accuracy control for evalf.

What you're seeing is a printed placeholder for a Matrix. The data is there but the dimension is larger than the default maximal size that gets displayed in full. You can control that.

interface(rtablesize=100):

If you issue that, then any Matrix, Vector, or Array whose dimensions are at most 100 in size will get its contents displayed in full. Anything larger gets displayed just as that kind of summary, although you can still access the entries in the usual way by using square-brackets.

As for the optimization problem, well, there may be many local minima, and the NonlinearFit command does not do multiparameter global optimization as far as I know.

In any case, couldn't you supply an upper bound for parameter b? You can practically fit just about any data as close as you want, if you allow a high enough oscillation. And if you need the amplitude to be enough to reach the highest and lowest values in list2 (or to satisfy some other constraint) then couldn't you bound a and d as well?

As mentioned, what you've been trying so far is symbolic integration since Theta is not yet numeric (and in any case it'd be evalf(Int(...)) for numeric integration).

It's not clear that the results from `simplify` are right, since the floats involved are across a wide scale, and you haven't increased Digits. As a very rough rule it's often better to not do such symbolic manipulations after you've assigned float values to some of the variables. I'd be tempted to leave the earlier integral in forming `expr1` also as inert Int, and then trying a numeric double integral.

But even supposing that `F` is OK, you've still got roundoff error to handle when plotting or integrating.

Look at the difference bewteen

plot(F(q)*sin(q)/q, q = 0 .. 50);

and

Digits:=100:
plot(F(q)*sin(q)/q, q = 0 .. 50);
Digits:=10;

So now you have this semi-infinite oscillatory integral to compute, which is something Maple is not strong at in general. And there are probably roundoff issues with the integrand, as it stands.

You could look at something like this

ans1 := Theta -> evalf(Int(q-> Re(F(q)*sin(q*Theta)/q),
                 0 .. 64*Pi, digits = 30, epsilon = 0.1e-1, method = _NoNAG));

Since you're only interested in plotting (vs Theta) then you don't need fine accuracy. But you likely do need higher working precision. (And Digits will likely have to be set even higher, as a requirement, as you push out the upper limit of `q`. You can't just keep the upper limit of 64*Pi as above -- you want the integral to infinity.)

Finally, you could get a rough view (vs Theta) with a point plot.

plot(ans1,4..5,numpoints=10,adaptive=false);

But the bigger problem of the integral's being semi-infinite and oscillatory is still there. You mentioned Fortran. How do you know that your scheme produced accurate results? What scheme was it?

I've heard hints that semi-infinite oscillatory integrals might be approached by computing the periodic segments (here, a nice constant 2*Pi/Theta in length), computing the numeric integrals for those separately, and then seeing if an acceleration technique could be done on the sequence of partial sums. Axel might know about this.

I doub't that the original inner integral of expr1 could be done purely symbolically, due to the abs of lns (but it'd be great for you if I were wrong on this).

Try this, if your F is a procedure or operator

Optimization:-Minimize(Y->evalf(Int(F(x,Y),x=0..Pi/2)), -0.1..0.1);

A slight change (using `eval`) would be needed if F is just an expression in x and y.

You might need to add `epsilon` tolerance option, or `method` option, to the evalf(Int(...)) depending on the numeric behaviour. You might also need to use a proc for the objective which returns a "large value" when the integral (objective function) fails to converge.

It's easier if you can provide the whole example.

Vdata:=Vector([0.2,0.4,0.6,1.0,1.2,1.7,1.9]):
Idata:=Vector([11.1,22.2,43.3,84.4,165.5,326.6,647.7]):

Imodfun:=proc(V)
  global METH;
  if not type(V,numeric) then
    'procname'(V);
  else
    CurveFitting:-ArrayInterpolation(Vdata,Idata,V,method=METH);
  end if;
end proc:

METH:=cubic:
Pcub:=plot(Imodfun(V), V=Vdata[1]..Vdata[-1], color=green):
METH:=spline:
Pspl:=plot(Imodfun(V), V=Vdata[1]..Vdata[-1], color=red):

plots:-display(Pcub,Pspl,
  plots:-pointplot(,symbol=circle,symbolsize=10) );

METH:=cubic:
Pdcub:=plot(D(Imodfun), Vdata[1]..Vdata[-1],color=green):
METH:=spline:
Pdspl:=plot(D(Imodfun), Vdata[1]..Vdata[-1],color=red):

plots:-display(Pdcub,Pdspl);
First 12 13 14 15 16 17 18 Last Page 14 of 48