mmcdara

First fix your errors to have a well po...

Answered: mmcdara 7037

May 18 2024

1 9

Next proceed as described in the attached file:
Ode_New_TWO_phase_mmcdara.mw

Advice...
Odes 2 and 4 both form a subsystem with unknown functions h(Theta) and H(Teta) that you should solve apart.
Note that it can be solved numerically or formally.
Once done, plug this soution into odes 1 and 3 (it's better to have solved the first substystem formally, but you can do the same using the option known [see help(dsolve[numeric]) ] if you solved the first subsystem numerically).

Errors to fix:

The first subsystem requires 3 ic/bc and you give 4.
The first subsystem requires 5 ic/bc and you give 4.

@Prakash J Look at this file D...

Answered: mmcdara 7037

May 17 2024

0 0

@Prakash J

Look at this file Difference_mmcdara.mw

@C_R Feel free to delete or reappropriate this comment if you consider that I have interfered inappropriately in your exchanges (I am no longer interested in this race for badges which seems to me to pollute Mapleprimes :-) )

To moderators: this post is deliberately a comment and not an answer, so please don't convert it into an answer. Thanks.

Bad use of D[2]...

Answered: mmcdara 7037

May 13 2024

0 1

First point: do not use D[2] (see red highlighted comment in the attached file).

Once corrected, your expression eq is not identically null (for this to happen the constant A, A2, B1, B2, lambda1 and lambda2 must verify some conditions.
This becomes clear if you consider that, for aany strictlyconfinupus function f, shift(f(n), n) - f(n) is a forward divided difference approximation of diff(f(z), z) at z=n.

verif_May24_mmcdara.mw

Note to the one who converted this comment into an answer:
Had I wanted to write a reply instead of a comment, I would have done so. It wasn't a mistake on my part, so thank you for not changing again the status of this post.

A proposal...

Answered: mmcdara 7037

May 09 2024

0 0

Download Proposal.mw

[added by moderator, using Maple 2022.2]
(Thanks to the moderator, whoever he/she is)

>

restart

>	with(plots):

>	eq := x^3-4xy+2*y^3 = 0

>	p1 := implicitplot(eq, x=-3..3, y=-3..3, grid=[100, 100], color=blue): p2 := plot([[1, -3], [1, 3]], color=red): display(p1, p2)

>

# Find points of coordinates (1, y) which belon to the curve

Y := [solve(eval(eq, x=1))];

# By default all the roots of a 3rd degree equations appear in complex form,
# even they are both real.
# Tp get their real representation transform them in trigonometric expressions:

Y := simplify(convert(Y,trig));

# Floating point approximations.

Yn := evalf(Y)

[(1/6)*(-54+(6*I)*303^(1/2))^(1/3)+4/(-54+(6*I)*303^(1/2))^(1/3), -(1/12)*(-54+(6*I)*303^(1/2))^(1/3)-2/(-54+(6*I)*303^(1/2))^(1/3)+((1/2)*I)*3^(1/2)*((1/6)*(-54+(6*I)*303^(1/2))^(1/3)-4/(-54+(6*I)*303^(1/2))^(1/3)), -(1/12)*(-54+(6*I)*303^(1/2))^(1/3)-2/(-54+(6*I)*303^(1/2))^(1/3)-((1/2)*I)*3^(1/2)*((1/6)*(-54+(6*I)*303^(1/2))^(1/3)-4/(-54+(6*I)*303^(1/2))^(1/3))]

[(2/3)*6^(1/2)*sin((1/3)*arctan((1/9)*303^(1/2))+(1/6)*Pi), -(1/3)*3^(1/2)*2^(1/2)*(3^(1/2)*cos((1/3)*arctan((1/9)*303^(1/2))+(1/6)*Pi)+sin((1/3)*arctan((1/9)*303^(1/2))+(1/6)*Pi)), (1/3)*3^(1/2)*2^(1/2)*(3^(1/2)*cos((1/3)*arctan((1/9)*303^(1/2))+(1/6)*Pi)-sin((1/3)*arctan((1/9)*303^(1/2))+(1/6)*Pi))]

>	# POI := Points Of Interest (I use float values for the sake of simplicity) # Of course if you prefer using exact values replace the command below by # POI := := map(t -> [1, t], Yn); POI := map(t -> [1, t], Yn);

>	# Formal expression of the gtadient at any (x, y) point Gradient := unapply(lhs~(diff~(eq, [x, y])), (x, y)); # Gradients at POI gradients := map(t -> Gradient(op(t)), POI);

>	# Tangent slopes at the POI slopes := map(t -> -t[1]/t[2], gradients)

>	display( p1, p2, pointplot(POI, symbol=circle, symbolsize=20, color=black), seq( plot(slopes[i]*(x-POI[i][1])+POI[i][2], x=0..2, color=black, linestyle=3) , i=1..3 ) )

Download Proposal_2022.2.mw

A slight improvement where unitary gradients are also displayed : Proposal_with_gradients.mw

Analytical study of Eq roots...

Answered: mmcdara 7037

May 06 2024

3 3

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.

Does that suit you?...

Answered: mmcdara 7037

May 05 2024

1 2

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)

