Advanced Search

_Maxim_

729 Reputation

12 Badges

8 years, 327 days
Contact _Maxim_

MaplePrimes Activity


These are Posts that have been published by _Maxim_

Hyperasymptotics...

_Maxim_ 729
hypergeom meijerg
February 14 2018
7 3

One case where an "expansion beyond all orders" may be needed is investigating the asymptotic behavior of the difference of two functions with coinciding dominant series.

We are interested in the asymptotic behavior of F(z) for large positive z:

> 

h1 := proc (z) options operator, arrow; hypergeom([1/2], [5/4, 3/2, 7/4], z) end proc; h2 := proc (z) options operator, arrow; (3/4)*sqrt(2*Pi)*hypergeom([1/4], [3/4, 5/4, 3/2], z)/(GAMMA(3/4)*z^(1/4)) end proc; F := proc (z) options operator, arrow; h1(z)-h2(z)+(3/8)*sqrt(Pi)/sqrt(z) end proc

series does not succeed:

> 

series(F(z), z = infinity, 20)

O((1/z)^(23/3))*exp(3/(1/z)^(1/3))

 (1)

The reason is that the dominant terms containing the factor exp(3*z^(1/3)) in the two hypergeometric functions cancel exactly, and we have to look for the subdominant terms.

The order of the leading terms can be found from DETools:-formal_sol:

> 

deq1 := FunctionAdvisor(DE, h1(z), f(z))[2, 1]

diff(diff(diff(diff(f(z), z), z), z), z) = -(15/2)*(diff(diff(diff(f(z), z), z), z))/z-(195/16)*(diff(diff(f(z), z), z))/z^2+(1/32)*(32*z-105)*(diff(f(z), z))/z^3+(1/2)*f(z)/z^3

 (2)
> 

DETools:-formal_sol(deq1, f(z), z = infinity, order = 0)

[(1/z)^(1/2), -exp(-3/(-1/z)^(1/3))/z, -exp(3*(-1)^(1/3)/(-1/z)^(1/3))/z, -exp(-3*(-1)^(2/3)/(-1/z)^(1/3))/z]

 (3)

As expected, one of the solutions (the third one for positive z) contains the exp(3*z^(1/3)) factor, the leading term being of the order exp(3*z^(1/3))/z.

Another, subdominant, solution is algebraic and, in fact, is a series containing only one term, as 1/z^(1/2) is an exact solution. It will turn out that the algebraic part in F(z) also cancels out.

Thus we have to look for the subsubdominant terms, which contain decaying exponentials. We will accomplish this by applying the steepest descent method to the integral representations of h1(z) and h2(z).

> 

ms := convert([h1(z), h2(z)], MeijerG)

[(3/32)*Pi*2^(1/2)*MeijerG([[1/2], []], [[0], [-1/4, -1/2, -3/4]], -z), (3/32)*2^(1/2)*Pi*MeijerG([[3/4], []], [[0], [1/4, -1/4, -1/2]], -z)/z^(1/4)]

 (4)
> 

m2g := proc (m, y) local a, b, c, d; a, b := op(op(1, m)); c, d := op(op(2, m)); -((1/2)*I)*mul(`~`[GAMMA](`~`[`-`](1+y, a)))*mul(`~`[GAMMA](`~`[`-`](c, y)))*op(3, m)^y/(Pi*mul(`~`[GAMMA](`~`[`-`](b, y)))*mul(`~`[GAMMA](`~`[`-`](1+y, d)))) end proc
> 

gs := applyrule(conditional(e::anything, _op(0, e) = MeijerG) = 'm2g(e, y)', ms)

[-((3/64)*I)*2^(1/2)*GAMMA(1/2+y)*GAMMA(-y)*(-z)^y/(GAMMA(5/4+y)*GAMMA(3/2+y)*GAMMA(7/4+y)), -((3/64)*I)*2^(1/2)*GAMMA(1/4+y)*GAMMA(-y)*(-z)^y/(z^(1/4)*GAMMA(3/4+y)*GAMMA(5/4+y)*GAMMA(3/2+y))]

 (5)
> 

gs[2] := combine(eval(gs[2], [1/z^(1/4) = exp(I*Pi*(1/4))/(-z)^(1/4), y = y+1/4]), power)

