mmcdara

6326 Reputation

17 Badges

7 years, 348 days

MaplePrimes Activity


These are answers submitted by mmcdara

The attached file contains a numeric exploration of some properties of Eq and an analytic study about its number of real roots.
rho-analysis_mmcdara.mw

As you repliedme once that you wanted analytic results and not results based on numeric simulation or graphic interpretation, this should give you some ideas about how analyzing your equation.

J := n -> int(cos(theta)*LegendreP(n,theta), theta=0..Pi):
J(1)
                               -2
J(2)
                             -3 Pi
J(3);
                               15   2
                          33 - -- Pi 
                               2     
J(n)
      int(cos(theta) LegendreP(n, theta), theta = 0 .. Pi)

 

Replace

evalf(2*Int(F[m,mu],t=0..7*B[m,mu],method = _Dexp)

by

evalf(2*Int(F[m,mu],t=0..7*B[m,mu],method = _d01ajc)

See the attached file (a print is added in the loop which failed to see the progression of the computations).

LoopError_mmcdara.mw

BTW: I don't understand what you are trying to achieve with your three last lines.

You have primitive pythagorean triples (PPT) and non primitive ones.
For instance [3, 4, 5] is a PPT from which some other triples can be found:

  • through permutations you get 6 triples,
  • by changing each number to minus its value (if you accept this) you get 8 triples (thus a total of 48 if you use permutations),
  • by multiplying each number by any positive integer > 1 you get a new triple.

The key is to find the PPT which are in some sense the generators for other pythagorean triples (PT).

Here is a cowhich build PPTs and some developments to generate induced PT.
PithagoreanTriples.mw

My analysis:

Line 1 , is the equivalent of 

Distrib := unapply(PDF(Normal(average, sqrt(var)), 1.1*x-0.1), (x, average, var))

Line 4  creates a sample (size n) of a random variable with a Truncated Normal Distribution (mean=1, variance=0.399, support=0.8..3).
Statistics doesn't have this kind of distribution but it is very easy to create one.
Sampling it is not a problem neither.
Illustration (code written in a hurry): TruncatedNormal.mw

I don't know what Total means?
Could it be

Mean(func~(RV) /~ Distrib~(RV, 1, 0.399))
# where RV is the name of the sample generated at line (6) of the attached file

If it is so, then the correponding Maple code could be this one:
 

restart

with(Statistics):

TruncatedNormal := proc(mu, var, a, b)
   
   local phi   := unapply(PDF(Normal(0, 1), x), x):
   local Phi   := unapply(CDF(Normal(0, 1), x), x):
   local sigma := sqrt(var):
   local xi    := (x-mu)/sigma:
   local alpha := (a-mu)/sigma:
   local beta  := (b-mu)/sigma:
   local Z     := Phi(beta) - Phi(alpha):

   Distribution(
        PDF = unapply( phi(xi)/(sigma*Z), x)
      , CDF = unapply( (Phi(xi) - Phi(alpha))/Z, x)
      , Support = a..b
      , Mode     = piecewise(mu < a, a, mu > b, b, mu)
      , Mean     = mu + (phi(alpha) - phi(beta))/Z*sigma
      , Median   = `if`(
                      a::numeric and b::numeric,
                      mu + Quantile(Normal(0, 1), (Phi(alpha) + Phi(beta))/2, numeric) * sigma,
                      `non numeric`
                   #   mu + (Phi@@-1)(Normal(0, 1), (Phi(alpha) + Phi(beta))/2) * sigma
                   )
      , Variance = sigma^2
                   * (
                       1
                       -
                       (beta*phi(beta) - alpha*phi(alpha)) / Z
                       -
                       ( (phi(alpha) - phi(beta)) / Z )^2
                     )

   )
end proc:
 

X := RandomVariable(TruncatedNormal(1, 0.399, 0.8, 3))

_R2

(1)

supp := Support(X, 'output = range')

.8 .. 3.

(2)

n := 10:

RV := Sample(X, n, method=[envelope, range=supp])

RV := Vector[row](10, {(1) = .897362485588384, (2) = 1.02691084339848, (3) = .935420121873961, (4) = 1.65069973578381, (5) = 1.02359559220548, (6) = .841779027249969, (7) = .928454628140565, (8) = 1.30070117114770, (9) = 1.00227755594279, (10) = 1.78564814703848})

(3)

func := x -> 1/(1 + sinh(2*x)*log(x)^2);

proc (x) options operator, arrow; 1/(1+sinh(2*x)*log(x)^2) end proc

(4)

Distrib := unapply( PDF(Normal(average, sqrt(var)), 1.1*x - 0.1), (x, average, var));

proc (x, average, var) options operator, arrow; (1/2)*2^(1/2)*exp(-(1/2)*(1.1*x-.1-average)^2/var)/(Pi^(1/2)*var^(1/2)) end proc

(5)

aux_1 := Mean(func~(RV) /~ Distrib~(RV, 1, 0.399))

HFloat(1.3470675276101085)

(6)

aux_2 := int(Distrib(x, 1, 0.399), x = supp)

.5781266537

(7)

_Int := evalf(aux_1*aux_2)

HFloat(0.7787756420451645)

(8)

# A larger X sample

n  := 10^5:
RV := Sample(X, n, method=[envelope, range=supp]):

# Check we correctly sampled X

plots:-display(
  Histogram(RV, minbins=100, style=polygon, transparency=0.4),
  plot(PDF(X, x), x=supp, thickness=3, color=red)
);

 

aux_1 := Mean(func~(RV) /~ Distrib~(RV, 1, 0.399)):
aux_2 := int(Distrib(x, 1, 0.399), x = supp):
_Int := evalf(aux_1*aux_2)

HFloat(0.6580071733702313)

(9)

ExactFormula    := Int(func(x), x=supp);
QuasiExactValue := value(%);

Int(1/(1+sinh(2*x)*ln(x)^2), x = .8 .. 3.)

 

.6768400757

(10)

 


Download TruncatedNormal_2.mw

Here is an example of the uncolding of the prism.
All the steps up to the construction of the Spanning Tree should be generic (at least for a convex polyhedron).
The last one (the drawing of the pattern of this prism) is done by hand.
I didn't try ro make this figure by using the true coordinates of your prism (for the method to apply see Here for instance).

As I don't have time to go any further, I hope that someone else here will be able to take advantage of this code to produce something more substantial.

restart

with(GraphTheory):
with(plots):
with(plottools):

T := polygon([[0, 0], [2, 1], [1, 3]])

POLYGONS([[0., 0.], [2., 1.], [1., 3.]])

(1)

solid := display(prism(T, height = 3)):
faces := map(u -> select(type, u, Matrix)[], [getdata(solid)]);
NF := numelems(faces)

faces := [Matrix(3, 3, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (2, 1) = 2.0, (2, 2) = 1.0, (2, 3) = .0, (3, 1) = 1.0, (3, 2) = 3.0, (3, 3) = .0}, datatype = float[8]), Matrix(4, 3, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = 3.0, (3, 1) = 2.0, (3, 2) = 1.0, (3, 3) = 3.0, (4, 1) = 2.0, (4, 2) = 1.0, (4, 3) = .0}, datatype = float[8]), Matrix(4, 3, {(1, 1) = 2.0, (1, 2) = 1.0, (1, 3) = .0, (2, 1) = 2.0, (2, 2) = 1.0, (2, 3) = 3.0, (3, 1) = 1.0, (3, 2) = 3.0, (3, 3) = 3.0, (4, 1) = 1.0, (4, 2) = 3.0, (4, 3) = .0}, datatype = float[8]), Matrix(4, 3, {(1, 1) = 1.0, (1, 2) = 3.0, (1, 3) = .0, (2, 1) = 1.0, (2, 2) = 3.0, (2, 3) = 3.0, (3, 1) = .0, (3, 2) = .0, (3, 3) = 3.0, (4, 1) = .0, (4, 2) = .0, (4, 3) = .0}, datatype = float[8]), Matrix(3, 3, {(1, 1) = 1.0, (1, 2) = 3.0, (1, 3) = 3.0, (2, 1) = 2.0, (2, 2) = 1.0, (2, 3) = 3.0, (3, 1) = .0, (3, 2) = .0, (3, 3) = 3.0}, datatype = float[8])]

 

