mmcdara

6279 Reputation

17 Badges

7 years, 325 days

MaplePrimes Activity


These are answers submitted by mmcdara

Use intat~ instead of int 

restart

kernelopts(version)

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

(1)

with(ScientificConstants):
with(Units):

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

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

(2)

int(g, t*Unit('s'));

Error, (in int) integration range or variable must be specified in the second argument, got t*Units:-Unit(s)

 

# Workaround

intat(g, tau=t*Unit('s'));
simplify(%, 'Unit');

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

 

9.80665*t*Units:-Unit(('m')/('s'))

(3)

# The vector case

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

a := `<,>`(0, 1, 0) *~ g

a := Vector(3, {(1) = 0., (2) = 9.80665*Units:-Unit(m/s^2), (3) = 0.})

(4)

v[0] := `<,>`(v__0, 0, 0) *~Unit('m'/'s')

v[0] := Vector(3, {(1) = `#msub(mi("v"),mi("0"))`*Unit('m'/'s'), (2) = 0, (3) = 0})

(5)

V := simplify(intat~(a, tau=t*Unit('s')));

dv := Vector(3, {(1) = 0., (2) = 9.80665*t*Units:-Unit(m/s), (3) = 0.})

(6)

Download Int_vector_Unit.mw

If you want to integrate twice to get the position X do this (remember that v__0 and t are dimensionless velocity and time) :

X := intat~(intat~(a, tau=t*Unit('s')), tau=t*Unit('s'))
     + 
     intat~(v[0], tau=t*Unit('s')):

simplify(X);

    v__0*t*Unit('m')*e[x] + 9.80665*t^2*Unit('m')*e[y] + 0.*e[z]