-((3/64)*I)*GAMMA(1/2+y)*((1/2)*2^(1/2)+((1/2)*I)*2^(1/2))*(-z)^y*GAMMA(-1/4-y)*2^(1/2)/(GAMMA(1+y)*GAMMA(3/2+y)*GAMMA(7/4+y))

 (6)

h1(z) and h2(z)are the integrals of gs[1] and of gs[2] over the same path, which is a loop encircling the poles ofGAMMA(-y) and of GAMMA(-1/4-y). Now a standard technique is to extend the integration contour far to the left, while still keeping both endpoints at "+infinity". Then the arguments of the gamma functions can be made large everywhere on the integration path, and the gamma functions can be replaced by their asymptotic approximations.

When moving the contour, we have to take into account the pole of the integrand at y = -1/2. The other poles of GAMMA(1/2+y) will be cancelled by the zeros of 1/GAMMA(3/2+y), which is why the algebraic part of the expansion will contain the single term of the order 1/z^(1/2).

This is the negative of the third term in F(z):

> 

`assuming`([simplify((2*Pi*I)*residue(gs[1]-gs[2], y = -1/2))], [z > 0])

-(3/8)*Pi^(1/2)/z^(1/2)

 (7)

Expanding the gamma functions produces terms containing exp(-I*Pi*y) and exp(I*Pi*y)

> 

`assuming`([simplify(convert(MultiSeries:-series((gs[1]-gs[2])/z^y, y = -infinity, 1), polynom))], [z > 0]); collect(convert(%, exp), exp)

(-3/64-(3/64)*I)*(-1/y)^(1/2)*exp(-3*(ln(-y)-1)*y)*exp(-I*Pi*y)/(y^3*Pi^(1/2))+(3/64-(3/64)*I)*(-1/y)^(1/2)*exp(-3*(ln(-y)-1)*y)*exp(I*Pi*y)/(y^3*Pi^(1/2))

 (8)

As we shall see, those terms have saddle points y0(z) = exp(`&+-`((1/3)*(2*Pi*I)))*z^(1/3) located in the left half-plane and contribute exponentially small factors exp(3*y0(z)). The terms for which the saddle point would be located at y = z^(1/3) have cancelled out, thus cancelling the exponentially large contributions. Another possible way to achieve the same result was to write h1(z)-h2(z) as a single Meijer G-function -(3/32)*MeijerG([[1/2], []], [[-1/4, 0], [-3/4, -1/2]], z).

We write the first term above in the form g(y)*exp(f(y)):

> 

f := proc (z, y) options operator, arrow; -3*y*(ln(-y)-1)-I*Pi*y+y*ln(z) end proc
> 

g := proc (y) options operator, arrow; (-3/64-(3/64)*I)*sqrt(-1/y)/(sqrt(Pi)*y^3) end proc
> 

diff(f(z, y), y)

-3*ln(-y)-I*Pi+ln(z)

 (9)

For this to become zero, we need argument(-y) = -(1/3)*Pi, and thus y = exp((1/3)*(2*I)*Pi)*z^(1/3). We can visualize the paths where the imaginary part of f(z, y) stays constant. The path of the steepest descent is the one that goes through the saddle point in the direction exp(I*Pi*(1/3)); the blue color indicates smaller values of the real part of f(z, y):

> 

y0 := proc (z) options operator, arrow; exp(((2/3)*I)*Pi)*z^(1/3) end proc
> 

(proc () local z; z := 2; plots:-display(plots:-contourplot(Re(f(z, u+I*v)), u = -5 .. 5, v = -5 .. 5, contours = ([seq])(Re(f(z, y0(z)))+i, i = -30 .. 6, 6), filledregions, coloring = [blue, red], grid = [100, 100]), plots:-implicitplot(Im(f(z, u+I*v)-f(z, y0(z))), u = -5 .. 5, v = -5 .. 5, gridrefine = 5, color = green), plot([cos((1/3)*Pi)*xi+Re(y0(z)), sin((1/3)*Pi)*xi+Im(y0(z)), xi = -3 .. 3], linestyle = dot, color = white), axes = boxed) end proc)()

 

The real part of f(z, y) has a maximum along this path at y0(z).

> 

`assuming`([(`@`(`@`(simplify, evalc), series))(f(z, y0(z)+exp(I*Pi*(1/3))*xi), xi = 0, 3)], [z > 0]); quad := convert(%, polynom)