5

(2)

faces := convert~(faces, listlist):

HasCommonEdge := proc(edge)
  local n, s, c:
  c := NULL:
  for n from 1 to NF do
    s := select(has, faces[n], edge):
    if s <> [] then
      if numelems(s) = 2 then c := c, n end if:
    end if:
  end do:
  return cat~(F, {c})
end proc:

ConnectivityByEdge := NULL:

for n from 1 to NF do
  L := [$1..numelems(faces[n]), 1]:
  for i from 1 to numelems(L)-1 do
    edge := faces[n][[L[i], L[i+1]]];
    ConnectivityByEdge := ConnectivityByEdge, HasCommonEdge(edge);
  end do:
end do:

ConnectivityByEdge := {ConnectivityByEdge}
    

{{F1, F2}, {F1, F3}, {F1, F4}, {F2, F3}, {F2, F4}, {F2, F5}, {F3, F4}, {F3, F5}, {F4, F5}}

(3)

G_Vertices := SpecialGraphs:-PrismGraph(3);
G_Faces    := Graph(ConnectivityByEdge):

DocumentTools:-Tabulate(
  [
    DrawGraph(G_Vertices, style=default, dimension=3),
    DrawGraph(G_Faces, style=default, dimension=3)
  ]
  , width=50
):

GRAPHLN(undirected, unweighted, [1, 2, 3, 4, 5, 6], Array(%id = 18446744078276146710), `GRAPHLN/table/1`, 0)

