sand15

1007 Reputation

15 Badges

10 years, 305 days

MaplePrimes Activity


These are answers submitted by sand15

restart

with(Statistics):

f[1] := t -> piecewise(t <= 0, 0, 0 < t and t < 2, 1/(Pi*sqrt(1-(1-t)^2)), t >= 2, 0);
f[2] := t -> piecewise(t <= 0, 0, 0 < t and t < 1, 2*t, t >= 1, 0);

proc (t) options operator, arrow; piecewise(t <= 0, 0, 0 < t and t < 2, 1/(Pi*sqrt(1-(1-t)^2)), 2 <= t, 0) end proc

 

proc (t) options operator, arrow; piecewise(t <= 0, 0, 0 < t and t < 1, 2*t, 1 <= t, 0) end proc

(1)

Dist1 := Distribution(PDF = (t -> f[1](t))):
Dist2 := Distribution(PDF = (t -> f[2](t))):

X := RandomVariable(Dist1):
Y := RandomVariable(Dist2):

Z := X*Y:

`PDF(Z)` := PDF(Z, t);
`CDF(Z)` := CDF(Z, t);

`PDF(Z)` := piecewise(t < 0, 0, t <= 2, -(2/3)*(t^2-t-2)/(sqrt(t)*sqrt(2-t)*Pi), 2 < t, 0)

 

piecewise(t <= 0, 0, t < 2, (1/6)*(2*t^(3/2)*(2-t)^(1/2)*(-t*(t-2))^(1/2)+2*t^(1/2)*(2-t)^(1/2)*(-t*(t-2))^(1/2)+6*t^(1/2)*(2-t)^(1/2)*arcsin(-1+t)+3*Pi*(-t*(t-2))^(1/2))/(Pi*(-t*(t-2))^(1/2)), 2 <= t, 1)

(2)

plot(`PDF(Z)`, t=0..2)

 
 

 

Download Inquiry_sand15.mw

It doesn't matter what name the integral has, Mellin or anything else
Let X and two independent random variables whose supports belong to [0, +oo) and with respective pdf fX and fY  then the pdf fZ of Z=XY is given by 


This formula can be easily obtained by reasoning in terms of probability and not of densities of probability. Indeed Prob( Z < z ) = Prob( XY < z ) = Prob( X < x & Y < z/x)

A useful trick if you don't know about this formula  consists in defining two auxiliary random variables A=ln( X ) and B=ln( Y ) and to observe that the pdf of C=ln( Z ) is the convolution product of the pdfs of  A and B.
Once you know the pdf of C it remains to compute the one of exp( C ) .

Unless Maple is wrong or if rewritting the system in a more concise form removes some simplifications or factorizations, it doesn't seem a solution exists:

restart

K2 := -(Pn-Pr)/(1-delta)+(-beta*((-2*delta+2)*tau^2+((-Pn+Pr-delta+1)*Cr+(-2+2*delta)*(w+i2))*tau+Cr*i2*(-1+delta+Pn-Pr))*upsilon/((4*delta-4)*tau+Cr*(-1+delta+Pn-Pr))+Pr)/delta+(Pr-w-Cr)*beta*(-2+2*delta)*tau*upsilon/(((4*delta-4)*tau+Cr*(-1+delta+Pn-Pr))*delta) = 0:

K3 := 1-(Pn-Pr)/(1-delta)-(Pn-Cn)/(1-delta)+(Pr-w-Cr)*(1/(1-delta)-(-beta*(Cr*i2-Cr*tau)*upsilon/((4*delta-4)*tau+Cr*(-1+delta+Pn-Pr))+beta*((-2*delta+2)*tau^2+((-Pn+Pr-delta+1)*Cr+(-2+2*delta)*(w+i2))*tau+Cr*i2*(-1+delta+Pn-Pr))*upsilon*Cr/((4*delta-4)*tau+Cr*(-1+delta+Pn-Pr))^2)/delta)+Ce*rho0/(1-delta) = 0:

K4 := (Pn-Cn)/(1-delta)+(Pn-Pr)/(1-delta)-(-beta*((-2*delta+2)*tau^2+((-Pn+Pr-delta+1)*Cr+(-2+2*delta)*(w+i2))*tau+Cr*i2*(-1+delta+Pn-Pr))*upsilon/((4*delta-4)*tau+Cr*(-1+delta+Pn-Pr))+Pr)/delta+(Pr-w-Cr)*(-1/(1-delta)-(-beta*(-Cr*i2+Cr*tau)*upsilon/((4*delta-4)*tau+Cr*(-1+delta+Pn-Pr))-beta*((-2*delta+2)*tau^2+((-Pn+Pr-delta+1)*Cr+(-2+2*delta)*(w+i2))*tau+Cr*i2*(-1+delta+Pn-Pr))*upsilon*Cr/((4*delta-4)*tau+Cr*(-1+delta+Pn-Pr))^2+1)/delta)-Ce*rho0/(1-delta) = 0:

sol := timelimit(20, solve({K2, K3, K4}, {Pn, Pr, w}))

Error, (in LinearAlgebra:-Determinant) time expired

 


A good strategy almost always consits in rewritting the equations in a more concise form