series((3/2)*(-1+I*3^(1/2))*z^(1/3)-((3/2)/z^(1/3))*xi^2+O(xi^3),xi,3)

 (10)

Now we can compute the lead asymptotic term contributed by the saddle point y0(z):

> 

lt1 := `assuming`([(`@`(simplify, evalc))(g(y0(z))*exp(I*Pi*(1/3))*(int(exp(quad), xi = -infinity .. infinity)))], [z > 0])

-(1/64)*exp(-(3/2)*z^(1/3))*3^(1/2)*((-1+I)*cos((3/2)*z^(1/3)*3^(1/2))+(-1-I)*sin((3/2)*z^(1/3)*3^(1/2)))*2^(1/2)/z

 (11)

We repeat the same procedure for the second term of the integrand.

> 

f := proc (z, y) options operator, arrow; -3*y*(ln(-y)-1)+I*Pi*y+y*ln(z) end proc
> 

g := proc (y) options operator, arrow; (3/64-(3/64)*I)*sqrt(-1/y)/(sqrt(Pi)*y^3) end proc
> 

diff(f(z, y), y)

-3*ln(-y)+I*Pi+ln(z)

 (12)
> 

y0 := proc (z) options operator, arrow; exp(-((2/3)*I)*Pi)*z^(1/3) end proc

The direction should be chosen as exp((1/3)*(2*I)*Pi) to be consistent with the direction of the integration contour, which goes from the lower to the upper half-plane.

> 

lterm := proc (gy, fy, eq, dir) options operator, arrow; (eval(gy*exp(fy), eq))*dir*sqrt(-2*Pi/((eval(diff(fy, `$`(y, 2)), eq))*dir^2)) end proc
> 

lt2 := `assuming`([(`@`(simplify, evalc))(lterm(g(y), f(z, y), y = y0(z), exp((1/3)*(2*I)*Pi)))], [z > 0])

(1/64)*exp(-(3/2)*z^(1/3))*((1+I)*cos((3/2)*z^(1/3)*3^(1/2))+(1-I)*sin((3/2)*z^(1/3)*3^(1/2)))*3^(1/2)*2^(1/2)/z

 (13)

Combining the two results yields the leading term of F(z). The next terms can be obtained by expanding gs[1] and gs[2] to higher orders.

> 

Fasympt := unapply(simplify(lt1+lt2), z)

proc (z) options operator, arrow; (1/32)*exp(-(3/2)*z^(1/3))*3^(1/2)*2^(1/2)*(cos((3/2)*z^(1/3)*3^(1/2))+sin((3/2)*z^(1/3)*3^(1/2)))/z end proc

 (14)
> 

(proc () Digits := 50; plot(`~`[`*`](exp((3/2)*z^(1/3)), [F(z), Fasympt(z)]), z = 1000 .. 10000, linestyle = [solid, dot], thickness = [1, 5], axes = frame) end proc)()

 

Download steep.mw

An involute of y=x^3...

_Maxim_ 729
physics differential-geometry
November 12 2017
3 5

The problem is to analyze the behavior of an involute of the cubical parabola near the singular points and at infinity.

To do this efficiently, it is best to work with partially evaluated expressions.

We want to investigate the singular points of an involute of the cubical parabola y = x^3.

If the curve is given parametrically by "conjugate(r)(t)", an involute is given by

"conjugate(E)(t)=conjugate(r)(t)-(conjugate(r)'(t))/(|conjugate(r)'(t)|)*(∫)[a]^(t)|conjugate(r)'(tau)| ⅆtau"

What we want to do is to leave the integral unevaluated and compute only its derivatives. The simplest way is to define S := proc (a, t) options operator, arrow; Int(s(tau), tau = a .. t) end proc. Both diff and evalf will be handled automatically by Maple. We'll do it in a slightly more complicated way, leaving S(a, t) completely inert with the exception of the condition S(a, a)=0 and implementing a more efficient numerical evaluator.

> 

forget(diff); forget(series); r := proc (t) options operator, arrow; [t, t^3] end proc; s := unapply(sqrt((diff(r(t)[1], t))^2+(diff(r(t)[2], t))^2), t); S := proc (a, t) options operator, arrow; `if`(t = a, 0, 'S(a, t)') end proc; `diff/S` := proc (at, ft, t) options operator, arrow; s(ft)*(diff(ft, t))-s(at)*(diff(at, t)) end proc; E := unapply(r(t)-`~`[`*`](diff(r(t), t), S(a, t)/s(t)), [a, t])
> 