(4)

# How the faces are connected in the pattern of the prism?

ST := SpanningTree(G_Faces);
DrawGraph(ST)

GRAPHLN(undirected, unweighted, [F1, F2, F3, F4, F5], Array(%id = 18446744078367490526), `GRAPHLN/table/6`, 0)

 

 

# An illustration (note that face F5 cannot be connected to another edge of F2 because it
# shares no vertex with F1

face1 := display(polygon([[0, 0], [1, 0], [1/2, sqrt(3)/2]]), color=red):
face2 := display(polygon([[0, 0], [1, 0], [1, -1], [0, -1]]), color=cyan):
face3 := display(polygon([[1, 0], [1+sqrt(3)/2, 1/2], [1/2+sqrt(3)/2, sqrt(3)/2+1/2], [1/2, sqrt(3)/2]]), color=yellow):
face4 := display(polygon([[0, 0], [-sqrt(3)/2, 1/2], [1/2-sqrt(3)/2, sqrt(3)/2+1/2], [1/2, sqrt(3)/2]]), color=green):
face5 := display(polygon([[0, -1], [1, -1], [1/2, -1-sqrt(3)/2]]), color=blue):

display(
  face1, textplot([1/2, sqrt(3)/4, F1], color=white, font=[Tahoma, bold]),
  face2, textplot([1/2, -1/2, F2], color=black, font=[Tahoma, bold]),
  face3, textplot([(3+sqrt(3))/4, (1+sqrt(3))/4, F3], color=black, font=[Tahoma, bold]),
  face4, textplot([(1-sqrt(3))/4, (1+sqrt(3))/4, F4], color=black, font=[Tahoma, bold]),
  face5, textplot([1/2, -1-sqrt(3)/4, F5], color=white, font=[Tahoma, bold]),
  scaling=constrained,
  axes=none
)

 

# Another spanning tree

ST := SpanningTree(G_Faces, F2);
DrawGraph(ST, style=tree, root=F2)

GRAPHLN(undirected, unweighted, [F1, F2, F3, F4, F5], Array(%id = 18446744078367481854), `GRAPHLN/table/11`, 0)

 

 

face1 := display(polygon([[0, 0], [1, 0], [1/2, sqrt(3)/2]]), color=red):
face2 := display(polygon([[0, 0], [1, 0], [1, -1], [0, -1]]), color=cyan):
face3 := display(polygon([[0, 0], [-1, 0], [-1, -1], [0, -1]]), color=yellow):
face4 := display(polygon([[1, 0], [2, 0], [2, -1], [0, -1]]), color=green):
face5 := display(polygon([[0, -1], [1, -1], [1/2, -1-sqrt(3)/2]]), color=blue):

display(
  face1, textplot([1/2, sqrt(3)/4, F1], color=white, font=[Tahoma, bold]),
  face2, textplot([1/2, -1/2, F2], color=black, font=[Tahoma, bold]),
  face3, textplot([-1/2, -1/2, F3], color=black, font=[Tahoma, bold]),
  face4, textplot([3/2, -1/2, F4], color=black, font=[Tahoma, bold]),
  face5, textplot([1/2, -1-sqrt(3)/4, F5], color=white, font=[Tahoma, bold]),
  scaling=constrained,
  axes=none
)

 

# There are

NumberOfSpanningTrees(G_Faces);

# ways to unfold the prism (some are identical to within one symmetry)
 

75

(5)
 

 

Download A_first_sketch_for_a_code.mw

restart

expr := 3*G*(`&Delta;&gamma;`*H - `&sigma;y`(`&Delta;&gamma;`) + q)/(-q + `&sigma;y`(`&Delta;&gamma;`))^2;

3*G*(`&Delta;&gamma;`*H-`&sigma;y`(`&Delta;&gamma;`)+q)/(-q+`&sigma;y`(`&Delta;&gamma;`))^2

(1)

indets(expr, function)[];
rw := % = Z+q;
rw:= isolate(op(1, denom(expr))=Z, `&sigma;y`(`&Delta;&gamma;`))

`&sigma;y`(`&Delta;&gamma;`)

 

`&sigma;y`(`&Delta;&gamma;`) = Z+q

 

`&sigma;y`(`&Delta;&gamma;`) = Z+q

(2)

algsubs(rw, expr)

3*G*(H*`&Delta;&gamma;`-Z)/Z^2

(3)

expand(%);

3*G*`&Delta;&gamma;`*H/Z^2-3*G/Z

(4)

res := eval(%, isolate(rw, Z)) assuming Z > 0

3*G*`&Delta;&gamma;`*H/(-q+`&sigma;y`(`&Delta;&gamma;`))^2-3*G/(-q+`&sigma;y`(`&Delta;&gamma;`))

(5)

sort~(res, q, ascending)

3*G*`&Delta;&gamma;`*H/(`&sigma;y`(`&Delta;&gamma;`)-q)^2-3*G/(`&sigma;y`(`&Delta;&gamma;`)-q)