rewrite := proc(e, v)
  local f1, f2, C, M, f3:
  # assuming that the roots of the numerator do not cancel out the denominator:
  f1 := numer(lhs(e));

  C  := [coeffs(f1, [Pn, Pr, w], 'M')];
  M  := [M];
  f2 := [seq([seq(degree(t, i), i in [Pn, Pr, w])], t in M)];
  f3 := add(v[f2[i]]*M[i], i=1..nops(M));
  return f3=0, [seq(v[f2[i]] = C[i], i=1..nops(M)) ];
end proc:

newK2, rw2 := rewrite(K2, U):
newK3, rw4 := rewrite(K3, V):
newK4, rw4 := rewrite(K4, W):

Example: K2 rewritten and the corresponding values of the U coefficients:

newK2;
print():
print~(rw2):

Pn^2*U[[2, 0, 0]]+Pn*Pr*U[[1, 1, 0]]+Pr^2*U[[0, 2, 0]]+Pn*U[[1, 0, 0]]+Pr*U[[0, 1, 0]]+w*U[[0, 0, 1]]+U[[0, 0, 0]] = 0

 

 

U[[2, 0, 0]] = Cr*delta

 

U[[1, 1, 0]] = -Cr*delta-Cr

 

U[[1, 0, 0]] = -Cr*beta*delta*i2*upsilon+Cr*beta*delta*tau*upsilon+Cr*beta*i2*upsilon-Cr*beta*tau*upsilon+Cr*delta^2+4*delta^2*tau-Cr*delta-4*delta*tau

 

U[[0, 2, 0]] = Cr

 

U[[0, 1, 0]] = Cr*beta*delta*i2*upsilon-Cr*beta*delta*tau*upsilon+2*beta*delta^2*tau*upsilon-Cr*beta*i2*upsilon+Cr*beta*tau*upsilon-4*beta*delta*tau*upsilon+2*beta*tau*upsilon-Cr*delta-4*delta*tau+Cr+4*tau

 

U[[0, 0, 1]] = -4*beta*delta^2*tau*upsilon+8*beta*delta*tau*upsilon-4*beta*tau*upsilon

 

U[[0, 0, 0]] = -Cr*beta*delta^2*i2*upsilon-Cr*beta*delta^2*tau*upsilon-2*beta*delta^2*i2*tau*upsilon+2*beta*delta^2*tau^2*upsilon+2*Cr*beta*delta*i2*upsilon+2*Cr*beta*delta*tau*upsilon+4*beta*delta*i2*tau*upsilon-4*beta*delta*tau^2*upsilon-Cr*beta*i2*upsilon-Cr*beta*tau*upsilon-2*beta*i2*tau*upsilon+2*beta*tau^2*upsilon

(1)

Here is the rewritten system

sys := {newK2, newK3, newK4}:
for e in sys do e; print(); end do;
 

Pn^2*U[[2, 0, 0]]+Pn*Pr*U[[1, 1, 0]]+Pr^2*U[[0, 2, 0]]+Pn*U[[1, 0, 0]]+Pr*U[[0, 1, 0]]+w*U[[0, 0, 1]]+U[[0, 0, 0]] = 0

 

 

Pn^3*V[[3, 0, 0]]+Pn^2*Pr*V[[2, 1, 0]]+Pn^2*w*V[[2, 0, 1]]+Pn*Pr^2*V[[1, 2, 0]]+Pn*Pr*w*V[[1, 1, 1]]+Pr^3*V[[0, 3, 0]]+Pr^2*w*V[[0, 2, 1]]+Pn^2*V[[2, 0, 0]]+Pn*Pr*V[[1, 1, 0]]+Pn*w*V[[1, 0, 1]]+Pr^2*V[[0, 2, 0]]+Pr*w*V[[0, 1, 1]]+w^2*V[[0, 0, 2]]+Pn*V[[1, 0, 0]]+Pr*V[[0, 1, 0]]+w*V[[0, 0, 1]]+V[[0, 0, 0]] = 0

 

 

Pn^3*W[[3, 0, 0]]+Pn^2*Pr*W[[2, 1, 0]]+Pn^2*w*W[[2, 0, 1]]+Pn*Pr^2*W[[1, 2, 0]]+Pn*Pr*w*W[[1, 1, 1]]+Pr^3*W[[0, 3, 0]]+Pr^2*w*W[[0, 2, 1]]+Pn^2*W[[2, 0, 0]]+Pn*Pr*W[[1, 1, 0]]+Pn*w*W[[1, 0, 1]]+Pr^2*W[[0, 2, 0]]+Pr*w*W[[0, 1, 1]]+w^2*W[[0, 0, 2]]+Pn*W[[1, 0, 0]]+Pr*W[[0, 1, 0]]+w*W[[0, 0, 1]]+W[[0, 0, 0]] = 0

 

(2)

newK2 being linear wrt w just solve it wrt w:

infolevel[solve] := 4:

wsol_1 := solve(newK2, w);

Main: Entering solver with 1 equation in 1 variable

Transformer:   solving the uncoupled linear subsystem in w

Main: solving successful - now forming solutions
Main: Exiting solver returning 1 solution

 

-(Pn^2*U[[2, 0, 0]]+Pn*Pr*U[[1, 1, 0]]+Pr^2*U[[0, 2, 0]]+Pn*U[[1, 0, 0]]+Pr*U[[0, 1, 0]]+U[[0, 0, 0]])/U[[0, 0, 1]]