forget(evalf); `evalf/S` := (proc () local acur, Ssol; acur := undefined; Ssol := subs(dsolve({diff(f(t), t) = s(t), f(a) = 0}, f(t), numeric, parameters = [a], output = listprocedure), f(t)); proc (a, t) if [a, t]::(list(numeric)) then if a <> acur then Ssol(parameters = [a]); acur := a end if; Ssol(t) else 'S(a, t)' end if end proc end proc)()
> 

plot([[op(r(t)), t = -1 .. 3], [op(E(1, t)), t = -1.5 .. 4]], view = -1 .. 3, scaling = constrained)

 
> 

simplify(diff(E(a, t), t))

[18*S(a, t)*t^3/(9*t^4+1)^(3/2), -6*t*S(a, t)/(9*t^4+1)^(3/2)]

 (1)

The singular points are easily determined now: we need either S(a, t)=0, which implies t=a, or t=0. We assume a>0.

Compute the Taylor series around t=a first.

> 

`~`[series](E(a, t+a), t = 0, 4); ser := collect(convert(%, polynom), t, normal)

[-18*a^2*(6*a^4-1)*t^3/(9*a^4+1)^2+9*a^3*t^2/(9*a^4+1)+a, 2*(36*a^4-1)*t^3/(9*a^4+1)^2-3*a*t^2/(9*a^4+1)+a^3]

 (2)

The center of the expansion E(a, a) is the point [a, a^3] on the cubical parabola.

There are two quadratic terms, which give us the slope of the tangent line. The tangent coincides with the normal to the cubical parabola at this point.

Rotate the axes to make the tangent vertical.

> 

rot := proc (e, a) options operator, arrow; convert(Student:-LinearAlgebra:-RotationMatrix(a).Vector(e), list) end proc
> 

rot(ser-E(a, a), (`@`(arctan, op))(`~`[coeff](ser, t, 2))); ser := `assuming`([collect(%, t, normal)], [a > 0])

[-12*a^2*t^3/(9*a^4+1)^(3/2), -2*(162*a^8+9*a^4-1)*t^3/(9*a^4+1)^(5/2)+3*a*t^2/(9*a^4+1)^(1/2)]

 (3)

In the new coordinates, the leading terms are t^3 and t^2, and the involute has the shape of a semicubical parabola. When t increases, the involute moves from the first to the second quadrant. In terms of x and y, `~`[x]*alpha*y^(3/2), and the coefficient alpha is found from the two leading terms. For t>0:

> 

`assuming`([simplify(coeff(ser[1], t, 3)/coeff(ser[2], t, 2)^(3/2))], [a > 0])

-(4/3)*a^(1/2)*3^(1/2)/(9*a^4+1)^(3/4)

 (4)

Now consider the second singular point t=0 in the original coordinate system.

> 

ser := convert(`~`[series](E(a, t), t = 0, 6), polynom)

[-S(a, 0)+(9/2)*S(a, 0)*t^4+(18/5)*t^5, -3*S(a, 0)*t^2-2*t^3]

 (5)

The center is the point [-S(a, 0), 0], which corresponds to the point [0, 0] on the cubical parabola. The leading terms t^4 and t^2 indicate that the tangent line is vertical and that this is a return point, with both branches lying in the second quadrant in the coordinate system shifted by [-S(a, 0), 0].

In terms of x and y, `~`[x]*alpha*y^2, and to determine the next term, we need to analyze the two series expansions. t^4 gives a t^4 term and fixes alpha. t^2*(t^3) gives a t^5 term, but the coefficient doesn't match the coefficient at t^5 in x. So the next term has to be y^(5/2) in order to match t^5. For t>0, (b*t^3+a*t^2)^(5/2) = t^5*(b*t+a)^(5/2), which gives a^(5/2)*t^5 plus higher order terms. No other terms contribute, and we can equate the coefficients at t^4 and t^5 in x and y:

> 

collect(alpha*ser[2]^2+beta*coeff(ser[2], t, 2)^(5/2)*t^5, t); (`@`(simplify, solve))(`~`[coeff](%-ser[1], t, [4, 5]), {alpha, beta})