(6)
 

 

Download step_by_step.mw


@acer  and I intertwined: this is basically the same answer but I suggest you to dsolve WITHOUT giving Cl a numeric value.
The solution then depends on t AND Cl.
 

restart;

st := time():

N := 4:

T := 5.0/60:

M := 6905:

dose := t -> M*sum(Heaviside(t - 24*k/N)*exp((-t + 24*k/N)/T), k = 0 .. 2*N - 1):

deq  := diff(C(t), t) = dose(t) - Cl*C(t):

dsolve({deq, C(0) = 0}):
value(%) assuming Cl > 0:
sol := unapply(rhs(%), (t, Cl));

time() - st;

proc (t, Cl) options operator, arrow; (6905*Heaviside(t)*exp(t*(Cl-12))/(Cl-12)-6905*Heaviside(t)/(Cl-12)+6905*Heaviside(t-6)*exp(Cl*t-12*t+72)/(Cl-12)-6905*Heaviside(t-6)*exp(6*Cl)/(Cl-12)+6905*Heaviside(t-12)*exp(Cl*t-12*t+144)/(Cl-12)-6905*Heaviside(t-12)*exp(12*Cl)/(Cl-12)+6905*Heaviside(t-18)*exp(Cl*t-12*t+216)/(Cl-12)-6905*Heaviside(t-18)*exp(18*Cl)/(Cl-12)+6905*Heaviside(t-24)*exp(Cl*t-12*t+288)/(Cl-12)-6905*Heaviside(t-24)*exp(24*Cl)/(Cl-12)+6905*Heaviside(t-30)*exp(Cl*t-12*t+360)/(Cl-12)-6905*Heaviside(t-30)*exp(30*Cl)/(Cl-12)+6905*Heaviside(t-36)*exp(Cl*t-12*t+432)/(Cl-12)-6905*Heaviside(t-36)*exp(36*Cl)/(Cl-12)+6905*Heaviside(t-42)*exp(Cl*t-12*t+504)/(Cl-12)-6905*Heaviside(t-42)*exp(42*Cl)/(Cl-12))*exp(-Cl*t) end proc

 

.100

(1)

plots:-display(
  plot(
    [sol(t, 0.32), sol(t, 0.45)]
    , t = 0 .. 48
    , color=[red, blue]
    , legend=[typeset(Cl=0.32), typeset(Cl=0.45)]
 )
);

 

 


 

Download Analytic_solution.mw

At the end you get three different sets of solutions:

restart

u(xi):=c[0]+c[1]*a^(f(xi))

c[0]+c[1]*a^f(xi)

(1)

Eq1 := -(alpha^2 + beta)*u(xi) + k^2*diff(u(xi), xi $ 2) + b*u(xi)^3 = 0

-(alpha^2+beta)*(c[0]+c[1]*a^f(xi))+k^2*(c[1]*a^f(xi)*(diff(f(xi), xi))^2*ln(a)^2+c[1]*a^f(xi)*(diff(diff(f(xi), xi), xi))*ln(a))+b*(c[0]+c[1]*a^f(xi))^3 = 0

(2)

 

NULL

NULL

NULL

D1 := diff(f(xi), xi) = (a^(-f(xi))*p + q + r*a^f(xi))/ln(a);

diff(f(xi), xi) = (a^(-f(xi))*p+q+r*a^f(xi))/ln(a)

(3)

D2 := diff(f(xi), xi$2) = eval(diff(rhs(D1), xi), D1);

diff(diff(f(xi), xi), xi) = (-a^(-f(xi))*(a^(-f(xi))*p+q+r*a^f(xi))*p+r*a^f(xi)*(a^(-f(xi))*p+q+r*a^f(xi)))/ln(a)

(4)

DEq1_0 := eval(Eq1, {D1, D2})

-(alpha^2+beta)*(c[0]+c[1]*a^f(xi))+k^2*(c[1]*a^f(xi)*(a^(-f(xi))*p+q+r*a^f(xi))^2+c[1]*a^f(xi)*(-a^(-f(xi))*(a^(-f(xi))*p+q+r*a^f(xi))*p+r*a^f(xi)*(a^(-f(xi))*p+q+r*a^f(xi))))+b*(c[0]+c[1]*a^f(xi))^3 = 0

(5)

# Setting Z = a^f(xi)

DEq1_1 := collect(eval(expand(DEq1_0), a^f(xi) = Z), Z)

(2*k^2*r^2*c[1]+b*c[1]^3)*Z^3+(3*k^2*q*r*c[1]+3*b*c[0]*c[1]^2)*Z^2+(2*k^2*p*r*c[1]+k^2*q^2*c[1]+3*b*c[0]^2*c[1]-alpha^2*c[1]-beta*c[1])*Z+k^2*c[1]*p*q+b*c[0]^3-alpha^2*c[0]-beta*c[0] = 0

