8 years, 99 days

## Replace if Norm(err) <= tol th...

Be aware of undesired effects produced by break, look to this example

```restart
for n from 1 to 4 do
for i from 1 to 2 do
if n < 2 then
printf("   Within if structure %d  %d\n", n, i)
else
break
end if;
printf("Exit if structure %d  %d\n\n", n, i)
end do;
end do;

Within if structure 1  1
Exit if structure 1  1

Within if structure 1  2
Exit if structure 1  2
```

So just replace

```        if Norm(err) <= tol then
h := 0.1 * h * (c + 1) * (tol/Norm(err))^(0.2);
else
break
end if;
```

by

```        if Norm(err) <= tol then
h := 0.1 * h * (c + 1) * (tol/Norm(err))^(0.2);
end if;```

Are you sure Digits:=30 is necessary and that  Digits:=10 (default value) would npt ne enough?

Given the large number of points I would prefer one of these two graphics (but it is personal opinion):

```numerical_plot_y1 := plot(time_t, numerical_array_y1, style = point, symbol=point, color = blue, legend = ["TFIBF"]):
exact_plot_y1 := plot(time_t, exact_array_y1, style = point, symbol=point, color = red, legend = ["EXACT"]):
display({numerical_plot_y1, exact_plot_y1}, size=[800, 400]);

# or

numerical_plot_y1 := plot(time_t, numerical_array_y1, style = line, color = blue, legend = ["TFIBF"]):
exact_plot_y1 := plot(time_t[1..-2], exact_array_y1[1..-2], style = line, color = red, legend = ["EXACT"]):
display({numerical_plot_y1, exact_plot_y1}, size=[800, 400]);
```

## printlevel or trace?...

In Maple 2015 (the version I use) printlevel is an argument of kernelopts, not of interface.
Did you check that?
Note that you can also write printlevel := n directly.

Just a question: Why do you want to change the value of printlevel within your procedure?
If you want to display intermediate results I recommend you to use trace instead:

```restart
fac := proc(n::integer)
if n = 0 then
1
else
n * thisproc(n - 1);
end if;
end proc;
proc(n::integer)  ...  end;
fac(5)
120
trace(fac):
fac(5)

{--> enter fac, args = 5
{--> enter fac, args = 4
{--> enter fac, args = 3
{--> enter fac, args = 2
{--> enter fac, args = 1
{--> enter fac, args = 0
1
<-- exit fac (now in fac) = 1}
1
<-- exit fac (now in fac) = 1}
2
<-- exit fac (now in fac) = 2}
6
<-- exit fac (now in fac) = 6}
24
<-- exit fac (now in fac) = 24}
120
<-- exit fac (now at top level) = 120}
120

```

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

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

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

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.

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

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

 > restart
 > with(plots):
 > eq :=  x^3-4*x*y+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)

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

## Analytical study of Eq roots...

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?...

```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,...

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

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

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(2*x)*log(x)^2);
 (4)
 > Distrib := unapply( PDF(Normal(average, sqrt(var)), 1.1*x - 0.1), (x, average, var));
 (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(%);
 (10)
 >

## A first step...

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)

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

 > restart
 > expr := 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)
 > expand(%);
 (4)
 > res := eval(%, isolate(rw, Z)) assuming Z > 0
 (5)
 > sort~(res, q, ascending)
 (6)

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

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

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

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^2*diff(u(xi), xi \$ 2) + b*u(xi)^3 = 0
 (2)
 >
 >
 >
 >
 > D1 := 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);
 (4)
 > DEq1_0 := eval(Eq1, {D1, D2})
 (5)
 > # Setting Z = a^f(xi) DEq1_1 := collect(eval(expand(DEq1_0), a^f(xi) = Z), Z)
 (6)
 > # I don't understand why yhis command doesn't work ??? coeffs(DEq1_1, Z);
 > # 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):
 (7)
 > # After removing solutions which contain b=0 one gets three sols := convert(remove((x -> member(b=0, x)), sols), list): print~(sols):
 (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), %);
 (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):
 (10)