C_R

3412 Reputation

21 Badges

5 years, 312 days

MaplePrimes Activity


These are answers submitted by C_R

Prefixnotation is especially useful in functional programming where a sequence of commands is grouped using the composition operator @. It is way of separating functions (i.e. Maple commands) from arguments. Otherwise, nested commands with allot of parentheses have to be used.

There is little documentation on that and it is IMO not self-explaining. I got a pretty good explanation here in MaplePrimes. And here is a nice example of a self-made “unary prefix operator” that simplifies a rational expression. (I have put operator in quotation marks since depending on what is does it can equally be called function or command in Maple. In this example it is rather a command.)

I agree that for only two operands infix notation is often preferable. But in functional programming the number of operands (or better arguments if interpreted from a procedure perspective) can change in the course of processing composed commands.

I think a bit more explanatory documentation on top of dharr's examples would not harm.

I hope this completes the picture why Maple sometimes provides infix and prefix notation for the same operation. Functional programming without some prefix variants would be less powerful.

Units_and_Integrals_reply.mw

I have fixed a few errors

I assume that you are familar with Maple styles and in your question you want, for example, to change only the font of all Maple input in an existing document.

What you have to do is:

  • hide all other screen content using view-> show/hide contents
  • select all remaining content with the mouse
  • change the font
  • show the hinden content using view-> show/hide contents