(6)

# I don't understand why yhis command doesn't work ???
coeffs(DEq1_1, Z);

Error, invalid arguments to coeffs

 

# Workaround

with(LargeExpressions):

collect(DEq1_1, Z, Veil[C] );
CoefficientNullity := [ seq( 0 = Unveil[C]( C[i] ), i=1..LastUsed[C] ) ];
sols := {solve(CoefficientNullity)}:

print(cat('_'$120)):
print~(sols):

Z^3*C[1]+3*Z^2*C[2]-Z*C[3]-C[4] = 0

 

[0 = 2*k^2*r^2*c[1]+b*c[1]^3, 0 = k^2*q*r*c[1]+b*c[0]*c[1]^2, 0 = -2*k^2*p*r*c[1]-k^2*q^2*c[1]-3*b*c[0]^2*c[1]+alpha^2*c[1]+beta*c[1], 0 = -k^2*p*q*c[1]-b*c[0]^3+alpha^2*c[0]+beta*c[0]]

 

________________________________________________________________________________________________________________________

 

{alpha = alpha, b = 0, beta = -alpha^2, k = 0, p = p, q = q, r = r, c[0] = c[0], c[1] = c[1]}

 

{alpha = alpha, b = 0, beta = -alpha^2, k = k, p = p, q = 0, r = 0, c[0] = c[0], c[1] = c[1]}

 

{alpha = alpha, b = 0, beta = k^2*q^2-alpha^2, k = k, p = p, q = q, r = 0, c[0] = c[1]*p/q, c[1] = c[1]}

 

{alpha = alpha, b = b, beta = beta, k = k, p = p, q = q, r = r, c[0] = 0, c[1] = 0}

 

{alpha = alpha, b = (alpha^2+beta)/c[0]^2, beta = beta, k = k, p = p, q = q, r = r, c[0] = c[0], c[1] = 0}

 

{alpha = alpha, b = -2*k^2*r^2/c[1]^2, beta = 2*k^2*p*r-(1/2)*k^2*q^2-alpha^2, k = k, p = p, q = q, r = r, c[0] = (1/2)*q*c[1]/r, c[1] = c[1]}

(7)

# After removing solutions which contain b=0 one gets three

sols := convert(remove((x -> member(b=0, x)), sols), list):
print~(sols):

{alpha = alpha, b = b, beta = beta, k = k, p = p, q = q, r = r, c[0] = 0, c[1] = 0}

 

{alpha = alpha, b = (alpha^2+beta)/c[0]^2, beta = beta, k = k, p = p, q = q, r = r, c[0] = c[0], c[1] = 0}

 

{alpha = alpha, b = -2*k^2*r^2/c[1]^2, beta = 2*k^2*p*r-(1/2)*k^2*q^2-alpha^2, k = k, p = p, q = q, r = r, c[0] = (1/2)*q*c[1]/r, c[1] = c[1]}

(8)

# Of course you need a few other assumptions to remove all he solutions for which b could be <=0:

map(x -> select(has, x, b), sols);
b_PositivityConditions := map(x -> solve(rhs(x[1]) > 0), %);
 

[{b = b}, {b = (alpha^2+beta)/c[0]^2}, {b = -2*k^2*r^2/c[1]^2}]

 

[RealRange(Open(0), infinity), {alpha = alpha, c[0] < 0, -alpha^2 < beta}, {alpha = alpha, 0 < c[0], -alpha^2 < beta}]

(9)

# Final solutions:

FinalSolutions := [
                    `union`(sols[1], {b in b_PositivityConditions[1]}),
                    `union`(sols[2], b_PositivityConditions[2]),
                    `union`(sols[3], b_PositivityConditions[3])
                  ]:

print~(FinalSolutions):

{`in`(b, RealRange(Open(0), infinity)), alpha = alpha, b = b, beta = beta, k = k, p = p, q = q, r = r, c[0] = 0, c[1] = 0}

 

{alpha = alpha, b = (alpha^2+beta)/c[0]^2, beta = beta, k = k, p = p, q = q, r = r, c[0] = c[0], c[1] = 0, c[0] < 0, -alpha^2 < beta}

 

{alpha = alpha, b = -2*k^2*r^2/c[1]^2, beta = 2*k^2*p*r-(1/2)*k^2*q^2-alpha^2, k = k, p = p, q = q, r = r, c[0] = (1/2)*q*c[1]/r, c[1] = c[1], 0 < c[0], -alpha^2 < beta}

(10)
 

 

``

Download maple_sheet_mmcdara.mw

The case of spherical coordinates is more complex (read the reference I provided you in y previous reply).

restart

with(plots):

f := (x^2+y^2+z^2+12)^2 - 64*(x^2+y^2)

(x^2+y^2+z^2+12)^2-64*x^2-64*y^2

(1)


PAPPUS' THEOREM