4*alpha*t^6+(beta*(-3*S(a, 0))^(5/2)+12*alpha*S(a, 0))*t^5+9*alpha*S(a, 0)^2*t^4

  

{alpha = (1/2)/S(a, 0), beta = -(4/45)*3^(1/2)/(-S(a, 0))^(5/2)}

 (6)

We have obtained the coefficients in `~`[x]*alpha*y^2+beta*y^(5/2). alpha and beta are negative, while for t<0, beta has the opposite sign, and the branch corresponding to t<0 is the one which is closer to the vertical tangent.

This could have been done by using Groebner:-Basis to eliminate t from the two series and then using algcurves:-puiseux. But it is still necessary to determine the truncation order.

There is also a separate case a=0.

> 

ser := `~`[series](E(0, t), t = 0, 6)

[series((18/5)*t^5+O(t^6),t,6), series(-2*t^3+O(t^7),t,7)]

 (7)

Now the singularity is an inflection point, and the involute has a vertical tangent and moves from the second to the fourth quadrant. In the fourth quadrant, `~`[x]*alpha*(-y)^(5/3), and we find alpha in the same way as before:

> 

coeff(ser[1], t, 5)/(-coeff(ser[2], t, 3))^(5/3)

(9/10)*2^(1/3)

 (8)

Finally let's consider the asymptotic behavior of an involute at infinity. Now we do have to investigate S(a, t). One possible way is to apply the binomial formula to s(t) = 3*t^2*sqrt(1+1/(9*t^4))  and integrate term by term:

> 

`assuming`([int(3*tau^2*binomial(1/2, k)*(9*tau^4)^(-k), tau = a .. t)], [0 < a and a < t])

3*binomial(1/2, k)*(9^(-k)*a^(-4*k+3)-9^(-k)*t^(-4*k+3))/(4*k-3)

 (9)

This is the power series expansion of S(a, t) for large t, valid for 9*a^4 > 1. Substituting it into E(a, t), we obtain the required asymptotics.

> 

E(a, t)

[-S(a, t)/(9*t^4+1)^(1/2)+t, -3*t^2*S(a, t)/(9*t^4+1)^(1/2)+t^3]

 (10)

In the x component, we get -(1/3)*t+o(1)+t. In the y component, we get -t^3-C+o(1)+t^3. Thus the involute approaches a horizontal asymptote when t goes to +infinity, and the constant term -C gives the position of the asymptote:

> 

`C__+` := `assuming`([unapply(sum(3*binomial(1/2, k)*9^(-k)*a^(-4*k+3)/(4*k-3), k = 0 .. infinity), a)], [a > 1/sqrt(3)])

proc (a) options operator, arrow; -a^3*hypergeom([-3/4, -1/2], [1/4], -(1/9)/a^4) end proc

 (11)

That in fact extends to all positive a. For the asymptotics at -infinity, first we compute the integral over [a, -a], this time using the power series for small t. For the integral from -a to -t, S(-a, -t) = -S(a, t), and we need to subtract the value found above:

> 

`C__-` := `assuming`([unapply(sum(int(binomial(1/2, k)*(9*tau^4)^k, tau = a .. -a), k = 0 .. infinity)-`C__+`(a), a)], [0 < a and a < 1/sqrt(3)])

proc (a) options operator, arrow; -2*a*hypergeom([-1/2, 1/4], [5/4], -9*a^4)+a^3*hypergeom([-3/4, -1/2], [1/4], -(1/9)/a^4) end proc

 (12)

For negative a, `C__+` and `C__-` are swapped.

This could have been done by converting S(a, t) into a difference of two integrals of hypergeometric type, for which Maple is able to find the asymptotics automatically.

> 

inv := proc (a, T) plots:-display(plot([op(r(t)), t = -1 .. 2], color = black), plot([op(E(a, t)), t = -1.5 .. T], color = red), plottools:-line(r(T), E(a, T), linestyle = dot), map(proc (y) options operator, arrow; plottools:-line([-1, y], [3, y], linestyle = spacedot) end proc, [-`C__+`(a), -`C__-`(a)]), scaling = constrained, view = [-1 .. 3, -1 .. 3]) end proc
> 

Explore(inv(a, t), a = -.3 .. 1.3, t = -1. .. 4., initialvalues = [a = 1., t = 3.])
> 

``

Download cubic.mw

The Newton polygon...

