Markiyan Hirnyk

Markiyan Hirnyk
10 years, 190 days


These are answers submitted by Markiyan Hirnyk

@Carl Love 

>P := varepsilon^3+varepsilon^2+varepsilon+1;
>convert(map(c ->piecewise(2 <= degree(c), 0, c), [op(P)]), `+`);
1 + varepsilon
>selectremove(c -> 2 <= degree(c) , P)[2];
1 + varepsilon

PS. I think that was asked and answered a lot.

In order to obtain a concrete solution, one should evaluate parameters. For example,


restart; x := -(1/2*(-lambda^2*sigma^2-2*lambda*mu+2*RootOf(-exp(_Z)*erf((1/4)*sqrt(2)*(lambda^2*sigma^2+2*_Z)/(lambda*sigma))+exp(_Z)+erf((1/4)*sqrt(2)*(-lambda^2*sigma^2+2*_Z)/(lambda*sigma))+2*y-1)))/lambda
``

evalf(allvalues(eval(x, [lambda = 2, sigma = 1, mu = 2, y = 3/10])))

1.903267698

(1)

allvalues(eval(x, [lambda = 2, sigma = 1, mu = 2, y = .3]))

1.903267698

(2)

``

If values of some parameters are given as floats, then evalf is not necessary.

Download evaluating.mw

This can be done as follows, making use of the shadebetween command.

plots:-shadebetween(-y^2, x^2, y = 0 .. x, x = 0 .. 1, axes = frame, filled, style = surface);

 

plots:-shadebetween(-y^2, x^2, y = 0 .. x, x = 0 .. 1, axes = frame, filled, style = surface)

 

NULL

Edit. Ranges.

Download shadebetween1.mw

 

By

plots:-implicitplot(abs(eval(z*exp(1-z), z = x+I*y)) = 1, x = -2 .. 2, y = -2 .. 2, gridrefine = 3)

plots:-implicitplot(abs(eval(z*exp(1-z), z = x+I*y)) = 1, x = -2 .. 2, y = -2 .. 2, gridrefine = 3)``

 

``

 

Download implicitplot.mw

Because the improper integral converges, its value equals its Cauchy principal value:

>int(exp(-I*k*x)/cosh(x), x = -infinity .. infinity, CPV);
Pi*sech((1/2)*k*Pi)

To install a package, proceed as follows.

(1) Create or identify a directory (folder) where you will keep the mla file(s) and the help file. Copy the mla and help files there.

(2) Open Maple, at the prompt type libname. This will display where maple currently looks for packages. You need to modify libname to include the path to the directory used in (1). This can be done by assigning libname := newpath, libname;.

(3) If you are installing the DifferentialGeometry package: Execute with(DifferentialGeometry) to load the package. At the prompt execute DGbuild; and verify that the build number matches the build displayed on this website. This step can be ignored if you are installing DGApplications.

(4) Step (2) can be automatically executed at the start of every Maple session by adding the libname assignment to your maple init file. In Maple Help, see worksheet,reference,initialization. Or go to the URL http://www.maplesoft.com/support/help/Maple/view.aspx?path=worksheet/reference/initialization .

This is written here http://digitalcommons.usu.edu/dg_downloads/5/

Partial answer

April 21 2016 Markiyan Hirnyk 6919

can be done as follows. We find the rational parametrization in t  of the parabola and note that it is enough t=9*n with an arbitrary integer n to this end.

 

with(algcurves):````

L := parametrization(4*x^2-12*x*y+9*y^2-7*x+9*y+369, x, y, t)

[(1/9)*t^2-(1/3)*t+367, (2/27)*t^2-(1/9)*t+244]

(1)

x := L[1]

(1/9)*t^2-(1/3)*t+367

(2)

y := L[2]

(2/27)*t^2-(1/9)*t+244

(3)

eval([x, y], t = 9*n)

[9*n^2-3*n+367, 6*n^2-n+244]

(4)

``

``

Unfortunately, the command

>isolve({(1/9)*t^2-(1/3)*t+367 = n, (2/27)*t^2-(1/9)*t+244 = k});

produces an incorrect answer