crossection_radius := 2:
crossection_area   := Pi*crossection_radius^2:
Rotation_radius    := 4:
Torus_volume       := crossection_area*(2*Pi*Rotation_radius)

32*Pi^2

(2)


TRIPLE INTEGRAL IN CYLINDRICAL COORDINATES


The cut of the torus (considered here as a volume) at a given height z is an annulus (in green in the figure below). whose inner radius Ri and outer
radius Ro both depend on z.
In polar coordinates the surface of this annulus is easy to write as a double integral in, lets say, rho and theta, where rho ranges from Ri(z) to Ro(z).
Think to this cut as a slice of height dz.
To get the volume of the torus just integrate the slice surface from z=-2 to z=2.

R := Rotation_radius:
r := z -> sqrt(crossection_radius^2-z^2):

slice := proc(Z)
  local L, H, G, a:
  uses plottools:
  L := 6:
  H := 2:
  G := 20:
  a := display(
         plottools:-annulus([0, 0], R-r(Z)..R+r(Z), color="Chartreuse"),
         plottools:-circle([0, 0], R-r(Z), color=red, thickness=3),
         plottools:-circle([0, 0], R+r(Z), color=red, thickness=3)
       ):

  display(
    implicitplot3d(f, x=-L..L, y=-L..L, z=-H..Z, grid=[G, G, G], style=surface, color=gray, scaling=constrained),
    translate(extrude(a, 0..0.001), 0, 0, Z)
  )
end proc:


Explore(slice(z), parameters=[z=-1.999..1.999], initialvalues=[z=-1.999])

 

SliceSurface := z -> Int(rho, theta=0..2*Pi, rho=R-r(z)..R+r(z)):
'SliceSurface(z)' = SliceSurface(z);

SliceSurface(z) = Int(rho, theta = 0 .. 2*Pi, rho = 4-(-z^2+4)^(1/2) .. 4+(-z^2+4)^(1/2))

(3)

TorusVolume := Int(SliceSurface(z), z=-2..2);

value(%)

Int(Int(rho, theta = 0 .. 2*Pi, rho = 4-(-z^2+4)^(1/2) .. 4+(-z^2+4)^(1/2)), z = -2 .. 2)

 

32*Pi^2

(4)


TRIPLE INTEGRAL IN CARTESIAN COORDINATES


Proceed as before but express the surface if the slice in cartesian coordinates.
Its better to use symmetry and integrate over que quadrant x >= 0, y >= 0.

The first integral in CartesianSlice expresses the surface of the disk of radius Ro(z) while the second measures the disk of radius Ri(z).

CartesianSlice := z -> 4*( Int(Int(1, y=0..sqrt((R+r(z))^2-x^2)), x=0..R+r(z)) - Int(Int(1, y=0..sqrt((R-r(z))^2-x^2)), x=0..R-r(z)) )

proc (z) options operator, arrow; 4*(Int(Int(1, y = 0 .. sqrt((R+r(z))^2-x^2)), x = 0 .. R+r(z)))-4*(Int(Int(1, y = 0 .. sqrt((R-r(z))^2-x^2)), x = 0 .. R-r(z))) end proc

(5)

TorusVolume := add(map(Int, [op(CartesianSlice(z))], z=-2..2));

Int(4*(Int(Int(1, y = 0 .. ((4+(-z^2+4)^(1/2))^2-x^2)^(1/2)), x = 0 .. 4+(-z^2+4)^(1/2))), z = -2 .. 2)+Int(-4*(Int(Int(1, y = 0 .. ((4-(-z^2+4)^(1/2))^2-x^2)^(1/2)), x = 0 .. 4-(-z^2+4)^(1/2))), z = -2 .. 2)

(6)


As Maple doesn't seem capable to compute this integral (unless changing the integration variables to get back to cylindrical coordinates,
which is obviously of  no interest here), one can proceed numericallly and check the result is the correct one.

# Volume comprised within Ro(z)

IntegrationTools:-Expand(op(1, TorusVolume));
IntegrationTools:-CollapseNested(op(2, %));
OuterVolume := 4*evalf(%);

4*(Int(Int(Int(1, y = 0 .. (-z^2+20-x^2+8*(-z^2+4)^(1/2))^(1/2)), x = 0 .. 4+(-z^2+4)^(1/2)), z = -2 .. 2))

 

Int(1, [y = 0 .. (-z^2+20-x^2+8*(-z^2+4)^(1/2))^(1/2), x = 0 .. 4+(-z^2+4)^(1/2), z = -2 .. 2])

 

392.4859219

(7)

# Volume comprised within Ri(z)

IntegrationTools:-Expand(op(2, TorusVolume));
IntegrationTools:-CollapseNested(op(2, %));
InnerVolume := 4*evalf(%);

-4*(Int(Int(Int(1, y = 0 .. (-z^2+20-x^2-8*(-z^2+4)^(1/2))^(1/2)), x = 0 .. 4-(-z^2+4)^(1/2)), z = -2 .. 2))

 