_Maxim_ 729
series rootof
November 08 2017
2 4

Inspired by this question. There is a simple iterated procedure that can generate the Puiseux expansion.

Note the use of convert(..., rational, exact) to preserve the digits of the floating-point numbers.

> 

(proc () global COcrit, P; Digits := 30; COcrit := 1/convert(2*0.115e-12, rational, exact); P := (`@`(`@`(rcurry(unapply, [NO2, O3]), numer), rcurry(convert, rational, exact)))(1.8027350279544425625*O3^8/10^49+(2.982942852010948125*CO/10^49+2.27611508925183215625*NO2/10^47+3.754849807690565625/10^37)*O3^7+(1.2339511797887015625*CO^2/10^49+2.4836920553034140625*CO/10^37+(-1)*5.9759150186390484375*NO2/10^35+6.3862221928648528125*NO2^2/10^46+1.88311680046014890625*CO*NO2/10^47+(-1)*9.69704144499346875/10^24)*O3^6+((-1)*1.71193039685859375*CO^2/10^37+(-1)*5.7098636544065625*CO/10^24+(-1)*1.277325786277575*NO2/10^21+(-1)*1.0801570017944671875*NO2^2/10^32+(-1)*2.9081778565815421875*CO*NO2/10^34+(-1)*3.66453227489203125/10^11)*O3^5+(1.9152220505625*CO^2/10^25+8.1035862984375*CO/10^13+1.40846609345625*NO2/10^10+(-1)*5.19285353257125*NO2^2/10^21+(-1)*1.55036925507375*CO*NO2/10^21+(-1)*1.844759695120875*CO^2*NO2/10^34+(-1)*1.876842472427325*CO*NO2^2/10^32-7.201634275625)*O3^4+(1.793091274970625*NO2^2/10^7+8618.14231275*NO2+(-1)*2.298266460675*CO^2*NO2/10^22+(-1)*9.5902239009375*CO*NO2/10^10+(-1)*1.685705248305*CO*NO2^2/10^20+9.2666503797075*NO2^3/10^19)*O3^3+((-1)*2.5638555726*10^6*NO2^2+6.894799382025*NO2^2*CO^2/10^20+2.7563954788125*CO*NO2^2/10^7+3.5073544682475*NO2^3*CO/10^18+(-1)*0.604340578881e-4*NO2^3+(-1)*3.5683519372605*NO2^4/10^16)*O3^2+((-1)*8.75499651*10^6*NO2^3+0.482686765875e-5*NO2^3*CO+(-1)*0.98216263665e-4*NO2^4)*O3+5.4066609*10^7*NO2^4) end proc)()

We're interested in the asymptotics of RootOf(P(NO2, _Z)) for small NO2.

> 

plot3d(RootOf(P(NO2, _Z)), CO = .5*COcrit .. 2*COcrit, NO2 = 10^(-3) .. 10^(-6))

Start with P(NO2, Z) and look for the expansion for Z when NO2 is small.

Expand P(NO2, Z) into monomials and convert each monomial a*Typesetting:-mi("NO2", italic = "true", mathvariant = "italic")^Typesetting:-mi("p", italic = "true", mathvariant = "italic")*Typesetting:-mi("Z", italic = "true", mathvariant = "italic")^Typesetting:-mi("q", italic = "true", mathvariant = "italic") into the point [p, q].

> 

([op])(collect(P(NO2, Z), [NO2, Z], distributed, normal)); map(proc (e) options operator, arrow; `~`[degree](e, [NO2, Z]) end proc, %)

[[4, 2], [4, 1], [4, 0], [3, 3], [3, 2], [3, 1], [2, 6], [2, 5], [2, 4], [2, 3], [2, 2], [1, 7], [1, 6], [1, 5], [1, 4], [1, 3], [0, 8], [0, 7], [0, 6], [0, 5], [0, 4]]

 (1)

The Newton polygon is the convex hull of the points.

> 

newtonPoly := proc (P, xy) local terms, pts; terms := ([op])(collect(P, xy, distributed, normal)); pts := map(proc (e) options operator, arrow; `~`[degree](e, xy) end proc, terms); plots:-display(plottools:-polygon(simplex:-convexhull(pts)), plottools:-point(pts), symbol = solidcircle, symbolsize = 15, style = line, thickness = 3, color = [magenta, black], axis = [gridlines = spacing(1)]) end proc
> 