(3)

Form a 2-by-2 subsystem while plugging the solution above into {newK3 and newK4}

sys34 := eval({newK3, newK4}, w=wsol_1):

newnewK3, newrw3 := rewrite(sys34[1], X):
newnewK4, newrw4 := rewrite(sys34[2], Y):

newsys := {newnewK3, newnewK4}:
for e in newsys do e; print(); end do;

Pn^4*X[[4, 0, 0]]+Pn^3*Pr*X[[3, 1, 0]]+Pn^2*Pr^2*X[[2, 2, 0]]+Pn*Pr^3*X[[1, 3, 0]]+Pr^4*X[[0, 4, 0]]+Pn^3*X[[3, 0, 0]]+Pn^2*Pr*X[[2, 1, 0]]+Pn*Pr^2*X[[1, 2, 0]]+Pr^3*X[[0, 3, 0]]+Pn^2*X[[2, 0, 0]]+Pn*Pr*X[[1, 1, 0]]+Pr^2*X[[0, 2, 0]]+Pn*X[[1, 0, 0]]+Pr*X[[0, 1, 0]]+X[[0, 0, 0]] = 0

 

 

Pn^4*Y[[4, 0, 0]]+Pn^3*Pr*Y[[3, 1, 0]]+Pn^2*Pr^2*Y[[2, 2, 0]]+Pn*Pr^3*Y[[1, 3, 0]]+Pr^4*Y[[0, 4, 0]]+Pn^3*Y[[3, 0, 0]]+Pn^2*Pr*Y[[2, 1, 0]]+Pn*Pr^2*Y[[1, 2, 0]]+Pr^3*Y[[0, 3, 0]]+Pn^2*Y[[2, 0, 0]]+Pn*Pr*Y[[1, 1, 0]]+Pr^2*Y[[0, 2, 0]]+Pn*Y[[1, 0, 0]]+Pr*Y[[0, 1, 0]]+Y[[0, 0, 0]] = 0

 

(4)

Note that the lhs of each equations are complete polynomials of total degree 4 in the two indeterminates Pn and Pr

infolevel[solve] := 4:
sol := timelimit(20, solve(newsys, {Pn, Pr}));

Main: Entering solver with 2 equations in 2 variables
Main: attempting to solve as a linear system
Main: attempting to solve as a polynomial system
Main: polynomial system split into 1 parts under preprocessing
Main: trying resultant methods
PseudoResultant: 4093009 [2 7000140160 Pn] 2 2 649 0 3 0

PseudoResultant: factoring all equations length = 649
PseudoResultant: factoring done length = 649
LinearEq: computing pseudo sub-resultant in Pn of polynomials of length 322 and 322
LinearEq: computing subresultant determinant of dimension 6 and length 1512

LinearEq: subresultant determinant of length 92449 computed

LinearEq: computing pseudo sub-resultant in Pr of polynomials of length 322 and 322
LinearEq: computing subresultant determinant of dimension 6 and length 1512
LinearEq: subresultant determinant of length 92449 computed

PseudoResultant: computing pseudoresultant in Pn with lengths 322 and 322

Error, (in convert/sqrfree/sqrfree) time expired

 

timelimit(20, SolveTools:-PolynomialSystem(lhs~(newsys), {Pn, Pr}));

Main: polynomial system split into 1 parts under preprocessing

Main: trying resultant methods
PseudoResultant: 4445613 [2 7000142560 Pn] 2 2 809 0 3 0
PseudoResultant: factoring all equations length = 809

PseudoResultant: factoring done length = 809
LinearEq: computing pseudo sub-resultant in Pn of polynomials of length 402 and 402
LinearEq: computing subresultant determinant of dimension 6 and length 1850
LinearEq: subresultant determinant of length 127039 computed

LinearEq: computing pseudo sub-resultant in Pr of polynomials of length 402 and 402
LinearEq: computing subresultant determinant of dimension 6 and length 1850

LinearEq: subresultant determinant of length 127039 computed

PseudoResultant: computing pseudoresultant in Pn with lengths 402 and 402
SolutionsLost: setting solutions lost flag
PseudoResultant: 0 solutions found, now doing backsubstitution

SolutionsLost: setting solutions lost flag

 
 

 

Seems to be a dead line

Download Q_solve_sand15.mw

restart

aux := B12 = -6*(p1 + p2)/(p1 - p2)^2;

B12 = -6*(p1+p2)/(p1-p2)^2

(1)

F2  := (theta1*theta2*p1^2 + (-2*p2*theta1*theta2 - 6)*p1 + theta1*theta2*p2^2 - 6*p2)/(p1 - p2)^2;

(theta1*theta2*p1^2+(-2*p2*theta1*theta2-6)*p1+theta1*theta2*p2^2-6*p2)/(p1-p2)^2

(2)

collect(expand(F2), theta1):;
map(factor, %);

F2a := op(1, %) + simplify(op(2, %)+op(3, %));

theta1*theta2-6*p1/(p1-p2)^2-6*p2/(p1-p2)^2

 

theta1*theta2-6*(p1+p2)/(p1-p2)^2

(3)

F2b := eval(F2a, (rhs=lhs)(aux))