Int(1, [y = 0 .. (-z^2+20-x^2-8*(-z^2+4)^(1/2))^(1/2), x = 0 .. 4-(-z^2+4)^(1/2), z = -2 .. 2])

 

76.65858104

(8)

# Volume comprised between Ro(z) and Ri(z)

Volume := OuterVolume - InnerVolume;
identify(%);

315.8273409

 

32*Pi^2

(9)
 

 

Download torus_mmcdara.mw


First point: some stramge characters in the definition of P.
Second point: the variable x0 is never defined I guesses it was a type and arbitrarily replaced it by x.

Onece these corrections done:

  • Your design of experiment does not seem correctly constructed (if I correctly understood what you want to achieve)... I corrected it according to my understanding.
  • Fop can be integretated formally once for all. Thus there is no need to write Fop := unapply(...)and I feel it better to write 
    intFop := unapply(...)where intFopis the indefinite integral of Fopwrt x.


Note that the solution is trivial because IntFop is a sum of terms which are all indidually equal to 0 when x=0:

intFop(entries(titles, nolist));
eval(%, x=0)
                               0


Here is the code

restart

with(Student[Calculus1]):

P := -9*F*k*(2*sqrt(2)*(x^2+2)^2*((1/6)*x^2-1)*arctan((1/2)*x*sqrt(2))+Pi*(x^2+2)^2*((1/6)*x^2-1+epsilon*m*Br)*sqrt(2)+(4*(((1/6)*x^2-1)*x^2+5*x^2*(1/9)-14/9))*x*beta)/(16*(x^2+2)^2)+1/((2/3+(1/3)*x^2)*N) = 0;

-(9/16)*F*k*(2*2^(1/2)*(x^2+2)^2*((1/6)*x^2-1)*arctan((1/2)*x*2^(1/2))+Pi*(x^2+2)^2*((1/6)*x^2-1+epsilon*m*Br)*2^(1/2)+4*(((1/6)*x^2-1)*x^2+(5/9)*x^2-14/9)*x*beta)/(x^2+2)^2+1/((2/3+(1/3)*x^2)*N) = 0

(1)

titles := Vector[row]([k, epsilon,  Br, beta, m,F, N, result]);

titles := Vector[row](8, {(1) = k, (2) = epsilon, (3) = Br, (4) = beta, (5) = m, (6) = F, (7) = N, (8) = result})

(2)

data := table([k=[seq(.1 .. .9, .1)], epsilon=[0.5$9], Br=[0.1$9], beta=[1$9], m=[1.5$9],F=[1.5$9], N=[0.5$9]]);

table( [( m ) = [1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5], ( epsilon ) = [.5, .5, .5, .5, .5, .5, .5, .5, .5], ( F ) = [1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5], ( k ) = [.1, .2, .3, .4, .5, .6, .7, .8, .9], ( N ) = [.5, .5, .5, .5, .5, .5, .5, .5, .5], ( Br ) = [.1, .1, .1, .1, .1, .1, .1, .1, .1], ( beta ) = [1, 1, 1, 1, 1, 1, 1, 1, 1] ] )

(3)

DOE := `<|>`( seq(<data[i]>, i in titles[1..-2]) )

DOE := Matrix(9, 7, {(1, 1) = .1, (1, 2) = .5, (1, 3) = .1, (1, 4) = 1, (1, 5) = 1.5, (1, 6) = 1.5, (1, 7) = .5, (2, 1) = .2, (2, 2) = .5, (2, 3) = .1, (2, 4) = 1, (2, 5) = 1.5, (2, 6) = 1.5, (2, 7) = .5, (3, 1) = .3, (3, 2) = .5, (3, 3) = .1, (3, 4) = 1, (3, 5) = 1.5, (3, 6) = 1.5, (3, 7) = .5, (4, 1) = .4, (4, 2) = .5, (4, 3) = .1, (4, 4) = 1, (4, 5) = 1.5, (4, 6) = 1.5, (4, 7) = .5, (5, 1) = .5, (5, 2) = .5, (5, 3) = .1, (5, 4) = 1, (5, 5) = 1.5, (5, 6) = 1.5, (5, 7) = .5, (6, 1) = .6, (6, 2) = .5, (6, 3) = .1, (6, 4) = 1, (6, 5) = 1.5, (6, 6) = 1.5, (6, 7) = .5, (7, 1) = .7, (7, 2) = .5, (7, 3) = .1, (7, 4) = 1, (7, 5) = 1.5, (7, 6) = 1.5, (7, 7) = .5, (8, 1) = .8, (8, 2) = .5, (8, 3) = .1, (8, 4) = 1, (8, 5) = 1.5, (8, 6) = 1.5, (8, 7) = .5, (9, 1) = .9, (9, 2) = .5, (9, 3) = .1, (9, 4) = 1, (9, 5) = 1.5, (9, 6) = 1.5, (9, 7) = .5})