newtonPoly(P(NO2, Z), [NO2, Z])

 

The side closest to the origin will correspond to the asymptotics for small NO2.

Sum the corresponding monomials and solve for Z.

> 

sideAsympt := proc (P, xy, pts) local lead; lead := add(coeff(coeff(P, xy[1], p[1]), xy[2], p[2])*xy[1]^p[1]*xy[2]^p[2], `in`(p, pts)); ([solve])(lead, xy[2]) end proc
> 

sideAsympt(P(NO2, Z), [NO2, Z], [[0, 4], [1, 3], [2, 2], [3, 1], [4, 0]])

[600*NO2, 600*NO2, -57000000000000000*NO2/(4071*CO-17795000000000000), -228000000000000000*NO2/(4071*CO+35020000000000000)]

 (2)

We have obtained the first term in the expansion. If CO>COcrit, the smallest positive root is the one asymptotic to 600*NO2. That will be the value of RootOf(P(NO2, _Z)).

The expansion can be continued in the same manner.

> 

newtonPoly(P(NO2, NO2*(600+Z)), [NO2, Z])

 

In this case we don't have to choose between sides with different slopes. (Compare to "P(NO2,600*NO2+Z).)"

> 

sideAsympt(P(NO2, NO2*(600+Z)), [NO2, Z], [[4, 2], [5, 0]]); `assuming`([simplify(`assuming`([simplify(%)], [NO2 > 0]))], [CO > COcrit])

[(531/430000)*354^(1/2)*NO2^(1/2)*(200000000000000+23*CO)^(1/2)/(23*CO-100000000000000)^(1/2), -(531/430000)*354^(1/2)*NO2^(1/2)*(200000000000000+23*CO)^(1/2)/(23*CO-100000000000000)^(1/2)]

 (3)

We get `~`[Z]*NO2^(1/2), so the next order in the expansion is NO2^(3/2).

Again, if we want the solution that corresponds to RootOf(P(NO2, _Z)), we should take the negative term, as the principal root will be the smaller one.

Find one more term. To avoid fractional powers, take NO2r = NO2^(1/2).

> 

approx := 600*NO2r^2+NO2r^3*(-(531/430000)*sqrt(354)*sqrt(200000000000000+23*CO)/sqrt(23*CO-100000000000000)+Z)
> 

newtonPoly(P(NO2r^2, approx), [NO2r, Z])

 
> 

sideAsympt(P(NO2r^2, approx), [NO2r, Z], [[10, 1], [11, 0]])

[(4779/36980000000000)*(4341503*CO^2+29548445000000000000*CO-209981000000000000000000000000000)*NO2r/(529*CO^2-4600000000000000*CO+10000000000000000000000000000)]

 (4)

We get `~`[Z]*NO2r, so we have obtained the coefficient at NO2r^3*NO2r = NO2^2.

All of this can be done using the algcurves:-puiseux command. The only difficulty is the unwieldy expressions that algcurves:-puiseux will generate for P.

Let's also find the NO2^2 term by the method of dominant balance. Suppose that we don't know p yet and are looking for the term alpha*NO2^p = alpha*NO2p.

> 

approx := 600*NO2r^2-(531/430000)*sqrt(354)*sqrt(200000000000000+23*CO)*NO2r^3/sqrt(23*CO-100000000000000)+alpha*NO2p
> 

terms := (`@`([op], collect))(P(NO2r^2, approx), [NO2r, NO2p], distributed, normal)

Find all exponents of NO2, coming from NO2r and from NO2p.

> 

pows := map(proc (t) options operator, arrow; (1/2)*degree(t, NO2r)+p*degree(t, NO2p) end proc, terms)
> 

plot(pows, p = 1 .. 2, view = 4 .. 6, annotation = pows)

 

At p=2, the two dominant terms NO2^(11/2) and NO2^(7/2+p) can balance each other.

> 

(`@`(normal, coeff))(subs(NO2p = NO2r^4, P(NO2r^2, approx)), NO2r, 11); ([solve])(%, alpha)

[(4779/36980000000000)*(4341503*CO^2+29548445000000000000*CO-209981000000000000000000000000000)/(529*CO^2-4600000000000000*CO+10000000000000000000000000000)]

 (5)
> 

``

Download np.mw

Rotations of an icosahedron...