Units in (4) are correct if t in (4) is evaluated with units. With units standard int must be mapped over a column vector (I don't know why that is).

restart

with(VectorCalculus); SetCoordinates('cartesian[x, y, z]')

with(ScientificConstants)

"with(Units[Simple]):"

NULL

g := evalf(Constant(g, units))

9.80665*Units:-Unit(m/s^2)

(1)

`#mover(mi("\`v__0\`"),mo("&rarr;"))` := `<,>`(v[0]*Unit('m'/'s'), 0*Unit('m'/'s'), 0*Unit('m'/'s'))

Vector[column](%id = 36893490633426056244)

(2)

`#mover(mi("a"),mo("&rarr;"))` := `<,>`(0*Unit('m'/'s'^2), g, 0*Unit('m'/'s'^2))

Vector[column](%id = 36893490633436725244)

(3)

int(`#mover(mi("a"),mo("&rarr;"))`, t)

Vector[column](%id = 36893490633436713196)

(4)

"eval(?,t=2&lobrk;s&robrk;)"

Vector[column](%id = 36893490633436699932)

(5)

restart;
with(Units:-Standard):

`#mover(mi("a"),mo("&rarr;"))` := <0*Unit(('m')/'s'^2), g, 0*Unit(('m')/'s'^2)>

Vector(3, {(1) = 0, (2) = g, (3) = 0})

(6)

map(int,(6),t)

Vector(3, {(1) = 0, (2) = g*t, (3) = 0})

(7)

int~((6),t)

Vector(3, {(1) = 0, (2) = g*t, (3) = 0})

(8)

Download Units_and_Integrals_reply.mw

I assume that you are able to generate plots for the exact and the approximate solution using the plot command.

Each of these plots have to be assigned to a name with the assignment operator :=

Then you can plot boths solutions in one plot using plots:- display command.

If you have these solutions in a worksheet and you can't manage to compare them, you could upload them using the green arrow.

To your question: How to calculate this integral?

My preliminary answer is: Don't use Maples sophisticated algorithms.

Zoomed in it looks like this

This could explain why Maples sophisticated integration methods have difficulties. They will try to follow all the oscillations.

Annother reason could be: The function you want to integrate is not always real valued also (see the attachment).

Could you tell a bit more about your integrand?

PS:: I have just seen Roubens answer. Zoomed it looks also jagged.

Integral_reply.mw

Do you want the new constant to have the highest index? If filling of gaps is acceptable

n:=1; 
do n++ until not(has(_C||n,myconstants)):
new_constant := _C ||(n)
 

Otherwise, in your example myconstants is a set and the element with the highest index is to the right.

n:=1; 
do n++ until is(_C||n=myconstants[-1]):
new_constant := _C ||(n)

I guess there are shorter ways

but you will get three bulky solutions

For clarity I assume in the following that you combine the coefficients into new ones.
Your equations will look like this

eq1:=cos(3*phi)-a*cos(phi)=b

cos(3*phi)-a*cos(phi) = b

(1)

eq2:=sin(3*phi)-a*sin(phi)=c

sin(3*phi)-a*sin(phi) = c

(2)

 

Solve for phi

expand(eq1,trig);
sols:=solve(%,[phi])

4*cos(phi)^3-3*cos(phi)-a*cos(phi) = b

 

[[phi = arccos((1/6)*(27*b+3*(-3*a^3-27*a^2+81*b^2-81*a-81)^(1/2))^(1/3)-6*(-(1/12)*a-1/4)/(27*b+3*(-3*a^3-27*a^2+81*b^2-81*a-81)^(1/2))^(1/3))], [phi = Pi-arccos((1/12)*(27*b+3*(-3*a^3-27*a^2+81*b^2-81*a-81)^(1/2))^(1/3)-3*(-(1/12)*a-1/4)/(27*b+3*(-3*a^3-27*a^2+81*b^2-81*a-81)^(1/2))^(1/3)-((1/2)*I)*3^(1/2)*((1/6)*(27*b+3*(-3*a^3-27*a^2+81*b^2-81*a-81)^(1/2))^(1/3)+6*(-(1/12)*a-1/4)/(27*b+3*(-3*a^3-27*a^2+81*b^2-81*a-81)^(1/2))^(1/3)))], [phi = arccos(-(1/12)*(27*b+3*(-3*a^3-27*a^2+81*b^2-81*a-81)^(1/2))^(1/3)+3*(-(1/12)*a-1/4)/(27*b+3*(-3*a^3-27*a^2+81*b^2-81*a-81)^(1/2))^(1/3)-((1/2)*I)*3^(1/2)*((1/6)*(27*b+3*(-3*a^3-27*a^2+81*b^2-81*a-81)^(1/2))^(1/3)+6*(-(1/12)*a-1/4)/(27*b+3*(-3*a^3-27*a^2+81*b^2-81*a-81)^(1/2))^(1/3)))]]

(3)

You get 3 solutions for for phi that you can substitute one after each other into eq2

NULL

simplify(subs(sols[1],eq2))

(1/18)*(-9*(((2/3)*a-2)*(27*b+3*(-3*a^3-27*a^2+81*b^2-81*a-81)^(1/2))^(2/3)+(3*b+(1/3)*(-3*a^3-27*a^2+81*b^2-81*a-81)^(1/2))*(27*b+3*(-3*a^3-27*a^2+81*b^2-81*a-81)^(1/2))^(1/3)+(a+3)^2)/(27*b+3*(-3*a^3-27*a^2+81*b^2-81*a-81)^(1/2))^(2/3))^(1/2)*((-a+3)*(27*b+3*(-3*a^3-27*a^2+81*b^2-81*a-81)^(1/2))^(2/3)+(9*b+(-3*a^3-27*a^2+81*b^2-81*a-81)^(1/2))*(27*b+3*(-3*a^3-27*a^2+81*b^2-81*a-81)^(1/2))^(1/3)+3*(a+3)^2)/(27*b+3*(-3*a^3-27*a^2+81*b^2-81*a-81)^(1/2))^(2/3) = c

(4)

 

Depending on the initial parameters you might be able to simplify further and select an equation that describes your problem

Download elliminate_phi.mw

Maple is capable to integrate up to the first 200 roots without any special methods.
 

evalf(Int(abs(sin(x^4))/(sqrt(x) + x^2), x = 0 .. 5))

 

50000 roots are reached at about x=20.

 

evalf(eval((Pi*i)^(1/4), i = 50000))

However, maxintervals=50000 is in any case not sufficient to place between all roots up to infinity at least one interval. But this is somehow  irrelevant, because Maple does not even get to the point where 50000 intervals are reached.

What happens in your case: Maple subdivides the integrations range applying numerical analysis to the integrand and stops at 200 which matches the 200 roots (see above).

Error, (in evalf/int) NE_QUAD_MAX_SUBDIV:
  The maximum number of subdivisions has
  been reached: max_num_subint = 200

Reducing the numerical accuracy to 2 significant digits does prevent Maple from detecting more than 200 roots. The evaluated integrand becomes too small and is rounded to zero. This limits the detected roots.

Increasing to 3 allows you to integrate up to 17

evalf[3](Int(abs(sin(x^4))/(sqrt(x) + x^2), x = 0 .. 16.999))

Increasing to 4, and you can go up to 13

evalf[4](Int(abs(sin(x^4))/(sqrt(x) + x^2), x = 0 .. 13.0001))

I think this challenging integral needs to be attacked from a different angle. Maybe something with limits or an adaptive method that integrates only up to the next 200 roots until a convergence criteria is reached.

Update:
The error messages the OP was observing were contradictory in the sense that in one instance maxintervals was ignored, in the other not. It is still unclear why forcing maxintervals did not work in one case. Anyway, to substantiate my answer, I have added a worksheet which highlights essential parts of Maples integration control. It can be seen how evalf[n] sets integration error tolerances, at which instance the integrand is set to zero, and how Maple iteratively subdivides an integration interval, tries to apply default methods to a new subinterval, switches to specific methods and changes the methods when numerically required.

adaptive_integration.mw

If your region/domain is any other than the coords option of plot offers, then you have to write your own plot routine. What I mean by this is that in (for example) cartesian coordinates only rectangular regions can be mapped.

In the example below I have punched a hole in the domain of the quadrant Re(z)>0, Im(z)>0.

restart;
f:=z^2;# try f:=z to see how the domain is restricted

z_map:=proc(f,re,im) 
  if(sqrt((re-5)^2+(im-3)^2)>=3)then
    eval(f,z=re+I*im);
  else
    NULL;
  end if;
end proc;

p_re:=plots:-display(seq(plot([Re('z_map(f,re,im)'),Im('z_map(f,re,im)'),im=0..10]),re=0..20)):
p_im:=plots:-display(seq(plot([Re('z_map(f,re,im)'),Im('z_map(f,re,im)'),re=0..20]),im=0..10),color=green):
plots:-display(p_re,p_im,scaling=constrained)

The above could be expanded to a routine where you could pass the restrictions as arguments.

Update:

A triangular domain can be found here.

 

X := Matrix([[1, -X3_3/2 - 1/2, 0, -X2_3], [-X3_3/2 - 1/2, -2*X3_4 - 1, X2_3, 0], [0, X2_3, X3_3, X3_4], [-X2_3, 0, X3_4, 1]]);
vars := [X3_4, X3_3, X2_3];
w := A^3 - A;
rootz := RootOf(w, A);
Pols := [(-A^2 + 1)/(3*A^2 - 1), (-A^2 - 1)/(3*A^2 - 1), A*(3*A^2 - 1)*1/(3*A^2 - 1)];
vals := {allvalues(eval(Pols, A = rootz))};
ecs := map(x -> convert(vars = x, listofequations)[], vals);
Xs := [seq(eval(X, ecs[x]), x = 1 .. numelems(vals))];
Xz := [seq(subs(ecs[x], X), x = 1 .. numelems(vals))];
with(LinearAlgebra);
Eigs := map(x -> Eigenvalues(x), Xs);
BolEigs := map~(type, Eigs, poszero);
evalf(Eigs);

For BolEigs there is probably a shorter version without seq.

Edit: Must be poszero not positive. I have corrected that in the above

Edit2: BolEigs replaced by something shorter in the above

 

@Kitonum 

Using allsolutions, the root selector of 0.25 indicates that Maple tackles part of the problem numerically

restart;
solve(x^x = 1/sqrt(2), allsolutions);
allvalues(%);
evalf~([%]);

However, Maple 2023 could derive the exact solution

restart;
solve(x^x = 1/sqrt(a), {x});
eval(%, a = 2);
solve(x^x = 1/sqrt(a), {x}, allsolutions);
eval(%, [a = 2, _Z1 = 0, _Z3 = 0]);
eval(`%%`, [a = 2, _Z1 = 0, _Z3 = -1]);

I can't check right now if this is possible in Maple 2018.

I was wondering how many rational valued pairs (a, x) exist and if its worth writing code for this.

To verify what the worksheet does, I have plotted the path a cube on an air baring (almost friction less) would take when it is released from the carrousel.

Position of "cube" in the inertia frame after release from the carrousel at t=0.

r = r__0+I*v__0*t

r = r__0+I*v__0*t

(1)

lhs(r = r__0+I*v__0*t) = convert(rhs(r = r__0+I*v__0*t), polar)

r = polar(abs(r__0+I*v__0*t), argument(r__0+I*v__0*t))

(2)

convert(r = polar(abs(r__0+I*v__0*t), argument(r__0+I*v__0*t)), exp)

r = abs(r__0+I*v__0*t)*exp(I*argument(r__0+I*v__0*t))

(3)

with the speed of the "cube" at t=0

v__0 = r__0*`&omega;__0`

v__0 = r__0*omega__0

(4)

subs(v__0 = r__0*omega__0, r = abs(r__0+I*v__0*t)*exp(I*argument(r__0+I*v__0*t)))

r = abs(r__0+I*r__0*omega__0*t)*exp(I*argument(r__0+I*r__0*omega__0*t))

(5)

The angular orientation of carousel with respect to the inertia frame

`&varphi;` = `&omega;__0`*t

varphi = omega__0*t

(6)

Describing r in the rotating carousel frame by rotating the position back by -`&varphi;` 

r__car = rhs(r = abs(r__0+I*r__0*omega__0*t)*exp(I*argument(r__0+I*r__0*omega__0*t)))*exp(-I*`&varphi;`)

r__car = abs(r__0+I*r__0*omega__0*t)*exp(I*argument(r__0+I*r__0*omega__0*t))*exp(-I*varphi)

(7)

`assuming`([simplify(subs(varphi = omega__0*t, r__car = abs(r__0+I*r__0*omega__0*t)*exp(I*argument(r__0+I*r__0*omega__0*t))*exp(-I*varphi)))], [`&omega;__0` > 0, r__0 > 0, v__0 > 0, t > 0])

r__car = r__0*(1+I*omega__0*t)*exp(-I*omega__0*t)

(8)

This is a logarithmic spiral superposed with the unit circle

subs(r__0 = 1, `&omega;__0` = evalf((1/12)*Pi), r__car = r__0*(1+I*omega__0*t)*exp(-I*omega__0*t))

r__car = (1+(.2617993878*I)*t)*exp(-(.2617993878*I)*t)

(9)

plots:-complexplot([10*exp(I*t), rhs(r__car = (1+(.2617993878*I)*t)*exp(-(.2617993878*I)*t))], t = 0 .. 37.76, view = [-10 .. 10, -10 .. 10], gridlines, axes = boxed, title = (r__car = (1+(.2617993878*I)*t)*exp(-(.2617993878*I)*t)))

 

NULL

Download carousel_kinematics.mw

This seems to match your simulation. Depending on the reference frame in which the path is defined it is either a spiral or a straight line.

To your question of placing thrusters: The ode “only” allows to add complex thrusters. For example:   

results in a straight line. The path highly depends on the amplitude, the signs of the exponent, the phase and the frequency. The above thrust terms is intersting to play with and shows the advantages of having a complex description.

Now, if you want to trace a specific curve, you are dealing with an inverse dynamics problem: For a given path/trace you have to find a thruster function that adjust the thrust in a way that the integration of the ode matches your curve.

This should be possible if complex Fourier series are used.

  • In a first step a Fourier description of your target curve in the carrousel frame has to be found
  • Secondly, this description has to be transformed to its counterpart in the inertia frame, which (to keep it “simple”) should be done at constant angular velocity of carousel.
  • In the inertia frame the counterpart of the curve is now given in a parametric form, depending on the parameter time. To follow this kinematic path, a complex Fourier force series has to be derived. This is the inverse dynamics part. If I am not mistaken, two times differentiation w.r.t. to the parameter time should yield the force function. This is easy because in an inertia frame we do not have to deal with corriolis terms.
  • Once the force function has been found, this function has to be converted back to the rotating reference frame of the carrousel (this is the easiest part: multiplication by exp(1, -omega*t)).
  • Adding the force function to the ode and integration it should trace the desired curve.

Maple is particularly well suited for such an approach since series order and numerical fidelity can be set parametrically, and there are optimization tools. Would be interesting to see how good this works.

Concering rolling without slipage: A rolling ball with the same mass as the cube will trace a different path since gyroscopic effects lead to external forces acting on the ball that change the balls impulse in the inertia frame. I will try to demonstrate this in a separate post because this is cannot to be solved with a description in the complex plance.

What I do for a horizontal split is copy&pasting the lower part of the table and then deleting the lower part of the table.

So far this also works with equation labels. I have not tried to cut and paste yet, because I like to compare before deleting.

This however does not work when in the body below the table equation labels of the table are used. These references are lost. This use case would justify a split option in the menues.

table_split.mw

I had no time to study the whole thread in detail.

Just in case that you are still looking for a way that uses the same input for real valued input or integer input as well as for definite and semi-definite integrals, you can use the method _Dexp this way

Digits := 32

32

(1)

"evalf(Pi *Int(2*299792458^2*662607015*10^(-8)*10^(-34)/((exp(299792458*662607015*10^(-8)*10^(-34)/(1380649*10^(-6)*10^(-23)*lambda*5772.0)) - 1)*lambda^5), lambda = 0 .. infinity ,method=_Dexp));"

62938592.470335950467548474587301

(2)

evalf(I*Pi*nt((2*299792458^2*662607015)*10^(-8)*10^(-34)/((exp((299792458*662607015)*10^(-8)*10^(-34)/(1380649*10^(-6)*10^(-23)*lambda*5772))-1)*lambda^5), lambda = 0 .. infinity, method = _Dexp))

62938592.470335950467548474587304

(3)

"evalf(Pi*Int(2*299792458^2*662607015*10^(-8)*10^(-34)/((exp(299792458*662607015*10^(-8)*10^(-34)/(1380649*10^(-6)*10^(-23)*lambda*5772.0)) - 1)*lambda^5), lambda = 380*10^(-9) .. 750*10^(-9),method=_Dexp));                                                   "

27551199.571700602410253741952570

(4)

evalf(Pi*(Int((2*299792458^2*662607015)*10^(-8)*10^(-34)/((exp((299792458*662607015)*10^(-8)*10^(-34)/(1380649*10^(-6)*10^(-23)*lambda*5772))-1)*lambda^5), lambda = 380*10^(-9) .. 750*10^(-9), method = _Dexp)))

27551199.571700602410253741952570

(5)

NULL


It is important to use the inert form of int which prevents symbolic evaluation and/or exact integration methods, which in your case leads Maple to treat the integral as a beeing complex.

Download Complex_radiant_excitance_-_integer_vs_real.mw

First 6 7 8 9 10 11 12 Last Page 8 of 17