theta1*theta2+B12

(4)
 

 

Download rewrite.mw

Main remark:
Your model y = P(x)/Q(x) where P and Q are two polynomials of x can be linearized under a convenient
reparameterization: basically you discard the regressor x and introduce new ones which make the model linear).
Once done you can use  Statistics:-Fit instead of Statistics:-NonlinearFit (in fact Optimization:-LSSolve) and thus discard any convergence issue.
More of this you avoid getting a solution which depends on the initial values
         (you already noticed that the closer to the values used in the generative model, the better the results [which supposes this
          generative model is noiseless])

and you get access to all the output statisctics than  Statistics:-Fit offers)

Awaiting your data Ycg and Sx, here is a notional example (same model than yours but data produced by a noisy process) which explains how to linearize the model and go back to its initial nonlinear form.
 

restart

with(Statistics):

Model__ini := y = (p1*x^2 + p2*x + p3)/(q1*x + x^2 + q2)

y = (p1*x^2+p2*x+p3)/(q1*x+x^2+q2)

(1)


It is more practical do rewrite Model0 this way

Model__0 := y = (a[2]*x^2 + a[1]*x + a[0])/(b[2]*x^2 + b[1]*x + 1)

y = (x^2*a[2]+x*a[1]+a[0])/(x^2*b[2]+x*b[1]+1)

(2)


Assuming the denominator is never null in some range of x, this model can be linearized.
Indeed:

Model__1 := expand(lhs(Model__0)*denom(rhs(Model__0))) = numer(rhs(Model__0))

x^2*y*b[2]+x*y*b[1]+y = x^2*a[2]+x*a[1]+a[0]

(3)

VarChange := {x^2*y*b[2] = Z__21*b[2], x*y*b[1] = Z__11*b[1], x^2*a[2] = Z__20*a[2], x*a[1] = Z__10*a[1]}

{a[1]*x = Z__10*a[1], a[2]*x^2 = Z__20*a[2], x*y*b[1] = Z__11*b[1], x^2*y*b[2] = Z__21*b[2]}

(4)

Model__2 := subs(VarChange, Model__1)

Z__11*b[1]+Z__21*b[2]+y = Z__10*a[1]+Z__20*a[2]+a[0]

(5)

Model__lin := isolate(Model__2, y)

y = Z__10*a[1]-Z__11*b[1]+Z__20*a[2]-Z__21*b[2]+a[0]

(6)


Example.
I use as a generative process for Sx and Ycg (given the order you use in your question Ycg is a vector of
x values and Sx a vector of y values).
For instance

Seed := 175860914374622;
randomize(Seed);

175860914374622

 

175860914374622

(7)

GenerativeProcess := proc(N, sigma::positive)
  local R, params, values, model, f, X, Y:

  R      := () -> rand(-10..10)() / rand(1..20)():

  params := convert(indets(Model__0) minus {x, y}, list):
  values := [seq(R(), n=1..numelems(params))]:
  model  := eval(Model__0, params =~ values);

  print(model + Noemal(0, sigma)):

  f := unapply(rhs(model), x):
  X := Sample(Uniform(-1, 1), N)^+:
  Y := f~(X) + Sample(Normal(0, sigma), N)^+:
  return X, Y, params =~ values
end proc:
 

N     := 40:
sigma := 0.1:

Ycg, Sx, params := GenerativeProcess(N, sigma);

y+Noemal(0, .1) = (-(10/9)*x^2+(10/17)*x-3/10)/(-(2/3)*x^2+(7/11)*x+1)+Noemal(0, .1)

 

Vector[column](%id = 18446744078182577086), Vector[column](%id = 18446744078182577686), [a[0] = -3/10, a[1] = 10/17, a[2] = -10/9, b[1] = 7/11, b[2] = -2/3]

(8)


Construct new regressors

Model__lin;

z__10 := Ycg:
z__11 := Ycg *~Sx:
z__20 := Ycg^~2:
z__21 := Ycg^~2 *~Sx:
 

y = Z__10*a[1]-Z__11*b[1]+Z__20*a[2]-Z__21*b[2]+a[0]

(9)


Fit linear model

res := Fit(
        rhs(Model__lin)
        , `<|>`(z__10, z__11, z__20, z__21)
        , Sx
        , [Z__10, Z__11, Z__20, Z__21]
        , output=[leastsquaresfunction, parametervalues]
      )
  

[-HFloat(0.2777298343064421)+HFloat(0.5262609535304735)*Z__10-HFloat(1.1869703402160423)*Z__20-HFloat(0.6608454314903524)*Z__11+HFloat(0.6374123072191089)*Z__21, [a[0] = HFloat(-0.2777298343064421), a[1] = HFloat(0.5262609535304735), a[2] = HFloat(-1.1869703402160423), b[1] = HFloat(0.6608454314903524), b[2] = HFloat(-0.6374123072191089)]]

(10)


Go back to original model with regressor x

Z_fit := res[1]

-HFloat(0.2777298343064421)+HFloat(0.5262609535304735)*Z__10-HFloat(1.1869703402160423)*Z__20-HFloat(0.6608454314903524)*Z__11+HFloat(0.6374123072191089)*Z__21

(11)

aux  := indets(Model__lin) minus indets(Model__0):
Z2xy := solve(VarChange, aux);

