C_R

1960 Reputation

19 Badges

5 years, 279 days

MaplePrimes Activity


These are answers submitted by C_R

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*`ω__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

`ϕ` = `ω__0`*t

varphi = omega__0*t

(6)

Describing r in the rotating carousel frame by rotating the position back by -`ϕ` 

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*`ϕ`)

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)))], [`ω__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, `ω__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

To your question why the solutions are different:

The solver analyses and classifies the problem and then applies various methods on the input. In your case, the input is different and the solver has (among other things) to remove a trivial solution when deriving p1.

However, the results are correct. It's more a question why is retruns false. I think it is the missing pole at c=0.

With Maple 2023.2 and Physics:-Version() = 1566, is returns true. Which versions do you use?

restart

interface(version)

`Standard Worksheet Interface, Maple 2023.2, Windows 10, November 24 2023 Build ID 1762575`

(1)

p__1 := solve(c = p^a*p^b, p); p__2 := solve(c = p^(a+b), p)

exp(-ln(1/c)/(a+b))

 

exp(ln(c)/(a+b))

(2)

NULL

is(p__1 = p__2)

true

(3)

Comparing your solution

 

-ln(1/c)

-ln(1/c)

(4)

ln(c)

ln(c)

(5)

plots:-complexplot3d(proc (c) options operator, arrow; -ln(1/c) end proc, -2-2*I .. 2+2*I, title = -ln(1/c))

 

plots:-complexplot3d(proc (c) options operator, arrow; ln(c) end proc, -2-2*I .. 2+2*I, title = ln(c))

 
 

 

NULL

plots:-complexplot3d(proc (c) options operator, arrow; -ln(1/c)-ln(c) end proc, -2-2*I .. 2+2*I, title = -ln(1/c)-ln(c))

 

NULL

Download solve_exp.mw

Could this work for you?

restart:

with(Statistics):

U = RandomVariable(Uniform(a, b)),V = RandomVariable(Uniform(c, d));

U = _R, V = _R0

(1)

NULL

Y = U*V^2

Y = U*V^2

(2)

NULL

map(Mean,(subs((1),(2))))

Y = (1/6)*(a+b)*(c^2+c*d+d^2)

(3)

NULL

Unit side note (if this is also of interest for the students):

Using equations to define physical quantities is my prefered way of using units. In the same way, equations are used in the above to declare random variables to derive a general solution for the mean. In a second step a,b,c,d could be replaced by physical quantities (again using equation expressions). Doing the replacement/substitution in a later stage does not expose library commands used in the derivation process to units which increases the success rate.  In the case of Mean this  precaution is not required since Mean seems to handle units and could be used in RandomVariable or with RandomVariable. Both uses would results in a general solution for the mean where a,b,c,d are dimensionless with the danger that wrong values are plugged in for numerical evaluation. For this reason I did not include these two variants (which I consider not a good practise).

In a nutsshell, random variables should be used before using units in a derivation process for statistical magnitudes.

Download Alg_with_RandomVariable.mw

For the x-axis, there is an ugly but simple solution: add spaces.

This works for the y-axis when the label can be rotated (acceptable for words but not for symbols).

For the y-axis in your example, textplot works better.

yLabel:=plots:-textplot([-1,3,"y"], 'align'={'below', 'right'},font=["times","roman",10]):
p:=plot((x -> 1/abs(x))(x),x=-4..4):
plots:-display(p,yLabel, view = [-4 .. 4, 0 .. 3],labels=["                                                                        x",""])

You can use type to check if an expression is of tye equation and use it in an if statement in combination with a not command

https://www.maplesoft.com/support/help/Maple/view.aspx?path=type/equation&cid=507

Same for other types

Maple assumes that k and x in your expression are complex.

To get the expected ouput, assumptions on k and x are required

assume(x::real, k::real);
abs(exp(k*x*I));

or after restart

abs(exp(k*x*I));
(simplify(%) assuming (k::real, x::real));

or

abs(exp(k*x*I));
evalc(%);

which assumes that all variables are real valued

Edit: If you want less ticks on the axes, add the option tickmarks

plot(1/x, x = 1 .. 10, tickmarks = [5, 5])

More convenient is to select the plot with a right-click and go to axes, properties, tickmarks

 

If you want to plot points only:

  • With the plot command you can add the option numpoints=5 and change the the plotstyle to point (either by an option or the plot menu when the plot is selected)
plot(1/x, x = 1 .. 10, numpoints = 5, style = point)
  • Or: To plot points you can use the command plots:-pointplot in combiantion with a list of points you have to generate

The number of activations depends on the license bought.

For your tests, you can either contact support for a temporary additional license or install a test version from a website. Via support you might get an extended test period if Maplesoft is interested in your tests. This could be if the tests are extended and systematic (which seems to be the case).

Another way is to use an older version of Maple for the test but this might not be conclusive. 

Or subscribe for beta testing to get an addional license for tests.

Regarding your question why Maple is not providing an explicit solution:


Maple classifies this ode as separable

and applies

DEtools:-separablesol(ode)

with the solution

This solution does not satisfy the initial condition and therefore no solution is returned. The IC represents a singularity for this solution (not for your solution and not for vv's last solution).

Only if called with the implict option, c__1 can be determined

solve(dsolve(ode,implicit),{c__1})

The separation method leads to a dead end without a proposal where to go from there (e.g. series).
I had no time to check if DEtools commands could lead to the solutions from nn and vv. 


 

1 2 3 4 5 6 7 Last Page 2 of 10