(4)

#DOE := <seq(Vector(9, data[i]), i in titles[1 .. -2])>

Fop := lhs(P):
intFop := unapply(int(Fop, x), [entries(titles, nolist)]);

proc (k, epsilon, Br, beta, m, F, N, result) options operator, arrow; -(1/16)*F*k*2^(1/2)*arctan((1/2)*x*2^(1/2))*x^3+(9/8)*F*k*2^(1/2)*arctan((1/2)*x*2^(1/2))*x+(1/16)*F*k*x^2-(5/4)*F*k*ln((1/2)*x^2+1)-(9/16)*F*k*2^(1/2)*Pi*epsilon*m*Br*x-(1/32)*F*k*2^(1/2)*Pi*x^3-(3/16)*F*k*x^2*beta+(5/4)*F*k*ln((1/2)*x^2+1)*beta+(9/16)*F*k*2^(1/2)*Pi*x+(3/2)*2^(1/2)*arctan((1/2)*x*2^(1/2))/N end proc

(5)

V := Vector(9, i -> intFop(entries(DOE[i], nolist))):

NewtonsMethod~(V, x=1)

Vector[column]([[0.], [0.3e-21], [0.2e-21], [0.1e-21], [0.4e-22], [0.1e-22], [0.1e-22], [0.], [0.1e-23]])

(6)

  


 

Download Help_Newton_Method_mmcdara.mw


I mention Maple 2015 in the title for Maple 2019 and later use the notion of DataFrame which could answer your problem in a simpler way (or maybe not?).

To fix the idea I simulated a N by Q data set where N=100 individuals answered Q=4 questions:

  1. the first has 2 possible outcomes,
  2. so has the second,
  3. the third has 3 possible outcomes,
  4. and the last one has 2 possible outcomes.
     

I generate a 4 dimensional array, named PivotArray, whose each index contains the number of individuals matching this index.
For instance if 10 individuals answered 1 to quastion 1, 1 to question 2, 3 to question 3 and 2 to question 4, then PivotArray[1, 1, 3, 2] = 10.
Once done it remains you several possibilities to build a 2D array which represents the pivot table you need.
For instance you can "group" questions 1 and 2 on one side and questions 3 and 4 on the other side (then your pivot table will have dimensions (4, 6)).
Or group questions 1, 2 and 4 and left question 3 apart  (and then your pivot table will have dimensions (8, 3)).

The attached file build the 2D array given the 4D array and a particular "grouping".

A tool to extract all the indicuals wich give predefined answers is proposed.

Finally a "smart" presentation of the 2D pivot array is proposed (be extremely careful because this part is not completely generic).
More of this, if you use Maple >= 2108, replace

InsertContent(Worksheet(Group(Input( T )))):

by

InsertContent(Worksheet( T )):


If you need more help I advice you to provide a more detailed example of your data (how many outomes has each question).

Last by not least: I never use Statistics:-ChiSquareIndependenceTest (I use to use R instead for this kind of computations [more poreful and versatile]), so I don't know how Maple behaves.

PivotTable0.mw


As my Maple 2015 doesn't accept polygonbyname("rectangle", ...) I replaced it by rectangle(...).
But this should work with your newer Maple's version:
 

restart;
with(plots); with(plottools);
L := point([2, 2], color = black, symbol = solidcircle, symbolsize = 30);
K := line([2, 2], [1, 2], color = red, linestyle = dash, thickness = 4);
K1 := line([2, 2], [2+1.5*((1/2)*sqrt(2)), 2+1.5*((1/2)*sqrt(2))], color = blue, linestyle = dash, thickness = 4);
K3 := rectangle([0, 0], [4, 4], color = "Red");
K2 := annulus([2, 2], 1 .. 3/2, color = blue, style = 'polygon', transparency = .8); 
K4 := disk([2, 2], 3/2, color = white, style = 'polygon', transparency = 0);
display(L, K, K1, K2, K4, K3, size = [300, 300], axes = none);

Couronne_mmcdara.mw


 


 

restart
kernelopts(version)
     Maple 2015.2, APPLE UNIVERSAL OSX, Dec 20 2015, Build ID 1097895

h := (t, x) -> (m*t^2 + 6*t - 2*x)^2/(36*g*t^2):  # Your syntax is not supported by Maple 2015

# Answer to your question
T := [1,5,10,15,20]:
Vector(h~(T, x)):


# More generally: evaluation of h(t, x) on the grid T times X:
X := [$1..5]:
TX := Matrix(numelems(T), numelems(X), (i, j) -> (T[i], X[j])):
h~(TX):


A_solution.mw

..., this could answer your problem  (content of the file can't be uploaded)

proposal_mmcdara.mw

The trick is that its better to solve A=0 wrt to lambda[1] (which expresses lambda[1] as a rational function of x) and next to fit this relation by a simple model that you can easily "invert" to get x as a function of lambda[1] if you prefer.

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