Basically,...

Answered: mmcdara 7037

April 22 2024

1 3

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.

Algebric generation of Pythagorean tripl...

Answered: mmcdara 7037

April 22 2024

0 1

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

Could you explain in a few words what th...

Answered: mmcdara 7037

April 18 2024

3 3

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))

(1)

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

(2)

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

(3)

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

(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))

(6)

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

(7)

>	_Int := evalf(aux_1*aux_2)

(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)

(9)

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

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

(10)

>

Download TruncatedNormal_2.mw

A first step...

Answered: mmcdara 7037

April 18 2024

2 0

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]])

(1)

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

(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}

(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 ):

(4)

>	# How the faces are connected in the pattern of the prism? ST := SpanningTree(G_Faces); DrawGraph(ST)

>

# 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)

>

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)

(5)

Download A_first_sketch_for_a_code.mw

A step-by-step explanation on how to do ...

Answered: mmcdara 7037

April 18 2024

1 0

>

restart

>	expr := 3G(`Δγ`*H - `σy`(`Δγ`) + q)/(-q + `σy`(`Δγ`))^2;

3*G*(`Δγ`*H-`σy`(`Δγ`)+q)/(-q+`σy`(`Δγ`))^2

(1)

>	indets(expr, function)[]; rw := % = Z+q; rw:= isolate(op(1, denom(expr))=Z, `σy`(`Δγ`))

(2)

>	algsubs(rw, expr)

3*G*(H*`Δγ`-Z)/Z^2

(3)

>	expand(%);

3*G*`Δγ`*H/Z^2-3*G/Z

(4)

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

3*G*`Δγ`*H/(-q+`σy`(`Δγ`))^2-3*G/(-q+`σy`(`Δγ`))

(5)

>	sort~(res, q, ascending)

3*G*`Δγ`*H/(`σy`(`Δγ`)-q)^2-3*G/(`σy`(`Δγ`)-q)

(6)

Download step_by_step.mw

@acer and I intertwined: this is b...

Answered: mmcdara 7037

April 17 2024

1 0

@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 -> Msum(Heaviside(t - 24k/N)exp((-t + 24k/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;

(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

If I understood correctly, here is an an...

Answered: mmcdara 7037

April 16 2024

3 3

At the end you get three different sets of solutions:

>

restart

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

(1)

>

>	Eq1 := -(alpha^2 + beta)u(xi) + k^2diff(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)

>

>	D1 := diff(f(xi), xi) = (a^(-f(xi))p + q + ra^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):

[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 = (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 = (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):

${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

Cartesian and cylindrical coordinates...

Answered: mmcdara 7037

April 15 2024

1 2

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)

(1)

PAPPUS' THEOREM

>	crossection_radius := 2: crossection_area := Picrossection_radius^2: Rotation_radius := 4: Torus_volume := crossection_area(2PiRotation_radius)

(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)

(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])

(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])

(8)

>	# Volume comprised between Ro(z) and Ri(z) Volume := OuterVolume - InnerVolume; identify(%);

(9)

Download torus_mmcdara.mw

If I understood correctly......

Answered: mmcdara 7037

April 07 2024

1 4

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

>

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]);

(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]) )

(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

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

Ask a Question

Create a Post

7037 Reputation

18 Badges

MaplePrimes Activity

These are answers submitted by mmcdara

First fix your errors to have a well po...

@Prakash J Look at this file D...

Bad use of D[2]...

A proposal...

Analytical study of Eq roots...

Does that suit you?...

Basically,...

Algebric generation of Pythagorean tripl...

Could you explain in a few words what th...

A first step...

A step-by-step explanation on how to do ...

@acer and I intertwined: this is b...

If I understood correctly, here is an an...

Cartesian and cylindrical coordinates...

If I understood correctly......

Save this setting as your default sorting preference?

Ask a Question

Create a Post

Generating PDF…

Save this setting as your default sorting preference?
Note: You can change your preference any time in your account settings
Don't show this again