X_fit := isolate(y = eval(Z_fit, Z2xy), y);

{Z__10 = x, Z__11 = x*y, Z__20 = x^2, Z__21 = x^2*y}

 

y = (-HFloat(0.2777298343064421)+HFloat(0.5262609535304735)*x-HFloat(1.1869703402160423)*x^2)/(1+HFloat(0.6608454314903524)*x-HFloat(0.6374123072191089)*x^2)

(12)


Compare observations to model predictions

plots:-display(
  ScatterPlot(Ycg, Sx)
  , plot(rhs(X_fit), x=(min..max)(Ycg))
)

 


Compare the values of the parameters for the generative process and their estimations

`<|>`(
  `<,>`("True values", "Fitted values")
  , convert([evalf(params), res[2]], Matrix)
)

Matrix(2, 6, {(1, 1) = "True values", (1, 2) = a[0] = -.3000000000, (1, 3) = a[1] = .5882352941, (1, 4) = a[2] = -1.111111111, (1, 5) = b[1] = .6363636364, (1, 6) = b[2] = -.6666666667, (2, 1) = "Fitted values", (2, 2) = a[0] = -.277729834306442, (2, 3) = a[1] = .526260953530474, (2, 4) = a[2] = -1.18697034021604, (2, 5) = b[1] = .660845431490352, (2, 6) = b[2] = -.637412307219109})

(13)

 

 

Download Illustrative_example.mw


Add-on:
You could find instructive this comparison to NonlinearFitIllustrative_example_bad_NonlinearFit.mw

Based on Units but focuses only (as an example) on dimension LENGTH ans SI units (so "are" and "hectare" are not considered here).

 

restart

kernelopts(version)

`Maple 2015.2, APPLE UNIVERSAL OSX, Dec 20 2015, Build ID 1097895`

(1)

with(Units):

randomize(175847324662552);

175847324662552

(2)


For lengths, areas and volumes only

Ulength := [seq(cat(u, m), u in [m, c, d, ``, da, h, k])];

[mm, cm, dm, m, dam, hm, km]

(3)

Quizz := proc(N)
  local R, `Q&A`, n, r, u, U, e, x, X:

  R := proc()
    local  n, e:
    n := rand(0..100)():
    e := rand(-6..6)():
    if e <= 0 then
      return evalf[2](n*10^e)
    else
      return n*10^e
    end if:
  end proc:

  `Q&A` := NULL:
  for n from 1 to N do
    r     := R():
    u, U  := op(combinat:-randperm(Ulength)[1..2]):
    e     := rand(1..3)():
    x     := r*u^e:
    X     := convert(r*Unit(u^e), 'units', U^e):
    
    printf("\n%a = .......... %a\n", x, U^e):
    `Q&A` := `Q&A`, Printf("\n%a = %a %a\n", x, evalf[4](op(1, X)), U^e):
  end do:
  print():
  return [`Q&A`]
end proc:

answers := Quizz(5):


.14e-2*dam^3 = .......... mm^3

2.*dm^3 = .......... dam^3

.73*m^3 = .......... cm^3

3.*km^2 = .......... dam^2

89.*km^3 = .......... hm^3

 

(4)

# Correction

map(evalf, eval(answers, Printf=printf))


.14e-2*dam^3 = .1400e10 mm^3

2.*dm^3 = .2000e-5 dam^3

.73*m^3 = .7300e6 cm^3

3.*km^2 = .3000e5 dam^2

89.*km^3 = .8900e5 hm^3

 

[]

(5)

 

 

Download Units_Quizz_2.mw



An evolution which also enables quizz about velocities Units_Quizz_3.mw

 

 

 

 

discont(V, beta)
                            {-1, 0}

The disontinuity you see at x= 8 .1017 is an artifact due to number of Digits not large enough:

restart; with(plots)

C := 100:

alpha := .1:

m := 0.5e-1:

n := 0.9e-1:

k := 0.1e-1:

theta := (alpha*m+k)/(n-m):

lambda := 1.1*theta:

s := C*((1-theta/lambda)/(1+beta))^(1/beta):

V := (lambda*(n-m)*s*(1-(s/C)^beta)-k*s)/alpha;

16.50000000*(0.909090909e-1/(1+beta))^(1/beta)*(1-((0.909090909e-1/(1+beta))^(1/beta))^beta)-10.0*(0.909090909e-1/(1+beta))^(1/beta)

(1)

discont(V, beta)

{-1, 0}

(2)

``

plot(V, beta = 1 .. 10^18, color = [red]);

 

Digits := 20:
plot(V, beta = 1 .. 10^18)

 

# Or better:

plot(eval(V, beta=10^u), u=0..18)

 
 

``

Download MaplePrimes_Sep_19_sand15.mw

restart;

with(plots):
 

q := 0.9:
M := 1.01:

points := []:

for tau from 0.01 by 0.01 to 0.09 do

    eq1 := M^2*(1 - sqrt(1 - 2*x/M^2))
           + f*tau*(1 - exp(x/tau))
           + (2*(1-f)/(3*q-1))*(1 - (1+(q-1)*x)^((3*q-1)/(2*q-2))) = 0:

    eq2 := 1/sqrt(1 - 2*x/M^2)
           - f*exp(x/tau)
           - (1-f)*(1+(q-1)*x)^((q+1)/(2*q-2)) = 0:


    
    sol := fsolve({eq1, eq2}, {f, x}, {f=0.01..0.9, x=0.001..2});

    if sol::function then
        printf("tau=%.2f  ->  no solution found\n", tau);
    else
        fval := eval(f, sol);
        xval := eval(x, sol);
        points := [op(points), [fval, xval]];
        printf("tau=%.2f  ->  f=%.6f, x=%.6f\n", tau, fval, xval);
    end if:

od:
 

tau=0.01  ->  no solution found

tau=0.02  ->  no solution found

tau=0.03  ->  no solution found

tau=0.04  ->  no solution found

tau=0.05  ->  no solution found

tau=0.06  ->  no solution found

tau=0.07  ->  no solution found

tau=0.08  ->  no solution found

tau=0.09  ->  no solution found

 

# Change search ranges for fsolve or simply ignore them

points := []:

for tau from 0.01 by 0.01 to 0.09 do

    eq1 := M^2*(1 - sqrt(1 - 2*x/M^2))
           + f*tau*(1 - exp(x/tau))
           + (2*(1-f)/(3*q-1))*(1 - (1+(q-1)*x)^((3*q-1)/(2*q-2))) = 0:

    eq2 := 1/sqrt(1 - 2*x/M^2)
           - f*exp(x/tau)
           - (1-f)*(1+(q-1)*x)^((q+1)/(2*q-2)) = 0:


    
    sol := fsolve({eq1, eq2}, {f, x});

    if sol::function then
        printf("tau=%.2f  ->  no solution found\n", tau);
    else
        fval := eval(f, sol);
        xval := eval(x, sol);
        points := [op(points), [fval, xval]];
        printf("tau=%.2f  ->  f=%.6f, x=%.6f\n", tau, fval, xval);
    end if:

od:

tau=0.01  ->  f=0.851356, x=0.000000

tau=0.02  ->  f=0.856683, x=0.000000

tau=0.03  ->  f=0.001070, x=0.032805

tau=0.04  ->  f=0.001678, x=0.059488

tau=0.05  ->  f=0.002478, x=0.088687

tau=0.06  ->  f=0.003503, x=0.121100

tau=0.07  ->  f=0.004804, x=0.158068

tau=0.08  ->  f=0.006453, x=0.202525

tau=0.09  ->  f=0.008590, x=0.264387

 

points

[[.8513555201, 0.6964830318e-12], [.8566830887, 0.3912817906e-13], [0.1069988964e-2, 0.3280452951e-1], [0.1678278783e-2, 0.5948758155e-1], [0.2477635868e-2, 0.8868702371e-1], [0.3503054692e-2, .1211004031], [0.4803501518e-2, .1580677255], [0.6453490291e-2, .2025245703], [0.8590449799e-2, .2643866253]]

(1)

# You can select a posterior the points which suit you:

MyPoints := select(t -> is(t[1] > 0 and t[2] < 0.2), points);
print():
MyPoints := select(t -> is(t[1] > 0 and t[2] < 0.2 and t[2] > 0.01), points);

[[.8513555201, 0.6964830318e-12], [.8566830887, 0.3912817906e-13], [0.1069988964e-2, 0.3280452951e-1], [0.1678278783e-2, 0.5948758155e-1], [0.2477635868e-2, 0.8868702371e-1], [0.3503054692e-2, .1211004031], [0.4803501518e-2, .1580677255]]

 

 

[[0.1069988964e-2, 0.3280452951e-1], [0.1678278783e-2, 0.5948758155e-1], [0.2477635868e-2, 0.8868702371e-1], [0.3503054692e-2, .1211004031], [0.4803501518e-2, .1580677255]]

(2)

pointplot(points, symbol=solidcircle, color=blue,
          labels=["f","x"]):

NULL

Download solve_sand15.mw

The type depends on how you defne the distribution.
If some attributes are those of a discrete random variable than Maple understands that this random variable is discrete.

restart

with(Statistics):


PDF should be resereved to continuous random variable

dist := Distribution(
          PDF = unapply( piecewise(x > 0, exp(-6)*6^x/x!, 0), x)
        ):

X := RandomVariable(dist):

exports(attributes(X)[3]);

PDF(X, x)

Conditions, PDF

 

piecewise(0 < x, exp(-6)*6^x/factorial(x), 0)

(1)


For a discrete random variable use ProbabilityFunction instead

dist := Distribution(
          ProbabilityFunction = unapply( piecewise(x > 0, exp(-6)*6^x/x!, 0), x)
        ):

X := RandomVariable(dist):

exports(attributes(X)[3]);

(attributes(X)[3]):-Type;
 

Conditions, ProbabilityFunction, Type

 

discrete

(2)


Look how Maple interpretes the probability function in terms of PDF

:-ProbabilityFunction('X') = ProbabilityFunction(X, x);

:-PDF('X') = PDF(X, x)

ProbabilityFunction(X) = piecewise(0 < x, exp(-6)*6^x/factorial(x), 0)

 

PDF(X) = sum(exp(-6)*6^k*Dirac(x-k)/factorial(k), k = 1 .. infinity)

(3)
 

 

Download discrete.mw

BY THE WAY: A correct definition of the ProbabilityFunction is