_Maxim_ 729
rotation group icosahedron
September 07 2017
1 0

This is a toy example illustrating three possible ways to define a group.

An icosahedron:

with(GroupTheory): with(Student[LinearAlgebra]): with(geom3d):

icosahedron(ii): vv := vertices(ii):

PLOT3D(POLYGONS(op(evalf(faces(ii))), TRANSPARENCY(.75)),
  op(zip(TEXT, 1.1*evalf(vv), [`$`(1..12)])), AXESSTYLE(NONE));

The group of rotations is generated by the rotation around the diagonal (1,4) and the rotation around the line joining the midpoints of the edges (1,2) and (3,4), by the angles 2*Pi/5 and Pi respectively.

Define the group by how the two rotations permute the 6 main diagonals of the icosahedron:

gr := PermutationGroup({[[2, 5, 3, 4, 6]], [[1, 2], [5, 6]]});

IdentifySmallGroup(gr);
                             60, 5

Or define the group by the relations between the two rotations:

gr2 := FPGroup([a, b], [[a$5], [b$2], [a, b, a, b, a, b]]);

IdentifySmallGroup(gr2);
                             60, 5

Finally, define a group with elements that are rotation matrices:

m1 := RotationMatrix(2*Pi/5, Vector(op(4, vv) - op(1, vv))):
m2 := RotationMatrix(Pi, Vector(add(op(3..4, vv) - op(1..2, vv)))):
m1, m2 := op(evala(convert([m1, m2], radical))):

gr3 := CustomGroup([m1, m2], `.` = evala@`.`, `/` = evala@rcurry(`^`, -1), `=` = Equal);

IdentifySmallGroup(gr3);
                             60, 5

AreIsomorphic(SmallGroup(60, 5), AlternatingGroup(5));
                              true

One question is how to find out that the group is actually A5 without looking up the group (60, 5) elsewhere.

Also, it doesn't matter for this example, but adding `1`=IdentityMatrix(3) to the CustomGroup definition gives an error.

 

Directional derivatives of complex funct...

_Maxim_ 729
derivative complex directional-derivatives
September 06 2017
1 3

It can be interesting to consider a directional derivative of f(z) in the direction w:

%ddiff(f(z), z, w) = %limit((f(z + w*h) - f(z))/(w*h), h = 0);
ddiff := proc (fz0, z, dir) local rule, fz, dfz, ans;
  dfz := %ddiff(fz, z, dir)*dir;
  rule :=
   [abs(1, fz::anything) = (conjugate(fz)*dfz + fz*conjugate(dfz))/(2*abs(fz)*dfz),
    signum(1, fz::anything) = (conjugate(fz)*dfz - fz*conjugate(dfz))/(2*abs(fz)*conjugate(fz)*dfz)];
  ans := applyrule(rule, diff(fz0, z));
  ans := value(ans);
  ans := [ans, op(convert~(ans, [abs, argument, Re, Im, signum, conjugate]))];
  op(1, sort(simplify(ans, size), length)) end proc;

For analytic functions the derivative is the same in any direction:

ddiff(sqrt(Re(f(z))^2+Im(f(z))^2)*exp(I*argument(f(z))), z, w); # f(z) in disguise
                             d      
                            --- f(z)
                             dz

For non-analytic functions that's no longer the case:

ddiff(conjugate(z), z, w);
                               1     
                           ----------
                                    2
                           signum(w) 

ddiff(conjugate(f(z)), z, 1+I);
                             ________
                              d      
                          -I --- f(z)
                              dz     

ddiff(abs(z), z, z);
                               1    
                           ---------
                           signum(z)

ddiff(ln(abs(z)), z, z);
                               1
                               -
                               z

Some of those derivates have simple geometric interpretations: the derivative of argument(z) in the direction z is zero, since argument(z) doesn't change when moving in the direction z from the point z; the derivative of abs(z) in the direction I*z is zero, since the direction is tangent to the circle abs(z)=constant; since signum(z) is a function of argument(z) only, its derivative in the direction z is zero as well.

Interestingly, the derivative taken twice in the direction z is zero for each of the six basic functions:

map(fz -> ddiff(ddiff(fz, z, z), z, z), [abs, argument, Re, Im, signum, conjugate](z));
                       [0, 0, 0, 0, 0, 0]

Does that have some simple geometric interpretation as well?

 

1 2 Page 1 of 2