I was capable to derive a surprisingly simple expression based on the hypergeometric function but not the one you want (or more precisely I got a LaguerreL based expression but the result is not correct).
In the route to the proof I have used two properties of the HermiteH polynomials (physical form) that Maple "doesn't know" and I did some transformations by hand (for instance commuting the operators Int and Sum. I don't know if this is something that can be done using directly Maple?

Maybe someone will  reuse the material presented in the attached file and finalize the demonstration?
Goog luck.

restart;

with(IntegrationTools):

L := (m, n) -> exp(a^2)*a^(abs(n-m))*sqrt(2^(abs(n-m))*min(n, m)!/max(n, m)!)*LaguerreL(min(n, m), abs(n-m), -2*a^2);

proc (m, n) options operator, arrow; exp(a^2)*a^abs(n-m)*sqrt(2^abs(n-m)*factorial(min(n, m))/factorial(max(n, m)))*LaguerreL(min(n, m), abs(n-m), -2*a^2) end proc

(1)

H := 2^(-m/2-n/2)*exp(lambda^2*sigma^2/4)/sqrt(n!)/sqrt(m!*Pi)*Int(HermiteH(m,s+lambda*sigma/2)*HermiteH(n,s+lambda*sigma/2)*exp(-s^2), s=-infinity..infinity);


# Set a = lambda*sigma/2

H := unapply(eval(H, lambda=a*2/sigma), (m, n));
 

2^(-(1/2)*m-(1/2)*n)*exp((1/4)*lambda^2*sigma^2)*(Int(HermiteH(m, s+(1/2)*lambda*sigma)*HermiteH(n, s+(1/2)*lambda*sigma)*exp(-s^2), s = -infinity .. infinity))/(factorial(n)^(1/2)*(factorial(m)*Pi)^(1/2))

 

proc (m, n) options operator, arrow; 2^(-(1/2)*m-(1/2)*n)*exp(a^2)*(Int(HermiteH(m, s+a)*HermiteH(n, s+a)*exp(-s^2), s = -infinity .. infinity))/(factorial(n)^(1/2)*(factorial(m)*Pi)^(1/2)) end proc

(2)

# I think this is the key property to use
#
# Basic Property 1:

BP1m := HermiteH(m, s+a) = Sum(binomial(m, i)*HermiteH(i, s)*(2*a)^(m-i), i=0..m);
BP1n := HermiteH(n, s+a) = Sum(binomial(n, j)*HermiteH(j, s)*(2*a)^(n-j), j=0..n);

HermiteH(m, s+a) = Sum(binomial(m, i)*HermiteH(i, s)*(2*a)^(m-i), i = 0 .. m)

 

HermiteH(n, s+a) = Sum(binomial(n, j)*HermiteH(j, s)*(2*a)^(n-j), j = 0 .. n)

(3)

# So:
SBP1 := eval(H(m, n), [BP1m, BP1n]);

2^(-(1/2)*m-(1/2)*n)*exp(a^2)*(Int((Sum(binomial(m, i)*HermiteH(i, s)*(2*a)^(m-i), i = 0 .. m))*(Sum(binomial(n, j)*HermiteH(j, s)*(2*a)^(n-j), j = 0 .. n))*exp(-s^2), s = -infinity .. infinity))/(factorial(n)^(1/2)*(factorial(m)*Pi)^(1/2))

(4)

#  Commute Sum operators and the Int operator

OutInt   := mul(remove(has, [op(SBP1)], Int)[]);
GetI     := GetIntegrand(SBP1):
SumTerms := op([1, 1], GetI)*op([2, 1], GetI):

S, R := selectremove(has, [op(%)], HermiteH):

SBP1 := OutInt * Sum(Sum(mul(R[]) * Int(mul(S[])*exp(-s^2), s=-infinity..+infinity), op([1, 2], GetI)), op([2, 2], GetI));

2^(-(1/2)*m-(1/2)*n)*exp(a^2)/(factorial(n)^(1/2)*(factorial(m)*Pi)^(1/2))

 

2^(-(1/2)*m-(1/2)*n)*exp(a^2)*(Sum(Sum(binomial(m, i)*(2*a)^(m-i)*binomial(n, j)*(2*a)^(n-j)*(Int(HermiteH(i, s)*HermiteH(j, s)*exp(-s^2), s = -infinity .. infinity)), i = 0 .. m), j = 0 .. n))/(factorial(n)^(1/2)*(factorial(m)*Pi)^(1/2))

(5)

# If we expand SSum we will get (m+1)*(n+1) terms of the form
# T[i, j] = int(C[i, j]*HermiteH(i, x)*Hermite(j, x)*exp(-x^2), x=-infinity..+infinity).
#
# But we know that
#     1)  T[i, j<>i] = 0
#     2)  T[i, j=i] = sqrt(Pi)*2^i*i!
#
# It is easy to see that all terms of the form T[i, i] are got for i = 0.. min(m, n)
# Then the product of the two dpuble sums reduces to

# Basic Property 2:

BP2 := Int(HermiteH(i, s)*HermiteH(j, s)*exp(-s^2), s = -infinity .. infinity) = piecewise(j=i, sqrt(Pi)*2^i*i!, 0);

BP2 := Int(HermiteH(i, s)*HermiteH(j, s)*exp(-s^2), s = -infinity .. infinity) = piecewise(j = i, sqrt(Pi)*2^i*factorial(i), 0)

(6)

# Use basic property 2 into SBP1

SBP2 := eval(SBP1, BP2);

# Symplify by hand

SBP2 := OutInt * Sum(binomial(m, i)*binomial(n, i)*(2*a)^(m-i)*(2*a)^(n-i)*sqrt(Pi)*2^i*i!, i=0..min(m, n));

SBP2 := 2^(-(1/2)*m-(1/2)*n)*exp(a^2)*(Sum(Sum(binomial(m, i)*(2*a)^(m-i)*binomial(n, j)*(2*a)^(n-j)*piecewise(j = i, sqrt(Pi)*2^i*factorial(i), 0), i = 0 .. m), j = 0 .. n))/(sqrt(factorial(n))*sqrt(factorial(m)*Pi))

 

2^(-(1/2)*m-(1/2)*n)*exp(a^2)*(Sum(binomial(m, i)*binomial(n, i)*(2*a)^(m-i)*(2*a)^(n-i)*Pi^(1/2)*2^i*factorial(i), i = 0 .. min(m, n)))/(factorial(n)^(1/2)*(factorial(m)*Pi)^(1/2))

(7)

# Ask a little help (a huge one in fact) to Maple:

eval(SBP2, Sum=sum);

2^(-(1/2)*m-(1/2)*n)*exp(a^2)*Pi^(1/2)*((2*a)^(m+n)*hypergeom([-m, -n], [], (1/2)/a^2)-2*binomial(m, min(m, n)+1)*binomial(n, min(m, n)+1)*(2*a)^(m-2-2*min(m, n)+n)*factorial(min(m, n)+1)*hypergeom([1, -m+1+min(m, n), -n+1+min(m, n)], [2+min(m, n)], (1/2)/a^2)*2^min(m, n))/(factorial(n)^(1/2)*(factorial(m)*Pi)^(1/2))

(8)

# Two a priori cases... which appear to give the same result

mLEn := eval(SBP2, Sum=sum) assuming m <= n;
mGTn := eval(SBP2, Sum=sum) assuming m > n;

simplify(mLEn - mGTn);

2^(-(1/2)*m-(1/2)*n)*exp(a^2)*Pi^(1/2)*(2*a)^(m+n)*hypergeom([-n, -m], [], (1/2)/a^2)/(factorial(n)^(1/2)*(factorial(m)*Pi)^(1/2))

 

2^(-(1/2)*m-(1/2)*n)*exp(a^2)*Pi^(1/2)*(2*a)^(m+n)*hypergeom([-n, -m], [], (1/2)/a^2)/(factorial(n)^(1/2)*(factorial(m)*Pi)^(1/2))

 

0

(9)

# Define the expression of H(m, n) this way

Hermite2Hypergeom := mLEn

2^(-(1/2)*m-(1/2)*n)*exp(a^2)*Pi^(1/2)*(2*a)^(m+n)*hypergeom([-m, -n], [], (1/2)/a^2)/(factorial(n)^(1/2)*(factorial(m)*Pi)^(1/2))

(10)

# Check case m < n

simplify(L(1, 3), 'LaguerreL');
simplify(eval(H(1, 3), Int=int));
simplify(eval(Hermite2Hypergeom, [m=1, n=3]))

(1/3)*exp(a^2)*a^2*6^(1/2)*(2*a^2+3)

 

(1/3)*exp(a^2)*a^2*6^(1/2)*(2*a^2+3)

 

(1/3)*exp(a^2)*a^2*6^(1/2)*(2*a^2+3)

(11)

# Check case m = n

normal(simplify(L(3, 3), 'LaguerreL'));
simplify(eval(H(3, 3), Int=int));
simplify(eval(Hermite2Hypergeom, [m=3, n=3]))

(1/3)*exp(a^2)*(4*a^6+18*a^4+18*a^2+3)

 

(1/3)*exp(a^2)*(4*a^6+18*a^4+18*a^2+3)

 

(1/3)*exp(a^2)*(4*a^6+18*a^4+18*a^2+3)

(12)

# Check case m > n

simplify(L(3, 1), 'LaguerreL');
simplify(eval(H(3, 1), Int=int));
simplify(eval(Hermite2Hypergeom, [m=3, n=1]))

(1/3)*exp(a^2)*a^2*6^(1/2)*(2*a^2+3)

 

(1/3)*exp(a^2)*a^2*6^(1/2)*(2*a^2+3)

 

(1/3)*exp(a^2)*a^2*6^(1/2)*(2*a^2+3)

(13)

# Now let's start from Hermite2Hypergeom and try to convert it as LaguerreL.
#
# Roughly:

convert(Hermite2Hypergeom, `1F1`) assuming m::posint, n::posint:
convert~(%, LaguerreL):
Hermite2Laguerre := convert~(%, factorial);

2^(-(1/2)*m-(1/2)*n)*exp(a^2)*Pi^(1/2)*(2*a)^(m+n)*factorial(-n+m-1)*LaguerreL(m, n-m, -2*a^2)*factorial(m)*factorial(n-m)/(factorial(n)^(3/2)*(factorial(m)*Pi)^(1/2)*factorial(-n-1)*(-2*a^2)^m)

(14)

# That seems pretty, but there is a main problem: (-n-1)! = 0 for any positive integer n.
#
# (Note that the above expression does not account for the sign of m-n, which seems already
# quite suspect... but Hermite2Hypergeom did not neither andtheresults where correct.)

eval(Hermite2Laguerre, m=3);
eval(%, n=1);

-(1/8)*2^(-3/2-(1/2)*n)*exp(a^2)*6^(1/2)*(2*a)^(3+n)*factorial(-n+2)*LaguerreL(3, n-3, -2*a^2)*factorial(n-3)/(factorial(n)^(3/2)*factorial(-n-1)*a^6)

 

Error, numeric exception: division by zero

 

# I have spent a lot of time to try and understand how to get rid of this problem without any success.
#
# So I feel that the formula L(m, n) is not that far but I am not capable to get it.
# I hope someone will and that and my work will help.

Download Q1_mmcdara.mw


Writting "In the above integral, r and sigma are the random variables" is incorrect: R and Sigma are txo random variables whose realizations noted r and sigma.
Assuming that this is what you wanted to say, here is the solution
 

restart:

with(Statistics):

R    := RandomVariable(Normal(k*sigma, sigma)):
f__R := unapply(PDF(R, r), r);

proc (r) options operator, arrow; (1/2)*2^(1/2)*exp(-(1/2)*(-k*sigma+r)^2/sigma^2)/(Pi^(1/2)*sigma) end proc

(1)

Sigma     := RandomVariable(LogNormal(m__Sigma, s__Sigma)):
f__Sigma := unapply( (PDF(Sigma, sigma) assuming sigma > 0), sigma);

proc (sigma) options operator, arrow; (1/2)*2^(1/2)*exp(-(1/2)*(ln(sigma)-m__Sigma)^2/s__Sigma^2)/(sigma*s__Sigma*Pi^(1/2)) end proc

(2)

# Compute the inner integral (it could be done head-on)


InnerIntegral := Int(r*f__R(r), r=-infinity..+infinity);
InnerIntegral := Mean(R);

Int((1/2)*r*2^(1/2)*exp(-(1/2)*(-k*sigma+r)^2/sigma^2)/(Pi^(1/2)*sigma), r = -infinity .. infinity)

 

k*sigma

(3)

OuterIntegral := upperbound -> Int(sigma*f__Sigma(sigma)*InnerIntegral, sigma=0..upperbound);
value(OuterIntegral(+infinity));
eval(%) assuming s__Sigma > 0;

OuterIntegral := proc (upperbound) options operator, arrow; Int(sigma*`#msub(mi("f"),mi("&Sigma;",fontstyle = "normal"))`(sigma)*InnerIntegral, sigma = 0 .. upperbound) end proc

 

(1/2)*2^(1/2)*k*piecewise(csgn(1/s__Sigma^2) = 1, exp(2*s__Sigma^2+2*m__Sigma)*2^(1/2)*csgn(1/s__Sigma)*s__Sigma*Pi^(1/2), infinity)/(s__Sigma*Pi^(1/2))

 

k*exp(2*s__Sigma^2+2*m__Sigma)

(4)

# A simpler way to get the previous result:

Moment(Sigma, 2)*coeff(InnerIntegral, sigma):
simplify(%);


# (here the call to Statistics:-Moment automatically accounts for the conditions
# on the parameters of Sigma:)

print():
(attributes(Sigma)[3]):-Conditions;

k*exp(2*s__Sigma^2+2*m__Sigma)

 

 

[m__Sigma::real, 0 < s__Sigma]

(5)

# What happens if the integration is done for upperbound U < +infinity?

value(OuterIntegral(U)):
solution := simplify( (eval(%) assuming s__Sigma > 0, U > 0) )

(1/2)*k*exp(2*s__Sigma^2+2*m__Sigma)*(1+erf((1/2)*2^(1/2)*(-2*s__Sigma^2+ln(U)-m__Sigma)/s__Sigma))

(6)

# Now give to these parameters the values you want

indets(solution, name) minus {Pi}

{U, k, m__Sigma, s__Sigma}

(7)

 

 

Download integral.mw

Now, why is this semantic point important?

  • If r and sigma were INDEED  random variables, then the inner integral (and even the outer) would make  no sense. Even the simpler int(r, r=-infinity..+infinity) makes no sense.
     
  • If Sigma (I prefer using a capital do refer to a random variable) was a random variable, then
    R, defined by PDF(R, r) = Normal(k*sigma, sigma) where sigma is a realization of Sigma, would not be a random variable but a random process.
    A (1D) random proces is a (finite or not) collection of random variables indexed by a number (which can be itself the realization of a random variable).
    Thus, "R conditionned by Sigma=sigma", which I am going to note R[sigma] is a random variable (a gaussian one with mean k*sigma and standard deviation sigma), but not R which is an infinite collection of random variables R[sigma] with sigma=0..+infinity.

This kind of issues are widely adressed (and answered here).
To acoid them be EXTREMELY careful when using 2D input mode in document style ... or avoid it (my advice) and use 1D input mode and worksheet instead.

The explanation is in bold brown courier font in the attached file:
(I changed the construction of U, V and Upsilon for my Maple 2015 doesn't support yours... but this is anecdotal).

NULL

restart

U := proc (x, t) options operator, arrow; U0*exp(I*omega*t)*cos(m*Pi*x/L) end proc; -1; W := proc (x, t) options operator, arrow; W0*exp(I*omega*t)*sin(m*Pi*x/L) end proc; -1; Upsilon := proc (x, t) options operator, arrow; Upsilon0*exp(I*omega*t)*cos(m*Pi*x/L) end proc

proc (x, t) options operator, arrow; Upsilon0*exp(I*omega*t)*cos(m*Pi*x/L) end proc

(1)

m := 1; -1; nu := .38; -1; h := 10*l; -1; l := 17.6*10^(-6); -1; mu := E/(2*(1+nu)); -1; E := 1.44*10^9; -1; rho := 1220; -1; L := 10*h; -1; R := 10^10

10000000000

(2)
 

Phi := proc (z) options operator, arrow; z*(1-(4/3)*z^2/h^2) end proc

proc (z) options operator, arrow; z*(1-(4/3)*z^2/h^2) end proc

(3)

Q__11 := E/(-nu^2+1); 1; Q__55 := E/(2*(1+nu)); 1; G := E/(2*(1+nu)); 1; Q[110] := E*(int(z^0, z = -(1/2)*h .. (1/2)*h))/(-nu^2+1); 1; Q[111] := int(E*z/(-nu^2+1), z = -(1/2)*h .. (1/2)*h); 1; Q[112] := int(E*z^2/(-nu^2+1), z = -(1/2)*h .. (1/2)*h); 1; Q[550] := int(E*z^0/(2*(1+nu)), z = -(1/2)*h .. (1/2)*h); 1; Q[551] := int(E*z/(2*(1+nu)), z = -(1/2)*h .. (1/2)*h); 1; Q[552] := int(E*z^2/(2*(1+nu)), z = -(1/2)*h .. (1/2)*h)

1683029453.

 

521739130.4

 

521739130.4

 

296213.1837

 

0.

 

0.7646249649e-3

 

91826.08695

 

0.

 

0.2370337391e-3

(4)
 

Q[113] := int(Q[11]*Phi(z), z = -(1/2)*h .. (1/2)*h); 1; Q[114] := int(Q[11]*Phi(z)^2, z = -(1/2)*h .. (1/2)*h); 1; Q[115] := int(Q[11]*z*Phi(z), z = -(1/2)*h .. (1/2)*h)

0.

 

0.2942228318e-12*Q[11]

 

0.3634517333e-12*Q[11]

(5)

Q[553] := int(Q[55]*Phi(z), z = -(1/2)*h .. (1/2)*h); 1; Q[554] := int(Q[55]*Phi(z)^2, z = -(1/2)*h .. (1/2)*h); 1; Q[555] := int(Q[55]*z*Phi(z), z = -(1/2)*h .. (1/2)*h)

0.

 

0.2942228318e-12*Q[55]

 

0.3634517333e-12*Q[55]

(6)

Q[556] := int(Q[55]*(diff(Phi(z), z)), z = -(1/2)*h .. (1/2)*h); 1; Q[557] := int(Q[55]*z*(diff(Phi(z), z)), z = -(1/2)*h .. (1/2)*h); 1; Q[558] := int(Q[55]*Phi(z)*(diff(Phi(z), z)), z = -(1/2)*h .. (1/2)*h); 1; Q[559] := int(Q[55]*(diff(Phi(z), z))^2, z = -(1/2)*h .. (1/2)*h)

0.1173333333e-3*Q[55]

 

0.

 

0.

 

0.9386666667e-4*Q[55]

(7)

# WHAT DID YOU DEFINE?
#
# To see really what the variables you defined are do this

lprint(anames(user));

# Observe there are 3 "Q's":
#
#   on one side Q__11, Q__55
#   on the other one Q
#
# This means that what appears as Q subscriped 110, 550, 113, 553, 556, ... are the same variable Q
# with indices 110, 550, 113, 553, 556, ....
#
# To see this:

indices(Q);

# Q is a table:

eval(Q);

G, h, l, Phi, rho, Q__55, W, Upsilon, U, E, m, nu, Q, L, R, mu, Q__11

 

[111], [550], [110], [551], [557], [556], [559], [558], [553], [552], [555], [554], [114], [115], [112], [113]

 

table( [( 111 ) = 0., ( 550 ) = 91826.08695, ( 110 ) = 296213.1837, ( 551 ) = 0., ( 557 ) = 0., ( 556 ) = 0.1173333333e-3*Q[55], ( 559 ) = 0.9386666667e-4*Q[55], ( 558 ) = 0., ( 553 ) = 0., ( 552 ) = 0.2370337391e-3, ( 555 ) = 0.3634517333e-12*Q[55], ( 554 ) = 0.2942228318e-12*Q[55], ( 114 ) = 0.2942228318e-12*Q[11], ( 115 ) = 0.3634517333e-12*Q[11], ( 112 ) = 0.7646249649e-3, ( 113 ) = 0. ] )

(8)

# Note that neither [11] nor [55] indices of Q

beta[1] := int(mu*z^0, z = -(1/2)*h .. (1/2)*h); 1; beta[2] := int(mu*z, z = -(1/2)*h .. (1/2)*h); 1; beta[3] := int(mu*z^2, z = -(1/2)*h .. (1/2)*h)

91826.08697

 

0.

 

0.2370337392e-3

(9)

beta[4] := int(mu*Phi(z), z = -(1/2)*h .. (1/2)*h); 1; beta[5] := int(mu*Phi(z)^2, z = -(1/2)*h .. (1/2)*h); 1; beta[6] := int(mu*z*Phi(z), z = -(1/2)*h .. (1/2)*h); 1; beta[7] := int(mu*(diff(Phi(z), z)), z = -(1/2)*h .. (1/2)*h)

0.

 

0.1535075644e-3

 

0.1896269913e-3

 

61217.39131

(10)

beta[8] := int(mu*z*(diff(Phi(z), z)), z = -(1/2)*h .. (1/2)*h);

0.

 

0.

 

48973.91305

(11)

beta[11] := int(mu*(diff(Phi(z), z, z)), z = -(1/2)*h .. (1/2)*h);

0.

 

0.

 

0.1581027668e14

(12)

N__T := 0; -1; H__x := 0; -1; N__0 := 0; -1; K__2 := 0; -1; K__1 := 0; -1; DD := 0

0

(13)

I__0 := int(rho, z = -(1/2)*h .. (1/2)*h); 1; I__1 := int(rho*z, z = -(1/2)*h .. (1/2)*h); 1; I__2 := int(rho*z^2, z = -(1/2)*h .. (1/2)*h); 1; I__3 := int(rho*Phi(z), z = -(1/2)*h .. (1/2)*h); 1; I__4 := int(rho*Phi(z)^2, z = -(1/2)*h .. (1/2)*h); 1; I__5 := int(rho*z*Phi(z), z = -(1/2)*h .. (1/2)*h)

.2147200000

 

0.

 

0.5542638933e-9

 

0.

 

0.3589518547e-9

 

0.4434111147e-9

(14)

NULL

eq1 := -Q__110*(diff(U(x, t), x, x))+Q__111*(diff(W(x, t), x, x, x))-Q__113*(diff(Upsilon(x, t), x, x))-Q__110*(diff(W(x, t), x))/R-`#mi("\`Q__556\`")`*Upsilon(x, t)/R+`#mrow(mi("\`Q__550\`"),mo("&sdot;"),mi("U"),mfenced(mrow(mi("x"),mo("&comma;"),mi("t"))))`/R^2-`#mi("\`Q__551\`")`*(diff(W(x, t), x))/R^2+`#mrow(mi("\`Q__553\`"),mo("&sdot;"),mi("&Upsi;",fontstyle = "normal"),mfenced(mrow(mi("x"),mo("&comma;"),mi("t"))))`/R^2+l^2*beta__1*(diff(W(x, t), x, x, x))/(4*R)-beta__7*l^2*(diff(Upsilon(x, t), x, x))/(8*R)-beta__1*l^2*(diff(U(x, t), x, x))/(8*R^2)+l^2*beta__2*(diff(W(x, t), x, x, x))/(8*R^2)-beta__4*l^2*(diff(Upsilon(x, t), x, x))/(8*R^2)+I__0*(diff(U(x, t), t, t))-I__1*(diff(W(x, t), x, t, t))+I__3*(diff(Upsilon(x, t), t, t))

3186210.099*Q__110*U0*exp(I*omega*t)*cos(1784.995826*x)-5687371727.*Q__111*W0*exp(I*omega*t)*cos(1784.995826*x)+3186210.099*Q__113*Upsilon0*exp(I*omega*t)*cos(1784.995826*x)-0.1784995826e-6*Q__110*W0*exp(I*omega*t)*cos(1784.995826*x)-(1/10000000000)*`#msub(mi("\`Q"),mi("556\`"))`*Upsilon0*exp(I*omega*t)*cos(1784.995826*x)+(1/100000000000000000000)*`#mrow(msub(mi("\`Q"),mi("550\`")),mo("&sdot;"),mi("U"),mfenced(mrow(mi("x"),mo("&comma;"),mi("t"))))`-0.1784995826e-16*`#msub(mi("\`Q"),mi("551\`"))`*W0*exp(I*omega*t)*cos(1784.995826*x)+(1/100000000000000000000)*`#mrow(msub(mi("\`Q"),mi("553\`")),mo("&sdot;"),mi("&Upsi;",fontstyle = "normal"),mfenced(mrow(mi("x"),mo("&comma;"),mi("t"))))`-0.4404300665e-10*beta__1*W0*exp(I*omega*t)*cos(1784.995826*x)+0.1233700550e-13*beta__7*Upsilon0*exp(I*omega*t)*cos(1784.995826*x)+0.1233700550e-23*beta__1*U0*exp(I*omega*t)*cos(1784.995826*x)-0.2202150332e-20*beta__2*W0*exp(I*omega*t)*cos(1784.995826*x)+0.1233700550e-23*beta__4*Upsilon0*exp(I*omega*t)*cos(1784.995826*x)-.2147200000*U0*omega^2*exp(I*omega*t)*cos(1784.995826*x)

(15)

# WHAT ARE THE INDETERMINATES OF eq1?

lprint~(indets(eq1, name));

Q__110

Q__111
Q__113
U0
W0
beta__1
beta__2
beta__4
beta__7
omega
t
x
`#mi("\`Q__551\`")`
`#mi("\`Q__556\`")`
Upsilon0
`#mrow(mi("\`Q__550\`"),mo("&sdot;"),mi("U"),mfenced(mrow(mi("x"),mo("&comma;"),mi("t"))))`
`#mrow(mi("\`Q__553\`"),mo("&sdot;"),mi("&Upsi;",fontstyle = "normal"),mfenced(mrow(mi("x"),mo("&comma;"),mi("t"))))`

 

{}

(16)

# Compare to the table Q whose indices are remembered here

{indices(Q, nolist)}

{110, 111, 112, 113, 114, 115, 550, 551, 552, 553, 554, 555, 556, 557, 558, 559}

(17)

# So the answer to your question is obvious:
#
# You defined Q[110] and use Q__110 in eq1 (same thing for the other quantities)
#
# And the remedy is simple: be more carefull when you se 2D mode
# And a good advice (IMO), let's drop it and use 1D input mode instead to avoid this kind of issues.

Download TRASH_mmcdara.mw

Note that your problem does not limit tho the Qs but extends to the betas.

The problem with the 2D inputs is that what you see is not what you get.
Here you falsely believed that the variables in eq1 where those you previously defined.

You provide initial conditions for g[s] (equation (3)), one of them involving D(g[s]](1), but you solve an ODE which does not contain any derivative of g[s].

If g[s] must respect a differential eqution, write it and solve a 2-by-2 differential system instead of a single differential equation.

If g[s] is not governed by a differential equation, specify the expression of g[s](x) and remove all references to g[s] in the boundary conditions.
Here is an illustration of the "Known g[s] case":  Known_gs_mmcdara.mw

this proposal.mw ?

restart

 # Parameters always positive

assume(0 < sigma__d, 0 < sigma__v):
interface(showassumed=0):

Diff('lambda__1', sigma__v) = sqrt(5)/(5*sigma__d):
f__1 := rhs(%):
Diff('lambda__1', sigma__d) = -sqrt(5)*sigma__v/(5*sigma__d^2):
f__2 := abs(rhs(%)):

use plots in
  display(
    inequal(
      f__1-f__2 > 0
      , sigma__v=0.001..10, sigma__d=0.001..10
      , color="Chartreuse"
      , optionsexcluded=[color="Niagara DarkOrchid"]
    )
    ,
    # For a more complex case than yours
    # (here a simple plot(solve(f__1=f__2, sigma__d), sigma__d=0.001..10)
    # would be enough)
    implicitplot(
      f__1=f__2
      , sigma__v=0.001..10, sigma__d=0.001..10
      , color=blue
      , thickness=5
      , grid=[100, 100]
      , legend=typeset('f__1' = 'f__2')
    )
    , textplot([3, 7, typeset('f__1' > 'f__2')], font=[Times, bold, 14])
    , textplot([7, 3, typeset('f__1' < 'f__2')], font=[Times, bold, 14], color=white)
    , labels=[typeset(sigma__v), typeset(sigma__d)]
  )
end use;

 

Download proposal.mw

(customize it as you want) 

See also add_on.mw

@jalal 


To get this histogram you present:

  • Given the nature of your data I suppose they are already consist in a serie of elements of the type [C[k], N[k]] where
    • C[k] is the range of the kth class,
    • N[k] is the number of elements C[k] contains.

This is what I call synthetic data by opposition to raw data (see below).
 

  • If it is so you have two possibilities:
    1. To build an Statistics:-Histogram you must have raw data, not synthetic data.
      So either you have these raw data or you don't.
      • If you have them, let R the list/vector of incomes, just use Statistics:-Histogram with the ad hoc customization (first plot in the attached file).
         
      • If you do not have them, you can recreate fake raw data from synthetic data so that their histogram is the one you want (second plot in the attached file).
         
    2. If you only have synthetic data and bo not necessarily want to use Statistics:-Histogram, then you can proceed as explaine in the third plot in the attached file.



 

restart


CASE 1: YOU HAVE RAW DATA


Simulation of hypothetic raw data

R := Statistics:-Sample(Triangular(6, 97, 53), 10^3):


Raw data histogram

C := [seq(0..100, 10)]:
Statistics:-Histogram(
  R
  , frequencyscale=absolute
  , binbounds=C
  , axis[1]=[tickmarks=C]
  , view=[-1..101, default]
  , color="Pink"
  , transparency=0.5
  , labels=["Salaires (milliers de \$)", "Nombre d'employés"]
  , labeldirections=[default, vertical]
  , labelfont=[Tahoma, 10]
  , size=[700, 400]
)

 


CASE 2: YOU DO NOT HAVE RAW DATA BUT SYNTHETIC DATA ONLY

... AND YOU WANT TO USE  Statistics:-Histogram


Here are the synthetic data which correspond to the raw data R above

Classes       := [seq(C[k]..C[k+1], k=1..numelems(C)-1)]:
SyntheticData := Statistics:-TallyInto(R, Classes):
SyntheticData := Matrix(numelems(Classes), 2, (i, j) -> `if`(j=1, lhs(SyntheticData[i]), rhs(SyntheticData[i])));
 

SyntheticData := Matrix(10, 2, {(1, 1) = HFloat(0.0) .. HFloat(10.0), (1, 2) = 3, (2, 1) = HFloat(10.0) .. HFloat(20.0), (2, 2) = 43, (3, 1) = HFloat(20.0) .. HFloat(30.0), (3, 2) = 92, (4, 1) = HFloat(30.0) .. HFloat(40.0), (4, 2) = 141, (5, 1) = HFloat(40.0) .. HFloat(50.0), (5, 2) = 195, (6, 1) = HFloat(50.0) .. HFloat(60.0), (6, 2) = 218, (7, 1) = HFloat(60.0) .. HFloat(70.0), (7, 2) = 152, (8, 1) = HFloat(70.0) .. HFloat(80.0), (8, 2) = 95, (9, 1) = HFloat(80.0) .. HFloat(90.0), (9, 2) = 51, (10, 1) = HFloat(90.0) .. HFloat(100.0), (10, 2) = 10})

(1)



If you want to use Statistics:-Histogram you must create hypothetic raw data in a correct way.
The reconstruction below produces raw data that have the same histogram as above.

S := [seq((add(op(SyntheticData[k, 1]))/2)$(SyntheticData[k, 2]), k=1..numelems(Classes))]:

Statistics:-Histogram(
  S
  , frequencyscale=absolute
  , binbounds=C
  , axis[1]=[tickmarks=C]
  , view=[-1..101, default]
  , color="Pink"
  , transparency=0.5
  , labels=["Salaires (milliers de \$)", "Nombre d'employés"]
  , labeldirections=[default, vertical]
  , labelfont=[Tahoma, 10]
  , size=[700, 400]
)

 


CASE 3: YOU DO NOT HAVE RAW DATA BUT SYNTHETIC DATA ONLY

... AND YOU DO NOT NECESSARILLY WANT TO USE  Statistics:-Histogram

bars := [
          seq(
            plottools:-rectangle(
              [op(1, SyntheticData[k, 1]), 0]
              , [op(2, SyntheticData[k, 1]), SyntheticData[k, 2]]
              , color="Pink"
              , transparency=0.5
            )
            , k=1..numelems(Classes)
          )
        ]:

plots:-display(
  bars
  , labels=["Salaires (milliers de \$)", "Nombre d'employés"]
  , labeldirections=[default, vertical]
  , labelfont=[Tahoma, 10]
  , size=[700, 400]
)

 

 


 

Download Different_options.mw

It's likely that some typo or hidden character appeared in ode2 (bold red text).
As I never use  the document mode with 2D inputs, the attached file displays two types of fonts: so correct your own file the way I did if this annoys you.

restart

kernelopts(version)

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

(1)

with(PDETools):

with(LinearAlgebra):

with(PolynomialIdeals):

with(plots):

declare(f(x), prime = x);

f(x)*`will now be displayed as`*f

 

`derivatives with respect to`*x*`of functions of one variable will now be displayed with '`

 

theta(x)*`will now be displayed as`*theta

 

`derivatives with respect to`*x*`of functions of one variable will now be displayed with '`

(2)

NULL

ode1 := diff(f(x), x, x, x)+(1/2*(1+p*(diff(u[e](x), x))/u[e]))*f(x)*(diff(f(x), x, x))+p*(D(u[e]))(1-(diff(f(x), x))^2)/u[e]-p*((diff(f(x), x))*(diff(f(x), x, x))-(diff(f(x), x))*(diff(f(x), x, x))) = 0;

diff(diff(diff(f(x), x), x), x)+(1/2)*(1+p*(diff(u[e](x), x))/u[e])*f(x)*(diff(diff(f(x), x), x))+p*(D(u[e]))(1-(diff(f(x), x))^2)/u[e] = 0

 

(diff(diff(theta(x), x), x))/Pr+(1/2)*(1+p*(diff(u[e](x), x))/u[e])*f(x)*(diff(theta(x), x))+16*p*a*sigma*L*T[e]^3*(1-theta(x))/(rho*c[p]*u[e]*U[infinity])+U[infinity]^2*p*(diff(f(x), x))*u[e]*(diff(u[e](x), x))/(c[p]*(T[w]-T[e]))-U[infinity]^2*u[e]^2*(diff(diff(f(x), x), x))^2/(c[p]*(T[w]-T[e])) = 0

(3)

 

bcs := f(0) = 0, (D(f))(0) = 0, (D(f))(8) = 1;

f(0) = 0, (D(f))(0) = 0, (D(f))(8) = 1

 

theta(0) = 0, theta(8) = 1

 

U[infinity]*(1-p/L)

(4)

a1 := [a = 1/2, L = 5, U[infinity] = 1, Pr = .7, c[p] = 1004, T[w] = 290, T[e] = 300, abs(T[w]-T[e]) = 10, sigma = 5.67*10^(-8), rho = 1, p = .5];

[a = 1/2, L = 5, U[infinity] = 1, Pr = .7, c[p] = 1004, T[w] = 290, T[e] = 300, abs(-T[w]+T[e]) = 10, sigma = 0.5670000000e-7, rho = 1, p = .5]

 

[a = 1/2, L = 6, U[infinity] = 1, Pr = .7, c[p] = 1004, T[w] = 290, T[e] = 300, abs(-T[w]+T[e]) = 10, sigma = 0.5670000000e-7, rho = 1, p = .5]

 

[a = 1/2, L = 8, U[infinity] = 1, Pr = .7, c[p] = 1004, T[w] = 290, T[e] = 300, abs(-T[w]+T[e]) = 10, sigma = 0.5670000000e-7, rho = 1, p = .5]

(5)

b1 := eval(ode2, a1);
b2 := eval(ode2, a2);
b3 := eval(ode2, a3);

1.428571429*(diff(diff(theta(x), x), x))+.5000000000*f(x)*(diff(theta(x), x))+0.3388446214e-1-0.3388446214e-1*theta(x)+0.8067729084e-4*(diff(diff(f(x), x), x))^2 = 0

 

1.428571429*(diff(diff(theta(x), x), x))+.5000000000*f(x)*(diff(theta(x), x))+0.3992205722e-1-0.3992205722e-1*theta(x)+0.8369300576e-4*(diff(diff(f(x), x), x))^2 = 0

 

1.428571429*(diff(diff(theta(x), x), x))+.5000000000*f(x)*(diff(theta(x), x))+0.5204653387e-1-0.5204653387e-1*theta(x)+0.8754046315e-4*(diff(diff(f(x), x), x))^2 = 0

(6)

c1 := dsolve({bcs1, eval(b1, f(x) = x)}, numeric); c1(0);
c2 := dsolve({bcs1, eval(b2, f(x) = x)}, numeric); c2(0);
c3 := dsolve({bcs1, eval(b3, f(x) = x)}, numeric); c3(0);

proc (x_bvp) local res, data, solnproc, _ndsol, outpoint, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then outpoint := evalf[_EnvDSNumericSaveDigits](x_bvp) else outpoint := evalf(x_bvp) end if; data := Array(1..4, {(1) = proc (outpoint) local X, Y, YP, yout, errproc, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; X := Vector(12, {(1) = .0, (2) = .7110816593827365, (3) = 1.4252612944250316, (4) = 2.1450075553641894, (5) = 2.870698607458807, (6) = 3.6016439190865284, (7) = 4.335950433427024, (8) = 5.072381437739985, (9) = 5.809614980858862, (10) = 6.547193597737968, (11) = 7.281757576380688, (12) = 8.0}, datatype = float[8], order = C_order); Y := Matrix(12, 2, {(1, 1) = .0, (1, 2) = .4938484453872593, (2, 1) = .33592640765199755, (2, 2) = .43888481949800257, (3, 1) = .6120099871592314, (3, 2) = .32852292968068963, (4, 1) = .8036666351004551, (4, 2) = .20577920071677938, (5, 1) = .9150410854740033, (5, 2) = .10721924293183477, (6, 1) = .9688640080289379, (6, 2) = 0.4626167885323632e-1, (7, 1) = .990390459943389, (7, 2) = 0.165006856163208e-1, (8, 1) = .9975131515965905, (8, 2) = 0.4861478032432699e-2, (9, 1) = .999462252457013, (9, 2) = 0.11840633668385558e-2, (10, 1) = .9999041167369538, (10, 2) = 0.23851519861008805e-3, (11, 1) = .9999870438515986, (11, 2) = 0.4007921134972896e-4, (12, 1) = 1.0, (12, 2) = 0.58484376487594675e-5}, datatype = float[8], order = C_order); YP := Matrix(12, 2, {(1, 1) = .4938484453872593, (1, 2) = -0.237191234908843e-1, (2, 1) = .43888481949800257, (2, 2) = -.12498027451545507, (3, 1) = .32852292968068963, (3, 2) = -.1730836385803869, (4, 1) = .20577920071677938, (4, 2) = -.1591461343770899, (5, 1) = .10721924293183477, (5, 2) = -.10974309693499718, (6, 1) = 0.4626167885323632e-1, (6, 2) = -0.59054851436050264e-1, (7, 1) = 0.165006856163208e-1, (7, 2) = -0.25269084092256245e-1, (8, 1) = 0.4861478032432699e-2, (8, 2) = -0.8689730687886674e-2, (9, 1) = 0.11840633668385558e-2, (9, 2) = -0.24203881956517473e-2, (10, 1) = 0.23851519861008805e-3, (10, 2) = -0.5488360802490405e-3, (11, 1) = 0.4007921134972896e-4, (11, 2) = -0.10245379376869548e-3, (12, 1) = 0.58484376487594675e-5, (12, 2) = -0.16375625411613887e-4}, datatype = float[8], order = C_order); errproc := proc (x_bvp) local outpoint, X, Y, yout, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; Digits := 15; outpoint := evalf(x_bvp); X := Vector(12, {(1) = .0, (2) = .7110816593827365, (3) = 1.4252612944250316, (4) = 2.1450075553641894, (5) = 2.870698607458807, (6) = 3.6016439190865284, (7) = 4.335950433427024, (8) = 5.072381437739985, (9) = 5.809614980858862, (10) = 6.547193597737968, (11) = 7.281757576380688, (12) = 8.0}, datatype = float[8], order = C_order); Y := Matrix(12, 2, {(1, 1) = .0, (1, 2) = -0.45502255784570956e-7, (2, 1) = 0.9444480543660673e-8, (2, 2) = -0.32701855724739926e-7, (3, 1) = -0.13111776668557004e-6, (3, 2) = 0.7051671015686677e-7, (4, 1) = -0.10085891808028474e-6, (4, 2) = 0.12193667249947745e-6, (5, 1) = 0.13308574074702087e-6, (5, 2) = -0.7660319092840244e-7, (6, 1) = 0.16436590818059128e-6, (6, 2) = -0.18881505625378492e-6, (7, 1) = -0.20116971988762455e-7, (7, 2) = 0.7685426832740135e-9, (8, 1) = -0.9569009446204095e-7, (8, 2) = 0.13118291896195288e-6, (9, 1) = -0.31100640223178834e-7, (9, 2) = 0.4639384647499807e-7, (10, 1) = 0.1584209066120383e-7, (10, 2) = -0.35974084912581454e-7, (11, 1) = 0.11642475203736646e-7, (11, 2) = -0.2764855504318049e-7, (12, 1) = .0, (12, 2) = -0.1925248521887889e-8}, datatype = float[8], order = C_order); if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "right" then return X[12] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(1.8881505625378492e-7) elif outpoint = "errorproc" then error "this is already the error procedure" elif outpoint = "rawdata" then return [2, 12, [theta(x), diff(theta(x), x)], X, Y] else return ('procname')(x_bvp) end if end if; if outpoint < X[1] or X[12] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[12] end if; V := array([1 = 4, 2 = 0]); if Digits <= trunc(evalhf(Digits)) then L := Vector(4, 'datatype' = 'float'[8]); yout := Vector(2, 'datatype' = 'float'[8]); evalhf(`dsolve/numeric/lagrange`(12, 2, X, Y, outpoint, var(yout), var(L), var(V))) else L := Vector(4, 'datatype' = 'sfloat'); yout := Vector(2, 'datatype' = 'sfloat'); `dsolve/numeric/lagrange`(12, 2, X, Y, outpoint, yout, L, V) end if; [x = outpoint, seq('[theta(x), diff(theta(x), x)]'[i] = yout[i], i = 1 .. 2)] end proc; if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "method" then return "bvp" elif outpoint = "right" then return X[12] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(1.8881505625378492e-7) elif outpoint = "errorproc" then return eval(errproc) elif outpoint = "rawdata" then return [2, 12, "depnames", X, Y, YP] else error "non-numeric value" end if end if; if outpoint < X[1] or X[12] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[12] end if; if Digits <= trunc(evalhf(Digits)) and (_EnvInFsolve <> true or _EnvDSNumericSaveDigits <= trunc(evalhf(Digits))) then V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false)  ] ); L := Matrix(7, 2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0, (3, 1) = .0, (3, 2) = .0, (4, 1) = .0, (4, 2) = .0, (5, 1) = .0, (5, 2) = .0, (6, 1) = .0, (6, 2) = .0, (7, 1) = .0, (7, 2) = .0}, datatype = float[8], order = C_order); yout := Vector(2, {(1) = .0, (2) = .0}, datatype = float[8]); evalhf(`dsolve/numeric/hermite`(12, 2, X, Y, YP, outpoint, var(yout), var(L), var(V))) else if _EnvInFsolve = true then Digits := _EnvDSNumericSaveDigits end if; V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false)  ] ); L := Matrix(7, 2, {(1, 1) = 0., (1, 2) = 0., (2, 1) = 0., (2, 2) = 0., (3, 1) = 0., (3, 2) = 0., (4, 1) = 0., (4, 2) = 0., (5, 1) = 0., (5, 2) = 0., (6, 1) = 0., (6, 2) = 0., (7, 1) = 0., (7, 2) = 0.}, order = C_order); yout := Vector(2, {(1) = 0., (2) = 0.}); `dsolve/numeric/hermite`(12, 2, X, Y, YP, outpoint, yout, L, V) end if; [outpoint, seq(yout[i], i = 1 .. 2)] end proc, (2) = Array(0..0, {}), (3) = [x, theta(x), diff(theta(x), x)], (4) = 0}); solnproc := data[1]; if not type(outpoint, 'numeric') then if outpoint = "solnprocedure" then return eval(solnproc) elif member(outpoint, ["start", "left", "right", "errorproc", "rawdata", "order", "error"]) then return solnproc(x_bvp) elif outpoint = "sysvars" then return data[3] elif procname <> unknown then return ('procname')(x_bvp) else _ndsol := pointto(data[2][0]); return ('_ndsol')(x_bvp) end if end if; try res := solnproc(outpoint); [x = res[1], seq('[theta(x), diff(theta(x), x)]'[i] = res[i+1], i = 1 .. 2)] catch: error  end try end proc

 