ProbabilityFunction = unapply( piecewise(x >= 0, exp(-6)*6^x/x!, 0), x)

Indeed

# With your definition 

sum(PDF(X, x), x=-infinity..+infinity) = 1-e-6

# With my correction above

sum(PDF(X, x), x=-infinity..+infinity) = 1

Once correctly defined, the probability function is that of a Poisson random variable with parameter equal to 6.
I don't know if you were aware of this and if your question is to be interpreted in general terms as “How do you define a discrete random variable?” (the illustration provided being useless and only adding to the confusion), or if you didn't realize this, in which case there is no need to define a user distribution because the Poisson distribution already exists.

If you are interested to know how to define a parameterized distribution (assuming, of course, that you don't want to use Poisson 😉 , something I could understand given some general drawbacks concerning the way Maple defines discrete random variable, see below).

Look to this simplified version of the worksheet parameterized_discrete.mw, where a more complete definition of Distribution is made and which provides more details including a comparison to Statistics:-RandomVariable(Poisson(...)) :

restart

with(Statistics):


A parameterized Distribution

dist := proc(a)
  Distribution(
    ProbabilityFunction = unapply( piecewise(x >= 0, exp(-a)*a^x/x!, 0), x)
  )
end proc:

X := RandomVariable(dist(alpha)):

:-ProbabilityFunction('X') = ProbabilityFunction(X, x);

ProbabilityFunction(X) = piecewise(0 <= x, exp(-alpha)*alpha^x/factorial(x), 0)

(1)

sum(ProbabilityFunction(X, x), x=-infinity..+infinity)

1

(2)

X := RandomVariable(dist(6)):

:-ProbabilityFunction('X') = ProbabilityFunction(X, x);

ProbabilityFunction(X) = piecewise(0 <= x, exp(-6)*6^x/factorial(x), 0)

(3)
 

 

Funny:  you can use PDF (even if it is ugly) and lure Maple this way: discrete.mw

@AHSAN 

restart

with(plots):
with(plottools):
with(Statistics):

Z := Sample(Normal(0, 1), 10):

A := Array([seq(2+sin(10*Pi*i/15+2*Z[i]), i = 1..5)]):

B := Array([seq(2+sin(10*Pi*i/15+3*Z[i]), i = 1..5)]):

C := Array([seq(2+sin(10*Pi*i/15+4*Z[i]), i = 1..5)]):

p0 := ColumnGraph([A, B, C], legend = ["A", "B", "C"]):

NbGroups := 5:
NbConds  := 3:

dy := 0.5:

p1 := display(
  transform((x, y, z) -> [x, z, y])(extrude(p0, 0..dy))

  , seq(PLOT3D(CURVES([[0, dy, z], [NbGroups-1+NbConds/(NbConds+1), dy, z]], COLOR(RGB, 0.8$3))), z = 0.25..2.5, 0.25)
  , seq(PLOT3D(CURVES([[0, 0, z], [0, dy, z]], COLOR(RGB, 0.8$3))), z = 0.25..2.5, 0.25)
  , seq(textplot3d([i+(NbConds/2)/(NbConds+1), 0, 0, convert(i+1, string)], align=below, font=[Helvetica, bold, 15]), i=0..NbGroups-1)

  , scaling=constrained
  , orientation=[-66, 76, 2]
  , axis[1]=[tickmarks=[]]
  , axis[2]=[tickmarks=[]]
  , labels=["", "", "H"]
):

p1;

 

CondNames := ["A", "B", "C"]:

if NbConds::odd then
  dx := [seq(NbGroups-1+(NbConds/2)/(NbConds+1) + 1/(NbConds+1)*i, i=-(NbConds-1)/2..(NbConds-1)/2)]
else
  dx := [seq(NbGroups-1+(NbConds/2)/(NbConds+1) + 1/(NbConds+1)*i, i=-NbConds/2..NbConds/2)]
end if:

p2 := display(
  p1,
  seq(textplot3d([dx[i], 0, 0.25, CondNames[i]], align=below, font=[Helvetica, bold, 15], color=white), i=1..NbConds)
):

p2;

 

Download Example_3D.mw

lhs(Eq1) > 0 in the range you provide.

Beware_of_the_ranges.mw


The first problem you face is that LaTeX is a sotware system for typesetting documents and not a computer algebra system.
Maple  syntax relies upon maths while LaTeX syntax is most permissive: this means that if any Maple formula can be translated into a syntactically correct LaTeX expression the converse is not true.
For instance this LaTeX expression (where the "/" sign means, in the context of the paper, "conditionnally to") has no Maple meanings:


Another example taken at random from Goodle image.

So, assuming the LaTeX expression has a Maple sense, which would be the case if the author did some Maple work, converted some extression into LaTeX and imported it into its LaTeX development tool (Overleaf for instance), one may think that the LaTeX formula can be translated in Maple.


The LaTeX case

I assume below that you have the LaTeX source not only its compiled rendering.

The a priori best way:
(1) As you have Maple 2025 you must dispose of the MathML:-FromLaTeX function, see HERE.
As I am not that lucky I can't test it to tell you about its possible limitations.


(2) There exists some online translators (AI Chat for instance) but you must extremely very careful on the translation they provide (specifically if your expressions involve vectors, matrices, derivatives, integral).
To be confident in the online translator you could be tempted to use procced this way

  1. In Maple:
    1. write some expression expr
    2. Convert it to LaTeX.
    3. Copy the output.
  2. In the online translator:
    1. Paste the LaTeX expression.
    2. Ask for its Maple translation
    3. Copy it
  3. In Maple:
    1. Paste the experssion the translator provided
    2. compare to the original maple expression expr.

 

(3)  An old (2016) Mapleprimes question (note that SnuggleTeX no longer exists).

Nevertheless, even if you can use FromLaTeX, point (2) might be usefull depending on FromLaTeX limitations (if any).


The "no LaTeX source code" case

Without LaTeX source (pdf article for instance) a copy-paste of the expression into a Maple worksheet it possible but it is very likely that you will have to correct it.
For instance the illustration I provided at the top of this reply, copied-pasted in maple document mode gives this





The "hand-written" case

Is obviously even more complex than this latter because a first step to convert a handwritten note into LaTeX seems mandatory (see for instance NotesToLaTeX, MathWrite, or others)



In My Opinion,

Even i this work can seem tedious to you, I believe that converting handwritten notes or LaTeX compiled outputs (not a source code) into Maple expressions is much safer if you do it carefully without looking forthe help of some automatic translator.



Last point:

If you are interesting in writting something like LHS = RHS1 = RHS2 = ....RHSN in Maple, just browse the Mapleprimes site: there has been a recent question about whose @acer replied.

You should upload your worksheet so that we can verify there is no spurious character in the expression you want to expand.

restart

kernelopts(version)

`Maple 2015.2, APPLE UNIVERSAL OSX, Dec 20 2015, Build ID 1097895`

(1)

f := (m-1/2)^4 - (5/2+12*c)*(m-1/2)^2 + 1/4*(9/4+12*c)

(m-1/2)^4-(5/2+12*c)*(m-1/2)^2+9/16+3*c

(2)

fe := expand(f);

m^4-12*c*m^2-2*m^3+12*c*m-m^2+2*m

(3)

# check:

simplify(f-fe);

0

(4)

# In detail:

g := (m-1/2)^4 - (5/2+12*c)*(m-1/2)^2 + ``(1/4)*(9/4+12*c);
print():


opg  := [op(g)]:
eopg := map(expand, opg):
print~(`&equiv;`~(opg, eopg)):

(m-1/2)^4-(5/2+12*c)*(m-1/2)^2+``(1/4)*(9/4+12*c)

 

 

`&equiv;`((m-1/2)^4, m^4-2*m^3+(3/2)*m^2-(1/2)*m+1/16)

 

`&equiv;`(-(5/2+12*c)*(m-1/2)^2, -(5/2)*m^2+(5/2)*m-5/8-12*c*m^2+12*c*m-3*c)

 

`&equiv;`(``(1/4)*(9/4+12*c), 9/16+3*c)

(5)

`&equiv;`~(g, add(eopg));

`&equiv;`((m-1/2)^4-(5/2+12*c)*(m-1/2)^2+``(1/4)*(9/4+12*c), m^4-12*c*m^2-2*m^3+12*c*m-m^2+2*m)

(6)
 

 

Download expand.mw

The solution has a single branch if v__0 <> 0 and two branches instead.

branches.mw

I dont know if the recent Maple versions (I use Maple 2015) have something equivalent to solve/parametric for dsolve?
Something that could write

dsolve({ode, ic}, {v(x)}, 'parametric'='full', 'parameters'={v__0});

and would deliver all the branches depending on the value of v__0

You define f[i] as a function of two arguments p2 and pr but you use only one in the following command:

Cs := {seq(beta[i] = subs(pr = 0, f[i](pr)), i = 1 .. 2)};

Error, invalid input: f[1] uses a 2nd argument, pr, which is missing


Correct yourself the definition of Cs or f[i]'s before going further.

Have you ever considered that attaching an article might give the impression that you are asking the community to code and solve the equations it contains? Perhaps it is not surprising that no one then replies to you..

This put aside the error message is very clear: if a bvp problem is set over a compact subset of IR, let's say the interval [-1, +1], the the boundaries conditions have to be imposed at y=-1 and y=+1 and not aa some third point inside [-1, +1].
So your problem IS NOT a bvp problem.

This said there are essentially two ways to tackle it

  1. Solve independently two bvp problems:
    1. One on interval [-1, 0] with unfixed values at y=0
    2. One on interval [0, 1] with still unfixed calues at y=0... but the same of the previous problem
    3.  Connect the left and right solution in order that your "internal conditions at y=0" are fullfilled.
    4. To do this you will have to solve some minimization ptoblem.
       
  2. Solve a single bvp problem (the optionI choosed) over [-1, +1] with some bc unfixed and iterate over the solutions until your "internal conditions at y=0" are fullfilled.

The attached file explains you how to do this.

Last point: in OdeSys you have not given the parameter R some numeric value. So I arbitrarily choosed to take R=0.
If this choice does not suit yoy feel free to take another value.

The attached worksheet does not browse all the values of j as, for questions of time, I only considered j=1.
Here again it's up to you to complete the job if you agree with it.

symetry_paper_work_sand15.mw

1 2 3 4 5 6 7 Page 1 of 7