{k = (2/19683)*(4779+(1/1162261467)*_Z1)^2+673/3-(1/282429536481)*_Z1, n = (1/6561)*(4779+(1/1162261467)*_Z1)^2+308-(1/94143178827)*_Z1, t = 177+(1/31381059609)*_Z1}

Mathematica finds the complete answer

(C[1] \[Element] Integers && x == 369 - 9 C[1] + 9 C[1]^2 &&
y == 246 - 7 C[1] + 6 C[1]^2) || (C[1] \[Element] Integers &&
x == 367 + 3 C[1] + 9 C[1]^2 && y == 244 + C[1] + 6 C[1]^2)

 

Download integer_points_on_parabola.mw

General case

April 12 2016 Markiyan Hirnyk 6919

Following Igor Hlivka, the two-dimensional normal distribution is created and its characteristics are found. Specifying the parameters, a sample of the size 1000 is created and plotted. See section

 Drawing values from the distribution

in the Wiki article for math background.



restart; with(Statistics); with(LinearAlgebra); with(plots)

h := proc (i) options operator, arrow; x[i] end proc; g := proc (i) options operator, arrow; mu[i] end proc; w := proc (i) options operator, arrow; sigma[i] end proc; N := 2; X := Vector(N, h); M := Vector(N, g); S := Vector(N, w); CorrMat := Matrix([seq([seq(`if`(i = j, 1, rho), j = 1 .. N)], i = 1 .. N)]); C1 := simplify(S.Transpose(S)); CoVar := Zip(`*`, C1, CorrMat)

X := Vector(2, {(1) = x[1], (2) = x[2]})

 

M := Vector(2, {(1) = mu[1], (2) = mu[2]})

 

S := Vector(2, {(1) = sigma[1], (2) = sigma[2]})

 

CorrMat := Matrix(2, 2, {(1, 1) = 1, (1, 2) = rho, (2, 1) = rho, (2, 2) = 1})

 

CoVar := Matrix(2, 2, {(1, 1) = sigma[1]^2, (1, 2) = sigma[1]*sigma[2]*rho, (2, 1) = sigma[1]*sigma[2]*rho, (2, 2) = sigma[2]^2})

(1)

p1 := exp(-(1/2)*Typesetting:-delayDotProduct(Typesetting:-delayDotProduct(LinearAlgebra:-Transpose(X-M), 1/CoVar), X-M))/((2*Pi)^((1/2)*N)*sqrt(Determinant(CoVar))):

bipdf := simplify(p1, symbolic);

(1/2)*exp((1/2)*((x[2]-mu[2])^2*sigma[1]^2-2*rho*sigma[2]*(x[2]-mu[2])*(x[1]-mu[1])*sigma[1]+sigma[2]^2*(x[1]-mu[1])^2)/(sigma[1]^2*sigma[2]^2*(rho-1)*(rho+1)))/(sigma[1]*sigma[2]*(-rho^2+1)^(1/2)*Pi)

(2)

MarDist[X[1]] := `assuming`([int(bipdf, x[2] = -infinity .. infinity)], [-1 < rho and rho < 1, sigma[1] > 0, sigma[2] > 0]);

(1/2)*2^(1/2)*exp(-(1/2)*(mu[1]^2-2*mu[1]*x[1]+x[1]^2)/sigma[1]^2)/(Pi^(1/2)*sigma[1])

(3)

MarDist[X[2]] := `assuming`([int(bipdf, x[1] = -infinity .. infinity)], [-1 < rho and rho < 1, sigma[1] > 0, sigma[2] > 0]);

(1/2)*2^(1/2)*exp(-(1/2)*(mu[2]^2-2*mu[2]*x[2]+x[2]^2)/sigma[2]^2)/(sigma[2]*Pi^(1/2))

(4)

Conditional PDF of  "X|Y = y"          "f[X|Y]= (f[X,Y](x,y))/(f[Y](y))"
CondDist[XY] := simplify(bipdf/MarDist[X[2]], symbolic)

(1/2)*2^(1/2)*exp((1/2)*(sigma[1]*(x[2]-mu[2])*rho-sigma[2]*(x[1]-3))^2/(sigma[1]^2*sigma[2]^2*(rho-1)*(rho+1)))/((-rho^2+1)^(1/2)*Pi^(1/2)*sigma[1])

(5)

mu[1] := 3:``

mu[2] := 15: 