[x = 0., theta(x) = HFloat(0.0), diff(theta(x), x) = HFloat(0.4938484453872593)]

 

proc (x_bvp) local res, data, solnproc, _ndsol, outpoint, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then outpoint := evalf[_EnvDSNumericSaveDigits](x_bvp) else outpoint := evalf(x_bvp) end if; data := Array(1..4, {(1) = proc (outpoint) local X, Y, YP, yout, errproc, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; X := Vector(12, {(1) = .0, (2) = .7109194006521007, (3) = 1.4249868973749966, (4) = 2.1446957738695653, (5) = 2.8704181759615643, (6) = 3.6014311331509616, (7) = 4.335815327649273, (8) = 5.07232139321572, (9) = 5.809624985093265, (10) = 6.54727068912088, (11) = 7.281870021183384, (12) = 8.0}, datatype = float[8], order = C_order); Y := Matrix(12, 2, {(1, 1) = .0, (1, 2) = .49766174438535715, (2, 1) = .3375757993634647, (2, 2) = .44007354626176926, (3, 1) = .6138389834751503, (3, 2) = .32815856063407745, (4, 1) = .8049921933517767, (4, 2) = .2049224933509404, (5, 1) = .9157759787447678, (5, 2) = .1064979573437486, (6, 1) = .9691886486674092, (6, 2) = 0.45847467449322866e-1, (7, 1) = .990506554830634, (7, 2) = 0.16320440331941466e-1, (8, 1) = .9975470109277185, (8, 2) = 0.479979457345354e-2, (9, 1) = .9994703351484253, (9, 2) = 0.1167147574584417e-2, (10, 1) = .9999056808640605, (10, 2) = 0.23475927713373983e-3, (11, 1) = .9999872681820644, (11, 2) = 0.3939868164434801e-4, (12, 1) = 1.0, (12, 2) = 0.5746972025352357e-5}, datatype = float[8], order = C_order); YP := Matrix(12, 2, {(1, 1) = .49766174438535715, (1, 2) = -0.279454400456164e-1, (2, 1) = .44007354626176926, (2, 2) = -.12801162336374522, (3, 1) = .32815856063407745, (3, 2) = -.1744590166938988, (4, 1) = .2049224933509404, (4, 2) = -.15927332083416956, (5, 1) = .1064979573437486, (5, 2) = -.10934646266608739, (6, 1) = 0.45847467449322866e-1, (6, 2) = -0.58651810580892944e-1, (7, 1) = 0.16320440331941466e-1, (7, 2) = -0.25032143866202835e-1, (8, 1) = 0.479979457345354e-2, (8, 2) = -0.8589685100784845e-2, (9, 1) = 0.1167147574584417e-2, (9, 2) = -0.23880431153508012e-2, (10, 1) = 0.23475927713373983e-3, (10, 2) = -0.5405971765590906e-3, (11, 1) = 0.3939868164434801e-4, (11, 2) = -0.10076942378371007e-3, (12, 1) = 0.5746972025352357e-5, (12, 2) = -0.16091521666157516e-4}, datatype = float[8], order = C_order); errproc := proc (x_bvp) local outpoint, X, Y, yout, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; Digits := 15; outpoint := evalf(x_bvp); X := Vector(12, {(1) = .0, (2) = .7109194006521007, (3) = 1.4249868973749966, (4) = 2.1446957738695653, (5) = 2.8704181759615643, (6) = 3.6014311331509616, (7) = 4.335815327649273, (8) = 5.07232139321572, (9) = 5.809624985093265, (10) = 6.54727068912088, (11) = 7.281870021183384, (12) = 8.0}, datatype = float[8], order = C_order); Y := Matrix(12, 2, {(1, 1) = .0, (1, 2) = -0.4573949665850558e-7, (2, 1) = 0.6242874960547203e-8, (2, 2) = -0.3057749989248847e-7, (3, 1) = -0.12966272475504071e-6, (3, 2) = 0.7217465704512478e-7, (4, 1) = -0.940390491554525e-7, (4, 2) = 0.1178112154625901e-6, (5, 1) = 0.1357038240627767e-6, (5, 2) = -0.8040712939043494e-7, (6, 1) = 0.1607319279832156e-6, (6, 2) = -0.18568759788316905e-6, (7, 1) = -0.2282862282805592e-7, (7, 2) = 0.4786547870156748e-8, (8, 1) = -0.94843869899127e-7, (8, 2) = 0.13055152558459717e-6, (9, 1) = -0.2970037963750006e-7, (9, 2) = 0.4434338125388376e-7, (10, 1) = 0.16200082159941674e-7, (10, 2) = -0.36390277880281064e-7, (11, 1) = 0.11561105767115588e-7, (11, 2) = -0.27217683490459804e-7, (12, 1) = .0, (12, 2) = -0.16987712182963502e-8}, datatype = float[8], order = C_order); if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "right" then return X[12] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(1.8568759788316905e-7) elif outpoint = "errorproc" then error "this is already the error procedure" elif outpoint = "rawdata" then return [2, 12, [theta(x), diff(theta(x), x)], X, Y] else return ('procname')(x_bvp) end if end if; if outpoint < X[1] or X[12] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[12] end if; V := array([1 = 4, 2 = 0]); if Digits <= trunc(evalhf(Digits)) then L := Vector(4, 'datatype' = 'float'[8]); yout := Vector(2, 'datatype' = 'float'[8]); evalhf(`dsolve/numeric/lagrange`(12, 2, X, Y, outpoint, var(yout), var(L), var(V))) else L := Vector(4, 'datatype' = 'sfloat'); yout := Vector(2, 'datatype' = 'sfloat'); `dsolve/numeric/lagrange`(12, 2, X, Y, outpoint, yout, L, V) end if; [x = outpoint, seq('[theta(x), diff(theta(x), x)]'[i] = yout[i], i = 1 .. 2)] end proc; if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "method" then return "bvp" elif outpoint = "right" then return X[12] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(1.8568759788316905e-7) elif outpoint = "errorproc" then return eval(errproc) elif outpoint = "rawdata" then return [2, 12, "depnames", X, Y, YP] else error "non-numeric value" end if end if; if outpoint < X[1] or X[12] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[12] end if; if Digits <= trunc(evalhf(Digits)) and (_EnvInFsolve <> true or _EnvDSNumericSaveDigits <= trunc(evalhf(Digits))) then V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false)  ] ); L := Matrix(7, 2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0, (3, 1) = .0, (3, 2) = .0, (4, 1) = .0, (4, 2) = .0, (5, 1) = .0, (5, 2) = .0, (6, 1) = .0, (6, 2) = .0, (7, 1) = .0, (7, 2) = .0}, datatype = float[8], order = C_order); yout := Vector(2, {(1) = .0, (2) = .0}, datatype = float[8]); evalhf(`dsolve/numeric/hermite`(12, 2, X, Y, YP, outpoint, var(yout), var(L), var(V))) else if _EnvInFsolve = true then Digits := _EnvDSNumericSaveDigits end if; V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false)  ] ); L := Matrix(7, 2, {(1, 1) = 0., (1, 2) = 0., (2, 1) = 0., (2, 2) = 0., (3, 1) = 0., (3, 2) = 0., (4, 1) = 0., (4, 2) = 0., (5, 1) = 0., (5, 2) = 0., (6, 1) = 0., (6, 2) = 0., (7, 1) = 0., (7, 2) = 0.}, order = C_order); yout := Vector(2, {(1) = 0., (2) = 0.}); `dsolve/numeric/hermite`(12, 2, X, Y, YP, outpoint, yout, L, V) end if; [outpoint, seq(yout[i], i = 1 .. 2)] end proc, (2) = Array(0..0, {}), (3) = [x, theta(x), diff(theta(x), x)], (4) = 0}); solnproc := data[1]; if not type(outpoint, 'numeric') then if outpoint = "solnprocedure" then return eval(solnproc) elif member(outpoint, ["start", "left", "right", "errorproc", "rawdata", "order", "error"]) then return solnproc(x_bvp) elif outpoint = "sysvars" then return data[3] elif procname <> unknown then return ('procname')(x_bvp) else _ndsol := pointto(data[2][0]); return ('_ndsol')(x_bvp) end if end if; try res := solnproc(outpoint); [x = res[1], seq('[theta(x), diff(theta(x), x)]'[i] = res[i+1], i = 1 .. 2)] catch: error  end try end proc

 

[x = 0., theta(x) = HFloat(0.0), diff(theta(x), x) = HFloat(0.49766174438535715)]

 

proc (x_bvp) local res, data, solnproc, _ndsol, outpoint, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then outpoint := evalf[_EnvDSNumericSaveDigits](x_bvp) else outpoint := evalf(x_bvp) end if; data := Array(1..4, {(1) = proc (outpoint) local X, Y, YP, yout, errproc, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; X := Vector(12, {(1) = .0, (2) = .7105918488760237, (3) = 1.4244328611291581, (4) = 2.1440659423124293, (5) = 2.8698512980998347, (6) = 3.6010004098896413, (7) = 4.335541208202285, (8) = 5.072198845420454, (9) = 5.809644194733494, (10) = 6.5474258600340365, (11) = 7.282096976039094, (12) = 8.0}, datatype = float[8], order = C_order); Y := Matrix(12, 2, {(1, 1) = .0, (1, 2) = .5052552537006175, (2, 1) = .3408462755204489, (2, 2) = .442422215142914, (3, 1) = .6174560224034431, (3, 2) = .32742681953867314, (4, 1) = .8076074937875736, (4, 2) = .2032242438976165, (5, 1) = .9172227818187347, (5, 2) = .10507373726470326, (6, 1) = .9698264569016594, (6, 2) = 0.45031941250935854e-1, (7, 1) = .9907341917436213, (7, 2) = 0.15966437963481905e-1, (8, 1) = .9976132778134335, (8, 2) = 0.4678915550247515e-2, (9, 1) = .9994861265216441, (9, 2) = 0.11340645752511988e-2, (10, 1) = .9999087319608427, (10, 2) = 0.22742703929550063e-3, (11, 1) = .9999877052438856, (11, 2) = 0.38072228239622816e-4, (12, 1) = 1.0, (12, 2) = 0.5549311882673039e-5}, datatype = float[8], order = C_order); YP := Matrix(12, 2, {(1, 1) = .5052552537006175, (1, 2) = -0.364325736980702e-1, (2, 1) = .442422215142914, (2, 2) = -.13404823355722884, (3, 1) = .32742681953867314, (3, 2) = -.17717619408563023, (4, 1) = .2032242438976165, (4, 2) = -.15951351711336537, (5, 1) = .10507373726470326, (5, 2) = -.10855688752010623, (6, 1) = 0.45031941250935854e-1, (6, 2) = -0.5785531343159658e-1, (7, 1) = 0.15966437963481905e-1, (7, 2) = -0.2456567964351268e-1, (8, 1) = 0.4678915550247515e-2, (8, 2) = -0.8393290947591843e-2, (9, 1) = 0.11340645752511988e-2, (9, 2) = -0.2324700819301251e-2, (10, 1) = 0.22742703929550063e-3, (10, 2) = -0.5244967168305387e-3, (11, 1) = 0.38072228239622816e-4, (11, 2) = -0.9748390992631722e-4, (12, 1) = 0.5549311882673039e-5, (12, 2) = -0.15538073266825603e-4}, datatype = float[8], order = C_order); errproc := proc (x_bvp) local outpoint, X, Y, yout, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; Digits := 15; outpoint := evalf(x_bvp); X := Vector(12, {(1) = .0, (2) = .7105918488760237, (3) = 1.4244328611291581, (4) = 2.1440659423124293, (5) = 2.8698512980998347, (6) = 3.6010004098896413, (7) = 4.335541208202285, (8) = 5.072198845420454, (9) = 5.809644194733494, (10) = 6.5474258600340365, (11) = 7.282096976039094, (12) = 8.0}, datatype = float[8], order = C_order); Y := Matrix(12, 2, {(1, 1) = .0, (1, 2) = -0.4621216507042107e-7, (2, 1) = 0.18871763805575788e-9, (2, 2) = -0.26440976380442328e-7, (3, 1) = -0.12642870795085007e-6, (3, 2) = 0.7507398457807075e-7, (4, 1) = -0.8067853301949186e-7, (4, 2) = 0.10949148171767816e-6, (5, 1) = 0.1404802397309745e-6, (5, 2) = -0.8760075826073426e-7, (6, 1) = 0.15342379606531546e-6, (6, 2) = -0.1792563287988756e-6, (7, 1) = -0.28042403695503607e-7, (7, 2) = 0.12597322677227419e-7, (8, 1) = -0.9307097617997284e-7, (8, 2) = 0.12914164433307568e-6, (9, 1) = -0.26946194672378708e-7, (9, 2) = 0.4029811073645486e-7, (10, 1) = 0.1687782559135978e-7, (10, 2) = -0.37161040827639946e-7, (11, 1) = 0.11392020809816244e-7, (11, 2) = -0.2635259175010777e-7, (12, 1) = .0, (12, 2) = -0.12553181464882845e-8}, datatype = float[8], order = C_order); if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "right" then return X[12] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(1.792563287988756e-7) elif outpoint = "errorproc" then error "this is already the error procedure" elif outpoint = "rawdata" then return [2, 12, [theta(x), diff(theta(x), x)], X, Y] else return ('procname')(x_bvp) end if end if; if outpoint < X[1] or X[12] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[12] end if; V := array([1 = 4, 2 = 0]); if Digits <= trunc(evalhf(Digits)) then L := Vector(4, 'datatype' = 'float'[8]); yout := Vector(2, 'datatype' = 'float'[8]); evalhf(`dsolve/numeric/lagrange`(12, 2, X, Y, outpoint, var(yout), var(L), var(V))) else L := Vector(4, 'datatype' = 'sfloat'); yout := Vector(2, 'datatype' = 'sfloat'); `dsolve/numeric/lagrange`(12, 2, X, Y, outpoint, yout, L, V) end if; [x = outpoint, seq('[theta(x), diff(theta(x), x)]'[i] = yout[i], i = 1 .. 2)] end proc; if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "method" then return "bvp" elif outpoint = "right" then return X[12] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(1.792563287988756e-7) elif outpoint = "errorproc" then return eval(errproc) elif outpoint = "rawdata" then return [2, 12, "depnames", X, Y, YP] else error "non-numeric value" end if end if; if outpoint < X[1] or X[12] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[12] end if; if Digits <= trunc(evalhf(Digits)) and (_EnvInFsolve <> true or _EnvDSNumericSaveDigits <= trunc(evalhf(Digits))) then V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false)  ] ); L := Matrix(7, 2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0, (3, 1) = .0, (3, 2) = .0, (4, 1) = .0, (4, 2) = .0, (5, 1) = .0, (5, 2) = .0, (6, 1) = .0, (6, 2) = .0, (7, 1) = .0, (7, 2) = .0}, datatype = float[8], order = C_order); yout := Vector(2, {(1) = .0, (2) = .0}, datatype = float[8]); evalhf(`dsolve/numeric/hermite`(12, 2, X, Y, YP, outpoint, var(yout), var(L), var(V))) else if _EnvInFsolve = true then Digits := _EnvDSNumericSaveDigits end if; V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false)  ] ); L := Matrix(7, 2, {(1, 1) = 0., (1, 2) = 0., (2, 1) = 0., (2, 2) = 0., (3, 1) = 0., (3, 2) = 0., (4, 1) = 0., (4, 2) = 0., (5, 1) = 0., (5, 2) = 0., (6, 1) = 0., (6, 2) = 0., (7, 1) = 0., (7, 2) = 0.}, order = C_order); yout := Vector(2, {(1) = 0., (2) = 0.}); `dsolve/numeric/hermite`(12, 2, X, Y, YP, outpoint, yout, L, V) end if; [outpoint, seq(yout[i], i = 1 .. 2)] end proc, (2) = Array(0..0, {}), (3) = [x, theta(x), diff(theta(x), x)], (4) = 0}); solnproc := data[1]; if not type(outpoint, 'numeric') then if outpoint = "solnprocedure" then return eval(solnproc) elif member(outpoint, ["start", "left", "right", "errorproc", "rawdata", "order", "error"]) then return solnproc(x_bvp) elif outpoint = "sysvars" then return data[3] elif procname <> unknown then return ('procname')(x_bvp) else _ndsol := pointto(data[2][0]); return ('_ndsol')(x_bvp) end if end if; try res := solnproc(outpoint); [x = res[1], seq('[theta(x), diff(theta(x), x)]'[i] = res[i+1], i = 1 .. 2)] catch: error  end try end proc

 