sigma[1] := 1:

rho := -.9:NULL

CoVar

Matrix([[sigma[1]^2, sigma[1]*sigma[2]*rho], [sigma[1]*sigma[2]*rho, sigma[2]^2]])

(6)

A := LUDecomposition(CoVar, method = Cholesky)

Matrix([[1, 0], [-1.8, .8717797887]])

(7)

``

MU := Vector[row]([mu[1], mu[2]])

Vector[row]([3, 15])

(8)

sx := Sample(Normal(0, 1), 1000)

Vector(4, {(1) = ` 1 .. 1000 `*Vector[row], (2) = `Data Type: `*float[8], (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})

(9)

sy := Sample(Normal(0, 1), 1000)

Vector(4, {(1) = ` 1 .. 1000 `*Vector[row], (2) = `Data Type: `*float[8], (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})

(10)

sample := zip(proc (x, y) options operator, arrow; [x, y] end proc, sx, sy)

Vector(4, {(1) = ` 1 .. 1000 `*Vector[row], (2) = `Data Type: `*anything, (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})

(11)

B := map(proc (c) options operator, arrow; MU+convert(Multiply(A, convert(c, Vector[column])), Vector[row]) end proc, sample)

Vector(4, {(1) = ` 1 .. 1000 `*Vector[row], (2) = `Data Type: `*anything, (3) = `Storage: `*rectangular, (4) = `Order: `*Fortran_order})

(12)
pointplot(map(proc (c) options operator, arrow; convert(c, list) end proc, B), scaling = constrained);

 

pointplot(map(proc (c) options operator, arrow; convert(c, list) end proc, B), scaling = constrained)

 

``



Download _adjusted_sample_1000.mw

 

 

 

 



 

It appears there is a lot of such polynomials :

 

bc := {((D@@1)(g))(1)+(1/2)*((D@@1)(s))(1)^2 = 0, g(0) = 0, s(0) = 0, ((D@@1)(s))(0) = 0, ((D@@2)(s))(1) = 0, ((D@@3)(s))(1) = 0}

{(D(g))(1)+(1/2)*(D(s))(1)^2 = 0, g(0) = 0, s(0) = 0, (D(s))(0) = 0, ((D@@2)(s))(1) = 0, ((D@@3)(s))(1) = 0}

(1)

dsys := {diff(g(x), `$`(x, 4)) = 0, diff(s(x), `$`(x, 5)) = 0}

{diff(diff(diff(diff(diff(s(x), x), x), x), x), x) = 0, diff(diff(diff(diff(g(x), x), x), x), x) = 0}

(2)

sol := dsolve(`union`(dsys, bc), {g(x), s(x)})

{g(x) = (1/6)*(-_C1^2-2*_C3-2*_C4)*x^3+(1/2)*_C3*x^2+_C4*x, s(x) = (1/4)*_C1*x^4-_C1*x^3+(3/2)*_C1*x^2}

(3)

eval(sol, [_C1 = 1, _C3 = 1, _C4 = -1])

{g(x) = -(1/6)*x^3+(1/2)*x^2-x, s(x) = (1/4)*x^4-x^3+(3/2)*x^2}

(4)

``

 

Download polynomials.mw

 

Another way

April 06 2016 Markiyan Hirnyk 6919

 

match((a+b+1)*t = -a+b+1, t, 's')

true

(1)

s

{a = 0, b = -1}

(2)

``


Look in ?match for info.

Download another_way.mw

PS. The match command works in several dimensions too, whereas the solve,identity command works in one-dimensional case only.

Specifying the parameters, it can be solved with the DirectSearch package which should be downloaded from http://www.maplesoft.com/applications/view.aspx?SID=101333 and installed in your Maple >=13.

 

restart; A := eval(-(((beta*eta^2-(1/2)*beta+1)*p^2-(1/2)*beta^3)*KummerM(((-beta+2)*p^2-beta^3)/(4*p^2), 1, beta*eta^2)+(1/2)*KummerM(((-beta+6)*p^2-beta^3)/(4*p^2), 1, beta*eta^2)*((beta-2)*p^2+beta^3))*exp(-(1/2)*beta*eta^2)/(p^2*beta) = 0, [p = 1, eta = 1]);

-(((1/2)*beta+1-(1/2)*beta^3)*KummerM(-(1/4)*beta^3-(1/4)*beta+1/2, 1, beta)+(1/2)*KummerM(-(1/4)*beta^3-(1/4)*beta+3/2, 1, beta)*(beta^3+beta-2))*exp(-(1/2)*beta)/beta = 0

(1)

DirectSearch:-SolveEquations([Re(eval(A, beta = a+I*b)), Im(eval(A, beta = a+I*b))], {a = -5 .. 5, b = -5 .. 5}, AllSolutions, number = 300,solutions=10);

Matrix(10, 4, {(1, 1) = 0.1948977161e-26, (1, 2) = Vector(2, {(1) = HFloat(-1.0451202215029108e-14), (2) = HFloat(-4.289230156047799e-14)}), (1, 3) = [a = -4.39157583457543, b = 0.6122407479e-14], (1, 4) = 80, (2, 1) = 0.2468683635e-26, (2, 2) = Vector(2, {(1) = HFloat(4.76220423508406e-14), (2) = HFloat(-1.4171263788954476e-14)}), (2, 3) = [a = 0.9353300437e-13, b = .705280177757183], (2, 4) = 80, (3, 1) = 0.1107221646e-25, (3, 2) = Vector(2, {(1) = HFloat(1.0344270982184144e-13), (2) = HFloat(1.928269298431157e-14)}), (3, 3) = [a = 0.2043030414e-12, b = -.705280177757170], (3, 4) = 112, (4, 1) = 0.1498604938e-25, (4, 2) = Vector(2, {(1) = HFloat(-2.111504886120853e-14), (2) = HFloat(-1.2058276865205994e-13)}), (4, 3) = [a = 4.01763280323693, b = -0.1881620268e-13], (4, 4) = 82, (5, 1) = 0.2098144474e-25, (5, 2) = Vector(2, {(1) = HFloat(1.3962856649975636e-13), (2) = HFloat(-3.853969581723759e-14)}), (5, 3) = [a = -0.2730763486e-13, b = -3.24247382580224], (5, 4) = 87, (6, 1) = 0.5674455128e-25, (6, 2) = Vector(2, {(1) = HFloat(-7.700514667200426e-14), (2) = HFloat(2.2542129150220777e-13)}), (6, 3) = [a = 1.87588952396935, b = 0.7554868781e-13], (6, 4) = 81, (7, 1) = 0.8186683329e-25, (7, 2) = Vector(2, {(1) = HFloat(2.5763416746195455e-13), (2) = HFloat(-1.244647301461954e-13)}), (7, 3) = [a = 0.5969031808e-13, b = 2.71279028626207], (7, 4) = 54, (8, 1) = 0.9325892563e-25, (8, 2) = Vector(2, {(1) = HFloat(5.836258384818908e-14), (2) = HFloat(2.9975445690543977e-13)}), (8, 3) = [a = -1.87588952396936, b = 0.1004610334e-12], (8, 4) = 80, (9, 1) = 0.3667927203e-24, (9, 2) = Vector(2, {(1) = HFloat(-1.582465938541118e-13), (2) = HFloat(-5.845945054753967e-13)}), (9, 3) = [a = 0.4867515661e-13, b = -2.04633245164591], (9, 4) = 74, (10, 1) = 0.3928454083e-24, (10, 2) = Vector(2, {(1) = HFloat(-2.3898361217675557e-13), (2) = HFloat(5.794240600653144e-13)}), (10, 3) = [a = 0.4650368771e-13, b = 3.24247382580213], (10, 4) = 89})

(2)

``

Ten roots out of 29 roots in the square are found. In order to get more solutions, one should increase the value of the solutions option.

Download roots.mw

One should exploit the Newton-Leibnitz formula:


J := `assuming`([int(t^n*ln(t)^n, t = 0 .. x)], [n::posint, x > 0]);

-(n^2*GAMMA(n, -(n+1)*ln(x))*(-n-1)^(-n-1)+n*GAMMA(n, -(n+1)*ln(x))*(-n-1)^(-n-1)-x*(ln(x)*x)^n)/(n+1)

(1)

`assuming`([simplify(diff(J, x)-x^n*ln(x)^n)], [n::posint, x > 0])

0

(2)

``


Download direct_calculation.mw

Improvement

April 02 2016 Markiyan Hirnyk 6919

In fact, this is an improvement of the sum command: see ?sum for info. There is  _EnvFormal option to this end. For example,

 

restart; S := convert(1/(3*x+2), FPS, x)

Sum((1/2)*2^(-k)*(-3)^k*x^k, k = 0 .. infinity)

(1)

value(S)

sum((1/2)*2^(-k)*(-3)^k*x^k, k = 0 .. infinity)

(2)

_EnvFormal := true

true

(3)

value(S)

1/(3*x+2)

(4)

``

 

Download improvement.mw

Another way is 

>restart;sum((1/2)*2^(-k)*(-3)^k*x^k, k = 0 .. infinity)assuming abs(x)<1/ 6;

1/(3*x+2)

The DirectSearch package should be downloaded from http://www.maplesoft.com/applications/view.aspx?SID=101333 and installed in your Maple (>=13). Here is my solution done with the DirectSearch. A more powerful comp is required to obtain more points by increasing the number option.

 

restart; f1 := cos(t1)+1.35*cos(p1)-x; f2 := sin(t1)+1.35*sin(p1)-y; f3 := cos(t2)+1.35*cos(p2)-x+1.15; f4 := sin(t2)+1.35*sin(p2)-y; f5 := cos(t3)+1.35*cos(p3)-x+.575; f6 := sin(t3)+1.35*sin(p3)-y+.995; f7 := cos(t1)*cos(t2)*cos(t3)*sin(p1)*sin(p2)*sin(p3)-cos(t1)*cos(t2)*sin(t3)*sin(p1)*sin(p2)*cos(p3)-cos(t1)*sin(t2)*cos(t3)*sin(p1)*cos(p2)*sin(p3)+cos(t1)*sin(t2)*sin(t3)*sin(p1)*cos(p2)*cos(p3)-sin(t1)*cos(t2)*cos(t3)*cos(p1)*sin(p2)*sin(p3)+sin(t1)*cos(t2)*sin(t3)*cos(p1)*sin(p2)*cos(p3)+sin(t1)*sin(t2)*cos(t3)*cos(p1)*cos(p2)*sin(p3)-sin(t1)*sin(t2)*sin(t3)*cos(p1)*cos(p2)*cos(p3); f8 := cos(t1)*sin(p1)-sin(t1)*cos(p1) >= 0; f9 := cos(t2)*sin(p2)-sin(t2)*cos(p2) >= 0; f10 := cos(t3)*sin(p3)-sin(t3)*cos(p3) >= 0

sol := DirectSearch:-SolveEquations([seq(f || j, j = 1 .. 7)], {f10, f8, f9, p1 = -Pi .. Pi, p2 = -Pi .. Pi, p3 = -Pi .. Pi, t1 = -Pi .. Pi, t2 = -Pi .. Pi, t3 = -Pi .. Pi, x = -10 .. 10, y = -10 .. 10}, AllSolutions, method = quadratic, number = 200);

RTABLE(18446744074192502782, anything, Matrix, rectangular, Fortran_order, [], 2, 1 .. 182, 1 .. 4)

(1)

sol[1];

Vector[row](4, {(1) = 0.1071032753e-27, (2) = Vector(7, {(1) = HFloat(-2.55351295663786e-15), (2) = HFloat(6.661338147750939e-16), (3) = HFloat(4.884981308350689e-15), (4) = HFloat(-3.1086244689504383e-15), (5) = HFloat(7.105427357601002e-15), (6) = HFloat(3.9968028886505635e-15), (7) = HFloat(3.885780586188048e-16)}), (3) = [p1 = 2.81233098296781, p2 = 2.54814516295493, p3 = -2.62961132730027, t1 = 1.07095279467161, t2 = 2.54814516295493, t3 = 1.76837387160141, x = -.798191584527774, y = 1.31417257237213], (4) = 1032})

(2)

``

 

 

plots:-pointplot3d([seq([rhs(sol(1 .. 53, 3)[j][7]), rhs(sol(1 .. 53, 3)[j][8]), rhs(sol(1 .. 53, 3)[j][4])], j = 1 .. 53)], axes = frame)

 

``

 

Download By_DS.mw

 

 

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