[x = 0., theta(x) = HFloat(0.0), diff(theta(x), x) = HFloat(0.5052552537006175)]

(7)

:

p1 := odeplot(c1, [x, theta(x)], 0 .. 8, colour = green):

display(p1, p2, p3);

 

NULL

Download PGVTN_mmcdara.mw

PS 1: using PDETools:-declare doesn't ease the task of finding where the error is and you must be doubly careful using it.
PS 2: be consistent and declare u[e] too if you consider it as a function of x (which u[e]' seems to mean, even if u[e] is later defined as a constant function).


As explained in the attached file your code does not work with my Maple 2015?
So I copy-pasted your expression of eq5 and used it as the starting point to get eq6.

Your code is in grey, mine on bold brown; the successive transformations of F (the integral term in your eq6) enableobtaoining a standard form from which it clearly appears that the orthogonality condition of Hermite polynomials (Maple uses their physical form) can be used to simplify the expression of F:

Final remark: as _C2 intervenes through its square, there are two real and opposite values of _C2 such that F(n)=1.

restart;

kernelopts(version)

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

(1)


Here is your file NOT EXECUTED in Maple 2015 for PDEtools:-dchange(...) generates an error with Maple 2015

with(PDEtools):
#with(Maplets[Elements]):

eq0:= E1 = (n+1/2)*_h; param := m1 =_h*sigma^2;

E1 = (n+1/2)*_h

 

m1 = _h*sigma^2

(2)

os := expand(isolate((-_h^2/(2*m1))*diff(psi(x),x,x)+(m1*x^2/2)*psi(x)-E1*psi(x)=0,diff(psi(x),x,x)));

diff(diff(psi(x), x), x) = m1^2*x^2*psi(x)/_h^2-2*m1*E1*psi(x)/_h^2

(3)

eq1 := PDEtools:-dchange([x=sqrt(_h/m1)*zeta,psi(x)=f(zeta)*exp(-zeta^2/2)],os,[zeta,f(zeta)],params={m1,_h,E1});
eq2 := expand(isolate(eq1,diff(f(zeta),zeta$2)));
eq3 := simplify(eval(remove(has,convert(rhs(dsolve(eq2,f(zeta))),Hermite),KummerM),eq0)*exp(-zeta^2/2),radical)       assuming zeta::positive;
eq4 := subs(param,simplify(subs(zeta=sqrt(_h/m1)*x,eq3))) ;
eq5 := unapply(eq4,n,x,sigma);
eq6 := solve(int(eq5(n,x,sigma)^2, x=-infinity...infinity)=1,_C2,useassumptions) assuming _C2 :: real, sigma::positive,n::posint;
 

((diff(diff(f(zeta), zeta), zeta))*exp(-(1/2)*zeta^2)-2*(diff(f(zeta), zeta))*zeta*exp(-(1/2)*zeta^2)-f(zeta)*exp(-(1/2)*zeta^2)+f(zeta)*zeta^2*exp(-(1/2)*zeta^2))*m1/_h = m1*zeta^2*f(zeta)*exp(-(1/2)*zeta^2)/_h-2*m1*E1*f(zeta)*exp(-(1/2)*zeta^2)/_h^2

 

diff(diff(f(zeta), zeta), zeta) = -2*E1*f(zeta)/_h+2*(diff(f(zeta), zeta))*zeta+f(zeta)

 

_C2*HermiteH(n, zeta)*exp(-(1/2)*zeta^2)/2^n

 

_C2*HermiteH(n, (1/sigma^2)^(1/2)*x)*2^(-n)*exp(-(1/2)*x^2/sigma^2)

 

proc (n, x, sigma) options operator, arrow; _C2*HermiteH(n, (1/sigma^2)^(1/2)*x)*2^(-n)*exp(-(1/2)*x^2/sigma^2) end proc

 

(4)


From now on: I copy-paste the expression of "your" (that I cannot obtain with Maple 2015)

eq5 := proc (n, x, sigma) options operator, arrow; _C2*HermiteH(n, sqrt(1/sigma^2)*x)*2^(-n)*exp(-(1/2)*x^2/sigma^2) end proc

proc (n, x, sigma) options operator, arrow; _C2*HermiteH(n, sqrt(1/sigma^2)*x)*2^(-n)*exp(-(1/2)*x^2/sigma^2) end proc

(5)


Here is the integral contained in your relation "eq6"

After some elementary algebra I get the last relation above

F := Int(eq5(n,x,sigma)^2, x=-infinity...infinity) assuming n::posint;
F := combine(F, exp);
F := eval(F, sqrt(1/sigma^2) = 1/sigma);
F := simplify(IntegrationTools:-Expand(F));
F := IntegrationTools:-Change(F, x/sigma=z, z) assuming sigma > 0;
F := simplify(IntegrationTools:-Expand(F));

Int(_C2^2*HermiteH(n, (1/sigma^2)^(1/2)*x)^2*(2^(-n))^2*(exp(-(1/2)*x^2/sigma^2))^2, x = -infinity .. infinity)

 

Int(_C2^2*HermiteH(n, (1/sigma^2)^(1/2)*x)^2*(2^(-n))^2*exp(-x^2/sigma^2), x = -infinity .. infinity)

 

Int(_C2^2*HermiteH(n, x/sigma)^2*(2^(-n))^2*exp(-x^2/sigma^2), x = -infinity .. infinity)

 

_C2^2*(Int(HermiteH(n, x/sigma)^2*exp(-x^2/sigma^2), x = -infinity .. infinity))*4^(-n)

 

_C2^2*(Int(HermiteH(n, z)^2*exp(-z^2)*sigma, z = -infinity .. infinity))*4^(-n)

 

_C2^2*sigma*(Int(HermiteH(n, z)^2*exp(-z^2), z = -infinity .. infinity))*4^(-n)

(6)


From orthogonality of Hermite polynomials (Physical form used in Maple)
(see Wiki)

select(has, indets(F, function), Int)[];
OrthogonalityCondition := % = (2^n*n!*sqrt(Pi))

Int(HermiteH(n, z)^2*exp(-z^2), z = -infinity .. infinity)

 

Int(HermiteH(n, z)^2*exp(-z^2), z = -infinity .. infinity) = 2^n*factorial(n)*Pi^(1/2)

(7)

F := unapply(eval(F, OrthogonalityCondition), n)

proc (n) options operator, arrow; _C2^2*4^(-n)*sigma*2^n*factorial(n)*Pi^(1/2) end proc

(8)


Thus the expression of _C2 for any n

eq6 := unapply( [solve(F(n)=1, _C2)], n)

proc (n) options operator, arrow; [1/(factorial(n)*Pi^(1/2)*2^(-n)*sigma)^(1/2), -1/(factorial(n)*Pi^(1/2)*2^(-n)*sigma)^(1/2)] end proc

(9)

 


 

Download SolnforOS_dsolve_mmcdara.mw



The two derivatives give exactly the same thing a.

Nevertheless here is a way to do what you want derivatives.mw


I_would_do_that.mw

I declared prnt and clr as optional parameters (first one is boolean initialized to true, second one is string initialized to blue), I introduced a control of the color you pass, and I simplified your code a little bit.
Here are the modifications (complete code in the attached file)

spread2:=proc(
               ....
               {clr::string:="Blue"}, 
               {prnt::boolean:=true}
             )
           ....
           local col := substring(StringTools:-Capitalize(clr), 1..1);
           local COL := `if`(col="R", "Red", `if`(col="G", "Green", "Blue"))  
           ....
           if `not`(member(col, {"R", "G", "B"})) then
             error cat("allowed colors are red, green or blue, received ", clr)
           end if:

           if prnt then
             print(cat("Spread 2 [x,y] Points/Vectors wrt origin ",  COL));
           end if:
                 
           if col = "B" then
               return ....

           elif col="G" then
               return ....

           elif col="R" then
               return ....
          end if;
          end proc:



More controls could be necessary for a safer code (what if the denominators in the green and red cases are null?)

restart:

expr := h(arctan(y/x))+arctan(a/b)^2

h(arctan(y/x))+arctan(a/b)^2

(1)

eval(expr, arctan = (u -> arctan((numer, denom)(u))))

h(arctan(y, x))+arctan(a, b)^2

(2)

 

Download arctan.mw

Unfortunately I gave only Maple 2015 and WeibullPlot is not avaliable.

Look to this file and tell me if the plot I draw corresponds or not to what you would like to get?
I'm afraid I did not corerctly understood what is the issue you refer to. 
Here axis 1 is given 4 subticks and axis 2  3 subricks, you can see that the plot is correct.

restart

with(Statistics):

XX := RandomVariable(Weibull(1, .6)):

AA := Sample(XX, 100):

WeibullPlot(AA, reference = false, style = line, gridlines = true, size = [800, 500], axis = [tickmarks = [default, subticks = 4], color = darkgreen])

WeibullPlot(Vector[row](%id = 18446744078086407582), reference = false, style = line, gridlines = true, size = [800, 500], axis = [tickmarks = [default, subticks = 4], color = darkgreen])

(1)

kernelopts(version)

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

(2)

F  := unapply(CDF(Weibull(1, 0.6), z), z):

AA := sort(AA):
N  := numelems(AA):

interface(warnlevel=0):
WB := map(n -> [AA[n], solve(F(z)=n/N)], [$1..N-1]):

plots:-display(
  plot(
    WB
    , color="NavyBlue"
    , gridlines=true
    , axis[1]=[mode=log, color="DarkGreen", tickmarks=[subticks=3]]
    , axis[2]=[mode=log, color="DarkGreen", tickmarks=[subticks=4]]
    , size=[700, 500]
    , axes=boxed
  )
  , plot(Mean(Weibull(1, 0.6)), z=(min..max)(AA), color=red, linestyle=3)
  , plot([[Mean(AA), WB[1][2]], [Mean(AA), WB[-1][2]]], z=(min..max)(AA), color=red, linestyle=3)
)

 

Sample_data    := log[10]~(map2(op, 1, WB)):
Reference_data := log[10]~(map2(op, 2, WB)):
FIT := LinearFit([1, x], Reference_data, Sample_data, [x], output=solutionmodule):
print~(FIT:-Results()):

"residualmeansquare" = 0.483263502868184e-2

 

"residualsumofsquares" = .468765597782138

 

"residualstandarddeviation" = 0.695171563621660e-1

 

"degreesoffreedom" = 97

 

"parametervalues" = (Vector(2, {(1) = 0.566005759094454e-1, (2) = .996156636676824}))

 

"parametervector" = (Vector(2, {(1) = 0.566005759094454e-1, (2) = .996156636676824}))

 

"leastsquaresfunction" = 0.566005759094454e-1+.996156636676824*x

 

Typesetting:-mrow(Typesetting:-ms("standarderrors"), Typesetting:-mo("=", mathvariant = "normal", fence = "false", separator = "false", stretchy = "false", symmetric = "false", largeop = "false", movablelimits = "false", accent = "false", lspace = "0.2777778em", rspace = "0.2777778em"), Typesetting:-maction(Typesetting:-mfenced(Typesetting:-mtable(Typesetting:-mtr(Typesetting:-mtd(Typesetting:-mn("0.00770301199087664", mathvariant = "normal"), rowalign = "", columnalign = "", groupalign = "", rowspan = "1", columnspan = "1"), Typesetting:-mtd(Typesetting:-mn("0.00800405156247162", mathvariant = "normal"), rowalign = "", columnalign = "", groupalign = "", rowspan = "1", columnspan = "1"), rowalign = "", columnalign = "", groupalign = ""), align = "axis", rowalign = "baseline", columnalign = "center", groupalign = "{left}", alignmentscope = "true", columnwidth = "auto", width = "auto", rowspacing = "1.0ex", columnspacing = "0.8em", rowlines = "none", columnlines = "none", frame = "none", framespacing = "0.4em 0.5ex", equalrows = "false", equalcolumns = "false", displaystyle = "false", side = "right", minlabelspacing = "0.8em"), mathvariant = "normal", Typesetting:-msemantics = "RowVector", open = "[", close = "]", Typesetting:-msemantics = "RowVector"), actiontype = "rtableaddress", rtableid = "18446744078229114262"))

 

"confidenceintervals" = (Vector(2, {(1) = 0.413122294090191e-1 .. 0.718889224098718e-1, (2) = .980270809958746 .. 1.01204246339490}))

 

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

 

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

 

"variancecovariancematrix" = (Matrix(2, 2, {(1, 1) = 0.593363937315893e-4, (1, 2) = 0.259631230458869e-4, (2, 1) = 0.259631230458869e-4, (2, 2) = 0.640648414147044e-4}))

 

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

 

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

 

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

 

"AtkinsonTstatistic" = Vector[column](%id = 18446744078229089086)

(3)

p := plots:-display(
  plot([seq([Reference_data[n], Sample_data[n]], n=1..numelems(Sample_data))], color="NavyBlue")
  , plot(
      FIT:-Results("leastsquaresfunction"), x=(min..max)(Reference_data), color="NavyBlue", linestyle=3
    )
  , gridlines=true
):

p := plottools:-transform((x, y) -> [y, x])(p):

plots:-display(
  p
  , axis[1]=[tickmarks=[seq(k=10.0^k, k=-3..1)], color="DarkGreen"]
  , axis[2]=[tickmarks=[seq(k=10.0^k, k=-3..1)], color="DarkGreen"]
  , view=[(min..max)(Sample_data), (min..max)(Reference_data)]
  , size=[700, 500]
  , axes=boxed
)

 

 

Download WeibullPlot_2015.mw

but an interval.

Just do this to see what happens for a given y as x goes to 0.

limit(f(x, y), x=0)
                            /1\             /1\
                   -1 + sin| - | .. 1 + sin| - |
                            \y/             \y/

Maybe doing this will help you understand what happens

limit(f(1/x, 1/y), x=+infinity)
                   -1 + sin(y) .. 1 + sin(y)
limit(sin(x), x=+infinity)
                            -1 .. 1

You can also compute a directional limit 

limit(f(x, lambda*x), x=0)
                            -2 .. 2

In the attached file you will find an exploration of f(x,y) as you get closer and closer to the point (0, 0) (the value h of the slider means the neighborhood of (0, 0) has size 10-h). 
limit.mw


A quick workaround is

plots:-implicitplot(g=1e-8, x = 0 .. 80, y = 0 .. 60, scaling = constrained, grid=[400, 400])


Nevertheless do not gorget that g=0 on the domain  33 <= x <= 58,  0, <= y <= 49.99226668
gnull.mw

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