Advanced Search

sand15

1598 Reputation

16 Badges

11 years, 130 days
Contact sand15

MaplePrimes Activity


These are answers submitted by sand15

On one side xtest(t) is necessarily peri...

sand15 1598
April 02 2026
2 2



On one side xtest(t) is necessarily periodic with period 2*Pi.
On the other side it is obvious that x(t) is not 2*Pi periodic: It is therefore highly unlikely that xtest(t) is a good approximation of x(t).

Nevertheless...
A very simple way to find 'a' and 'b' is to do this:

T    := Array( [seq(t, t=0..60.1, 0.1)] ):
Tsol := dsolve({ics, Euler_eq}, numeric, output=T):

A := Tsol[2][1][.., 1]:
B := Tsol[2][1][.., 2]:

xtest := a0 + a1*cos(t) + b1*sin(t) + a2*cos(2*t) + b2*sin(2*t) + a3*cos(3*t) + b3*sin(3*t);

fit := Statistics:-Fit(xtest, A, B, t);

# compare x(t) and the best fit in the least squares sense

display(
  odeplot(Euler_sol, [t, x(t)], t=(min..max)(T), color=blue, legend="sol num")
  ,
  plot(fit, t=(min..max)(T), color=red, legend="fit")
)

(
Another simple way is to note that the sin(n*t) and cos(n*t) functions being mutually orthogonal over 2*Pi,  a[n] and b[n] are given by  a[n] = int(x(t)*cos(n*t), t=0..2*Pi) (up to a normalizing constant) and b[n] = int(x(t)*sin(n*t), t=0..2*Pi) (still up to a normalizing constant).
)


I also suggest you to compute the DFT (Discrete Fourier Transform) of the (rescaled) signal x(t) to observe what its frequencies are:

N := numelems(B);
C := SignalProcessing:-DFT(B);
display( 
  <
    < 
      listplot( B, 'title' = "Signal" );
      listplot( map( Re, C[1..N/2] ), 'title' = "Frequency" );
      listplot( map( Im, C[1..N/2] ), 'title' = "Phase Shift" ) 
    >
  > 
);

You will clearly see the frequency and shift patterns do not look like patterns of a periodic function(at least over the 0..60 range I chosed arbitrarily).


To investigate graphically the quasi-periodicity (introduction) of the solution, I suggest you to do this

SOL := dsolve({ics, Euler_eq}, numeric, method=rkf45, range=0..1000):  # 1000 is an arbitrary choice

odeplot(
  SOL
  , [t, x(t)]
  , t=0..200  # an upper bound not larger than 1000
  , refine=2
  , color=blue
  , legend="sol num"
  , axis[1]=[tickmarks=piticks]
  , gridlines=true
  , size=[1200, 400]
)

Drawing then 

SOL := dsolve({ics, Euler_eq}, numeric, method=rkf45, range=0..10000, maxfun=0):

odeplot(
  SOL
  , [x(t), diff(x(t), t)]
  , t=0..5000
  , refine=2
  , color=blue
  , legend="sol num"
  , axis[1]=[tickmarks=piticks]
  , gridlines=true
  , size=[1200, 400]
);

reveals a pattern close to the behaviour of a quasi-perriodic system (which could maybe be chaotic for other values of 'u' and 'v'). See full view and a zoom below:


Here is a zoom which depicts quasi-periodicity where the pattern  {x(t), t=0..21.811}  is translated by tau=63 and tau=125.7:



The long range fit of x(t) by xtest(t)  looks like this
The constant component (a[0]) should likely converge to 0 when the upper bound of the range (below t=0..1000) used to assess xtest(t) tends to infinity.
Because the harmonic components you consider both oscillate more rapidly (periods about 6, 3, 2) than the apparent periodicity (which seems to be arround 63, see above), those components must necessarily have a coefficient close to 0 too).


Indeed 


Finally x(t) is not periodic and its spectrum (see the DFT above) is continuous (not made of a countable number of rays like any periodic function spectrum).
I think quasi-periodicity or chaoticity of x(t) could be proved mathematically but I'm not qualified to do that (see for instance Fayad & Krikorian)

Using energy conservation...

sand15 1598
March 29 2026
0 3

@nm 

(Note that  _C1 disappears in relation (7) because x → ∞).

> 

restart
> 

ode_0 := diff(y(x), x$2) = sin(y(x))

diff(diff(y(x), x), x) = sin(y(x))

 (1)
> 

ode_1 := ode_0 * diff(y(x), x);

(diff(y(x), x))*(diff(diff(y(x), x), x)) = (diff(y(x), x))*sin(y(x))

 (2)
> 

ode_2 := map(int, ode_1, x);

(1/2)*(diff(y(x), x))^2 = -cos(y(x))

 (3)
> 

odes := map(allvalues, isolate(ode_2, diff(y(x), x)));

diff(y(x), x) = ((-2*cos(y(x)))^(1/2), -(-2*cos(y(x)))^(1/2))

 (4)


I investigate only the  sqrt(-2cos(y(x))) case, feel free to investigate the other one

> 

solimp := dsolve(lhs(odes) = op(1, [rhs(odes)]))

x-Intat(1/(-2*cos(_a))^(1/2), _a = y(x))+_C1 = 0

 (5)
> 

solexp := dsolve(lhs(odes) = op(1, [rhs(odes)]), explicit)

y(x) = RootOf(x-Intat(1/(-2*cos(_a))^(1/2), _a = _Z)+_C1)

 (6)
> 

limitsol := map(limit, solexp, x=+infinity)

limit(y(x), x = infinity) = limit(RootOf(x-Intat(1/(-2*cos(_a))^(1/2), _a = _Z)+_C1), x = infinity)

 (7)
> 

IntegralTerm := remove(has, select(has, indets(solexp, function), Intat), RootOf)[]

IntegralTerm := Intat(1/sqrt(-2*cos(_a)), _a = _Z)

 (8)
> 

integralTerm := value(IntegralTerm)

2^(1/2)*(4*cos((1/2)*_Z)^2-2)^(1/2)*InverseJacobiAM((1/2)*_Z, 2^(1/2))/(-4*cos((1/2)*_Z)^2+2)^(1/2)

 (9)
> 

limitsol := eval(limitsol, IntegralTerm = integralTerm)

limit(y(x), x = infinity) = RootOf(4*cos((1/2)*_Z)^2-2)

 (10)
> 

map(allvalues, limitsol);

limit(y(x), x = infinity) = (-Pi*_B1+(1/2)*Pi+4*Pi*_Z1, -3*Pi*_B2+(3/2)*Pi+4*Pi*_Z2)

 (11)
> 

BZ := convert(select(hasassumptions, indets((11), name)), list)

[_B1, _B2, _Z1, _Z2]

 (12)
> 

about~(BZ);

_B1:

  is assumed to be: OrProp(0,1)

_B2:
  is assumed to be: OrProp(0,1)

Originally _Z1, renamed _Z1~:
  is assumed to be: integer

Originally _Z2, renamed _Z2~:
  is assumed to be: integer

  

[]

 (13)
> 

eval((11), BZ =~1)

limit(y(x), x = infinity) = ((7/2)*Pi, (5/2)*Pi)

 (14)
> 

for b1 in [0, 1] do
  for b2 in [0, 1] do
    for z1 from -1 to 1 do
      for z2 from -1 to 1 do
        print(eval([b1, b2, z1, z2, (11)], BZ =~[b1, b2, z1, z2]))
      end do:
    end do:
  end do:
end do:

[0, 0, -1, -1, limit(y(x), x = infinity) = (-(7/2)*Pi, -(5/2)*Pi)]

  

[0, 0, -1, 0, limit(y(x), x = infinity) = (-(7/2)*Pi, (3/2)*Pi)]

  

[0, 0, -1, 1, limit(y(x), x = infinity) = (-(7/2)*Pi, (11/2)*Pi)]

  

[0, 0, 0, -1, limit(y(x), x = infinity) = ((1/2)*Pi, -(5/2)*Pi)]

  

[0, 0, 0, 0, limit(y(x), x = infinity) = ((1/2)*Pi, (3/2)*Pi)]

  

[0, 0, 0, 1, limit(y(x), x = infinity) = ((1/2)*Pi, (11/2)*Pi)]

  

[0, 0, 1, -1, limit(y(x), x = infinity) = ((9/2)*Pi, -(5/2)*Pi)]

  

[0, 0, 1, 0, limit(y(x), x = infinity) = ((9/2)*Pi, (3/2)*Pi)]

  

[0, 0, 1, 1, limit(y(x), x = infinity) = ((9/2)*Pi, (11/2)*Pi)]

  

[0, 1, -1, -1, limit(y(x), x = infinity) = (-(7/2)*Pi, -(11/2)*Pi)]

  

[0, 1, -1, 0, limit(y(x), x = infinity) = (-(7/2)*Pi, -(3/2)*Pi)]

  

[0, 1, -1, 1, limit(y(x), x = infinity) = (-(7/2)*Pi, (5/2)*Pi)]

  

[0, 1, 0, -1, limit(y(x), x = infinity) = ((1/2)*Pi, -(11/2)*Pi)]

  

[0, 1, 0, 0, limit(y(x), x = infinity) = ((1/2)*Pi, -(3/2)*Pi)]

  

[0, 1, 0, 1, limit(y(x), x = infinity) = ((1/2)*Pi, (5/2)*Pi)]

  

[0, 1, 1, -1, limit(y(x), x = infinity) = ((9/2)*Pi, -(11/2)*Pi)]

  

[0, 1, 1, 0, limit(y(x), x = infinity) = ((9/2)*Pi, -(3/2)*Pi)]

  

[0, 1, 1, 1, limit(y(x), x = infinity) = ((9/2)*Pi, (5/2)*Pi)]

  

[1, 0, -1, -1, limit(y(x), x = infinity) = (-(9/2)*Pi, -(5/2)*Pi)]

  

[1, 0, -1, 0, limit(y(x), x = infinity) = (-(9/2)*Pi, (3/2)*Pi)]

  

[1, 0, -1, 1, limit(y(x), x = infinity) = (-(9/2)*Pi, (11/2)*Pi)]

  

[1, 0, 0, -1, limit(y(x), x = infinity) = (-(1/2)*Pi, -(5/2)*Pi)]

  

[1, 0, 0, 0, limit(y(x), x = infinity) = (-(1/2)*Pi, (3/2)*Pi)]

  

[1, 0, 0, 1, limit(y(x), x = infinity) = (-(1/2)*Pi, (11/2)*Pi)]

  

[1, 0, 1, -1, limit(y(x), x = infinity) = ((7/2)*Pi, -(5/2)*Pi)]

  

[1, 0, 1, 0, limit(y(x), x = infinity) = ((7/2)*Pi, (3/2)*Pi)]

  

[1, 0, 1, 1, limit(y(x), x = infinity) = ((7/2)*Pi, (11/2)*Pi)]

  

[1, 1, -1, -1, limit(y(x), x = infinity) = (-(9/2)*Pi, -(11/2)*Pi)]

  

[1, 1, -1, 0, limit(y(x), x = infinity) = (-(9/2)*Pi, -(3/2)*Pi)]

  

[1, 1, -1, 1, limit(y(x), x = infinity) = (-(9/2)*Pi, (5/2)*Pi)]

  

[1, 1, 0, -1, limit(y(x), x = infinity) = (-(1/2)*Pi, -(11/2)*Pi)]

  

[1, 1, 0, 0, limit(y(x), x = infinity) = (-(1/2)*Pi, -(3/2)*Pi)]

  

[1, 1, 0, 1, limit(y(x), x = infinity) = (-(1/2)*Pi, (5/2)*Pi)]

  

[1, 1, 1, -1, limit(y(x), x = infinity) = ((7/2)*Pi, -(11/2)*Pi)]

  

[1, 1, 1, 0, limit(y(x), x = infinity) = ((7/2)*Pi, -(3/2)*Pi)]

  

[1, 1, 1, 1, limit(y(x), x = infinity) = ((7/2)*Pi, (5/2)*Pi)]

 (15)
 

 

Download Not_Only_Pi.mw

Hard to get a pretty animation...

sand15 1598
March 28 2026
2 1

Assuming (see @vv) the conditions for BSC to be isocele hold, the solution can be obtained elegantly (IMO) using package geom3d.

Nevertheless drawing a pretty animation is quite difficult.
So I propose you this bac23IllustrationEspace_sand15.mw worksheet and let you customize it as you want:

  • equation of plane P,
  • locations of B and C,
  • objects to draw/plot
  • colors and transparencies
  • point of view, 
  • add of textplot3d,
  • use of spacecurve instead of plot/PLOT
  • ... and so on.

Descriptions of the colord objects is given in the attached worksheet.

Typo... and formal expressions of PDF(d)...

sand15 1598
March 20 2026
1 1

You wrote d = ... instead of d := ...

This should answer your question.
Example of use:

Histogram(Sample(d, 10^4))



If you want to compute formally some statistics of d, for instance its PDF or CDF that is another story.

If you directly ask Maple for the PDF/CDF of sqrt(DX2+DY2) you won't get the expected result (I can't with my 2015 version).

To get this PDF (CDF) you have to proceed step by step in order to help Maple:

  • Build the PDF (or CDF) of DX2
     
  • Build the CDF (or PDF) of DY2
     
  • Compute the CDF of DX2+DY2, either from PDF(DX2) and CDF(DY2), or from CDF(DX2) and PDF(DY2) (much simpler than computing directly the PDF of DX2+DY2 from rhe convolution product of  PDF(DX2) by CDF(DY2)).
    Let K(t) the CDF of DX2+DY2.
     
  • Define a random variable Z2 this way 
    Z2 := RandomVariable(Distribution(CDF = (t -> K(t))));
  • Finally ask Maple what is the PDF of sqrt(Z2) (Maple 2015 does the job).
    Maybe your Maple version won't do it, if it is the case let me know.


You will get then a (very) lengthy expression for PDF(Z2) that you can compare to a direct sample of sqrt(DX2+DY2).
 

Here is the full code (final result at the end of the worksheet) TriangleEuclidean_sand15.mw

and a plot comparing PDF(Z2) to the histogram of a sample of sqrt(DX2+DY2):


The expression of PDF( sqrt(DX2+DY2) ) is that lengthy that it is almost unreadable (I tried a few ad hoc simplifications, mostly by hand, but the expression of PDF(Z2) remains complex, see the attached file).

To help you make an idea of how complex it is, here is its optimized Python code:

t3 = t ** 2
t4 = 4 * t3
t5 = t4 - 3
t6 = math.sqrt(t5)
t9 = t4 - 1
t10 = math.sqrt(t9)
t11 = 0.1e1 / t10
t13 = t3 < 1
t14 = t3 == 0.1e1 / 0.4e1
t15 = math.sqrt(3)
t17 = t3 + 0.1e1 / 0.3e1
t19 = 2 * t3
t20 = t19 - 1
t21 = 1 / t3
t24 = math.asin(t21 * t20 / 2)
t30 = t3 ** 2
t31 = math.sqrt(t3)
t32 = t31 * t30
t34 = t31 * t3
t36 = t30 * t3
t37 = t31 * t36
t39 = t3 * math.pi
t44 = t30 ** 2
t58 = 0.32e2 / 0.243e3 * (-27 * t24 * t17 * t3 * t10 * t15 + t10 * (t32 * (-0.144e3 / 0.5e1 * t15 + 0.152e3 / 0.5e1) + 72 * t34 - 0.192e3 / 0.35e2 * t37 + t15 * (-0.27e2 / 0.2e1 * t39 + 0.45e2 / 0.32e2 * math.pi + 0.9e1 / 0.10e2) + t44 + 3 * t36 - 0.159e3 / 0.2e1 * t30 - 0.34e2 / 0.5e1 * t3 - 0.29471e5 / 0.8960e4) + (0.402e3 / 0.5e1 * t30 + 0.384e3 / 0.5e1 * t36 - 0.231e3 / 0.10e2 * t3 - 0.9e1 / 0.20e2) * t15) * t11 if t14 else 0
t59 = t58 if t13 else 0
t60 = t6 * t59
t62 = 12 * t3 - 9
t63 = math.sqrt(t62)
t64 = t63 * t10
t67 = t3 - 0.2e1 / 0.3e1
t68 = 0 if t13 else 1
t70 = 0 if t3 < 0.3e1 / 0.4e1 else 1
t71 = t68 - t70
t72 = t71 * t67
t73 = t6 / 2
t74 = -t73 < 0
t75 = 0 if t74 else 1
t77 = t73 < 0
t78 = 0 if t77 else 1
t80 = t3 + 0.4e1 / 0.3e1
t82 = -t80 * t68 + t75 * t72 - t78 * t72
t83 = -0.3e1 / 0.4e1 + t3
t84 = t83 * t82
t85 = t10 * t3
t86 = t19 - 3
t89 = math.asin(t21 * t86 / 2)
t90 = t89 * t85
t94 = 0 if t3 < 0.1e1 / 0.4e1 else 1
t95 = t68 - t94
t96 = t15 * t95
t97 = t6 * t96
t98 = t24 * t17
t99 = t98 * t85
t105 = t15 * math.pi
t106 = t105 * t3 * t67
t108 = 9 * t36
t111 = 0.48e2 / 0.5e1 * t32 * t15 - 0.9e1 / 0.4e1 * t106 + t44 - t108 - 0.81e2 / 0.8e1 * t30 + 0.27e2 / 0.16e2 * t3 - 0.243e3 / 0.1280e4
t118 = t6 * t111 - 0.128e3 / 0.35e2 * t83 * (t36 - 0.26e2 / 0.3e1 * t30 + 0.433e3 / 0.128e3 * t3 + 0.333e3 / 0.256e3)
t119 = t118 * t71
t120 = t75 * t10
t122 = t78 * t10
t124 = t15 * t68
t125 = t6 * t124
t126 = t3 * t80
t127 = t3 - 2
t129 = math.asin(t21 * t127)
t131 = t129 * t10 * t126
t134 = t3 + 3
t135 = t3 - 1
t136 = math.sqrt(t135)
t137 = t136 * t135
t138 = t137 * t134
t141 = t135 ** 2
t142 = t136 * t141
t157 = 0.80e2 / 0.3e1 * t32 + 12 * t34 - 0.128e3 / 0.35e2 * t37 + t15 * (-0.3e1 / 0.2e1 * t39 - 0.9e1 / 0.2e1 * t30 * math.pi) - 0.53e2 / 0.2e1 * t30 + 0.4e1 / 0.3e1 * t44 - 8 * t36 + 0.17e2 / 0.448e3 - 0.34e2 / 0.15e2 * t3
t159 = 0.48e2 / 0.5e1 * t15
t160 = -0.152e3 / 0.15e2 + t159
t162 = t68 * t160 - t159 - 0.248e3 / 0.15e2
t165 = -24 * t68 + 12
t167 = t68 + 1
t169 = t3 - 0.10e2 / 0.3e1
t170 = t136 * t169
t178 = math.pi * t15 * t3
t179 = 0.9e1 / 0.2e1 * t178
t180 = -0.2e1 / 0.3e1 * t44 + 0.439e3 / 0.8e1 * t30 + 0.7493e4 / 0.240e3 * t3 - 2 * t36 - 0.2593e4 / 0.1792e4 + t179
t183 = -8 * t138 * t124 + 0.24e2 / 0.5e1 * t142 * t124 + t94 * t157 + t32 * t162 + t34 * t165 + 0.64e2 / 0.35e2 * t37 * t167 + 9 * t170 * t124 + t68 * t180 + 0.9e1 / 0.2e1 * t106 - t44 + t108
t185 = t83 * t68
t188 = t36 - 0.31e2 / 0.6e1 * t30 - 0.1919e4 / 0.128e3 * t3 - 0.423e3 / 0.256e3
t191 = t6 * t183 + 0.128e3 / 0.35e2 * t188 * t185
t194 = t30 + 0.83e2 / 0.64e2 * t3 + 0.3e1 / 0.128e3
t195 = t194 * t95
t197 = (-0.1e1 / 0.4e1 + t3) * t15
t198 = t6 * t197
t201 = 9 * t99 * t97 + t120 * t119 - t122 * t119 + 0.9e1 / 0.2e1 * t131 * t125 + t10 * t191 - 0.128e3 / 0.5e1 * t198 * t195
t203 = -0.81e2 / 0.32e2 * t64 * t60 - 54 * t90 * t84 + t63 * t201
t204 = 0.1e1 / t63
t205 = t204 * t203
t208 = 0.1e1 / t6
t214 = t11 * t208
t215 = undefined if t14 else 0
t230 = undefined if t3 == 1 else 0
t232 = undefined if t3 == 0.3e1 / 0.4e1 else 0
t233 = t230 - t232
t234 = t233 * t67
t236 = 0 if t74 else 0
t240 = 0 if t77 else 0
t248 = t89 * t10
t253 = t11 * t3
t258 = 1 / t30
t264 = math.sqrt(t258 * t62)
t269 = t230 - t215
t280 = t15 * t230
t342 = t15 * math.pi * t67
t344 = 4 * t36
t345 = 27 * t30
t346 = 0.64e2 / 0.35e2 * t37 * t230 + 0.32e2 / 0.5e1 * t32 * t167 + 9 * t170 * t280 + 9 * t136 * t124 + 0.9e1 / 0.2e1 / t136 * t169 * t124 + t68 * (-0.8e1 / 0.3e1 * t36 + 0.439e3 / 0.4e1 * t3 + 0.7493e4 / 0.240e3 - 6 * t30 + 0.9e1 / 0.2e1 * t105) + t230 * t180 + t179 + 0.9e1 / 0.2e1 * t342 - t344 + t345
t356 = 3 * t30
t365 = t118 * t233
t384 = (t6 * (24 * t34 * t15 - 0.9e1 / 0.4e1 * t178 - 0.9e1 / 0.4e1 * t342 + t344 - t345 - 0.81e2 / 0.4e1 * t3 + 0.27e2 / 0.16e2) + 2 * t208 * t111 - 0.128e3 / 0.35e2 * t83 * (t356 - 0.52e2 / 0.3e1 * t3 + 0.433e3 / 0.128e3) - 0.128e3 / 0.35e2 * t36 + 0.3328e4 / 0.105e3 * t30 - 0.433e3 / 0.35e2 * t3 - 0.333e3 / 0.70e2) * t71
t392 = 9 * t99 * t6 * t15 * t269 + 18 * t99 * t208 * t96 + 18 * t98 * t253 * t97 + 0.9e1 / 0.2e1 * t131 * t6 * t280 + 9 * t131 * t208 * t124 + 9 * t129 * t11 * t126 * t125 + t10 * (t6 * (-8 * t138 * t280 + 4 * t137 * t124 - 12 * t136 * t134 * t124 + 0.24e2 / 0.5e1 * t142 * t280 + t94 * (0.200e3 / 0.3e1 * t34 + 18 * t31 - 0.64e2 / 0.5e1 * t32 + t15 * (-0.3e1 / 0.2e1 * math.pi - 9 * t39) - 53 * t3 + 0.16e2 / 0.3e1 * t36 - 24 * t30 - 0.34e2 / 0.15e2) + t215 * t157 + t32 * t230 * t160 + 0.5e1 / 0.2e1 * t34 * t162 - 24 * t34 * t230 + 0.3e1 / 0.2e1 * t31 * t165 + t346) + 2 * t208 * t183 + 0.128e3 / 0.35e2 * t188 * t83 * t230 + 0.128e3 / 0.35e2 * t188 * t68 + 0.128e3 / 0.35e2 * (t356 - 0.31e2 / 0.3e1 * t3 - 0.1919e4 / 0.128e3) * t185) + 2 * t11 * t191 + t120 * t365 + t120 * t384 + t236 * t10 * t119 + 2 * t75 * t11 * t119 - t122 * t365
t407 = math.sqrt(4) * math.sqrt(t258 * t135)
t415 = t17 * t10
t420 = math.sqrt(t258 * t9)
t449 = -t122 * t384 - t240 * t10 * t119 - 2 * t78 * t11 * t119 - 0.128e3 / 0.5e1 * t6 * t15 * t195 + 0.9e1 / 0.2e1 / t407 * (-t258 * t127 + t21) * t85 * t80 * t6 * t124 + 18 / t420 * (t21 - t258 * t20 / 2) * t415 * t3 * t6 * t96 + 9 * t24 * t415 * t97 + 9 * t24 * t85 * t97 + 0.9e1 / 0.2e1 * t129 * t85 * t125 + 0.9e1 / 0.2e1 * t129 * t10 * t80 * t125 - 0.128e3 / 0.5e1 * t198 * t194 * t269 - 0.128e3 / 0.5e1 * t198 * (t19 + 0.83e2 / 0.64e2) * t95 - 0.256e3 / 0.5e1 * t208 * t197 * t195
t454 = -0.81e2 / 0.32e2 * t64 * t6 * t215 - 0.81e2 / 0.16e2 * t64 * t208 * t59 - 0.81e2 / 0.16e2 * t63 * t11 * t60 - 0.243e3 / 0.16e2 * t204 * t10 * t60 - 54 * t90 * t83 * (-t80 * t230 + t75 * t234 - t78 * t234 + t236 * t72 - t240 * t72 + t75 * t71 - t78 * t71 - t68) - 54 * t248 * t3 * t82 - 54 * t248 * t84 - 108 * t89 * t253 * t84 - 108 / t264 * (t21 - t258 * t86 / 2) * t10 * t3 * t84 + t63 * (t392 + t449) + 6 * t204 * t201
t466 = 0 if t <= 0 else 2 * (0.64e2 / 0.81e2 * t205 * t11 / t6 / t5 + 0.64e2 / 0.81e2 * t205 / t10 / t9 * t208 - 0.32e2 / 0.81e2 * t204 * t454 * t214 + 0.64e2 / 0.27e2 / t63 / t62 * t203 * t214) * t if 0 < t else 0


Or enlarge this image

A Maple2015 solution...

sand15 1598
March 16 2026
2 0

Customize as you want

NULL

 

NULL

``

> 

restart
> 

with(Statistics):
> 

X := RandomVariable(Normal(0, 1)):
> 

f := unapply(PDF(X, x), x):
> 

g := unapply(CDF(X, x), x):
> 

``
> 

h := proc(t)
if t::numeric then
  plots:-display(
    plot(f(x),  x = -5 .. t, gridlines = true, filled=true, color="Chartreuse")
    ,
    plot(g(x), x = -5..t, gridlines=false, color=blue)
    ,
    plot([[t, 0], [t, g(t)]], linestyle=3, color="DarkGray")
    , title=typeset(Int('f'(x), x=-5..evalf[2](t)) = evalf[4](g(t)))
  )
else
  'procname(_passed)'
end if:
end proc:

plots:-animate(h, [t], t=-5.0 .. 5.0, paraminfo=false)

 
> 

 

``

``

NULL
 

 

NULL

``

 

Download Q_CDF_PDF_sand15.mw


A variant using dualaxisplot: (look how to suppress the info about the parameter value [Maple 2015]) 

h := proc(t)
if t::numeric then
  plots:-dualaxisplot(
    plot(f(x),  x = -5 .. t, gridlines = true, filled=true, color="Chartreuse", legend=:-PDF('X'))
    ,
    plots:-display(
      plot(g(x), x = -5..t, gridlines=false, color=blue, legend=:-CDF('X'))
      ,
      plot([[t, 0], [t, g(t)]], linestyle=3, color="DarkGray")
      , axis[2]=[color=blue]
    )
    , title=typeset(Int('f'(x), x=-5..evalf[2](t)) = evalf[4](g(t)))
  )
else
  'procname(_passed)'
end if:
end proc:

plots:-animate(h, [t], t=-5.0 .. 5.0, paraminfo=false)

A solution...

sand15 1598
March 08 2026
1 0


With restrictions  y > 1 and _C1::real one gets this solution
Asolution.mw


Given we assumed y > 1, the solution is A_solution_gen[1].

Simply use the ad hoc outputs in Polynom...

sand15 1598
March 06 2026
0 4

I don't think that using the Horner representation is judicious.
What we generally want to know are the coefficients (that is the effects, and often the statistical significances) of xn terms, an information the Horner factorization does not provide directly, which makes it useless.

Note: be careful in using statistics returned by the output=solutionmodule option in the case you use PolynomialFit or ExponentialFit because some are wrong (I believe I posted something about that, either under this login or under the login @mmcdara, I don't remembe exactly).

Nevertheless "leastsquaresfunction" and "residuals" output, and a few more, are correct.

> 

restart;
> 

X   := [seq(0..7,0.5)]:
phi := x -> sin(x):
Y   := phi~(X):


Only necessary outputs

> 

fit, residuals := op(Statistics:-PolynomialFit(10, X, Y, x, output=["leastsquaresfunction", "residuals"])):

Predictor := unapply(fit, x):

plot(
  [phi, Predictor](x)
  , x=(min..max)(X)
  , color=[blue, red]
  , transparency=[0, 0.7]
  , thickness=[2, 11]
  , legend=[typeset('phi'(x)), typeset('Predictor'(x))]
);

plot(
  [ seq([X[i], residuals[i]], i=1..numelems(X)) ]
)

 

 


You might be interested to know a lot of other outpurts are avaliable

> 

# Yo get all of them:

P := Statistics:-PolynomialFit(10, X, Y, x, output=solutionmodule):
> 

# To see them

P:-Results():
> 

# Extract residuals alone

residuals := P:-Results("residuals"):

dataplot(residuals);

 
> 

Predictor := unapply(P:-Results("leastsquaresfunction"), x):

plot(
  [phi, Predictor](x)
  , x=(min..max)(X)
  , color=[blue, red]
  , transparency=[0, 0.7]
  , thickness=[2, 11]
  , legend=[typeset('phi'(x)), typeset('Predictor'(x))]
)

 
 

 

Download regr_fun_sand15.mw

Why do you say there is an error?...

sand15 1598
March 02 2026
1 5

I suggest you not to load simultanously Student[VectorCalculus] and LinearAlgebra to avoid conclicts (there are some in Maple 2015)
Once done (Maple 2015) I can't see any error.

Do you mean instead that the result Maple produces does not suit you?
If it is so maybe you should verify twice what you did by hand...

Question_1_sand15.mw

Fomal integration seems impossible...

sand15 1598
March 02 2026
2 6

EDITED

Because the formal integration of f,  not even to speak of the convolution product f ⊗ f, does not seem possible, I propose you two alternatives:

  1. A pointwise approximation of f followed by a discrete approximation (in fact a Riemann sum) of the integral which represents f ⊗ f
  2. A spline reconstruction of f  followed by an exact calculation of  f ⊗ f

Both approaches give equivalent results provided enough points are considered in both of them (I used 41 points in approach (1) and 21 in approach (2) but those numbers were chosen arbitrarily and it is likely that already good results could be obtained with smaller numbers).

Details are given in the attached worksheet, the animation below is the ice on the cake (if you like iced cakes of course).

Here is an updated version of my initial worksheet which contains a numerical estimation of  f ⊗ f  using evalf/Int (cpu times seems prohibitive).

lternatives_to_formal_integration_updated.mw

echoexpression=false...

sand15 1598
February 28 2026
2 0

add option echoexpression=false

Maple 2015 help:
help(Explore) > section "See Also" > Explore example worksheet > Exploring a mathematical expression > example 2

One way...

sand15 1598
February 26 2026
4 4

Your question can be highly simplified if you ask if the arguments of the log functions in G0 and G01 are equal.

NULL

> 

restart
> 

with(plots):
> 

H0 := -S1^2*eta1-S2^2*eta2-S1*gamma1-S2*gamma2;

-S1^2*eta1-S2^2*eta2-S1*gamma1-S2*gamma2

 (1)
> 

NULL
> 

Z0 := exp(-beta*H0);

exp(-beta*(-S1^2*eta1-S2^2*eta2-S1*gamma1-S2*gamma2))

 (2)
> 

Z0 := add(Z0, S1 = [-2, -1, 0, 1, 2]):
> 

# Your problem will be simpler to solve if you consider Z0 instead of G0
> 

> 

Z01 := (2*exp(4*beta*eta1)*cosh(2*beta*gamma1)+2*exp(beta*eta1)*cosh(beta*gamma1)+1)^((1/2)*N)*(2*exp(4*beta*eta2)*cosh(2*beta*gamma2)+2*exp(beta*eta2)*cosh(beta*gamma2)+1)^((1/2)*N)

(2*exp(4*beta*eta1)*cosh(2*beta*gamma1)+2*exp(beta*eta1)*cosh(beta*gamma1)+1)^((1/2)*N)*(2*exp(4*beta*eta2)*cosh(2*beta*gamma2)+2*exp(beta*eta2)*cosh(beta*gamma2)+1)^((1/2)*N)

 (3)
> 

# convert(Z01, exp) replaces hyperbolic trigs by exponentials

T := simplify(convert(Z01, exp));

(exp(2*beta*(2*eta1+gamma1))+exp(beta*(eta1+gamma1))+exp(2*beta*(2*eta1-gamma1))+exp(beta*(eta1-gamma1))+1)^((1/2)*N)*(exp(2*beta*(2*eta2+gamma2))+exp(beta*(eta2+gamma2))+exp(2*beta*(2*eta2-gamma2))+exp(beta*(eta2-gamma2))+1)^((1/2)*N)

 (4)
> 

# replace T = AN/2 * BN/2 by (A*B)N/2
# and discard the exponent already in common with G0

T0 := op([1, 1], T)* op([2, 1], T)

(exp(2*beta*(2*eta1+gamma1))+exp(beta*(eta1+gamma1))+exp(2*beta*(2*eta1-gamma1))+exp(beta*(eta1-gamma1))+1)*(exp(2*beta*(2*eta2+gamma2))+exp(beta*(eta2+gamma2))+exp(2*beta*(2*eta2-gamma2))+exp(beta*(eta2-gamma2))+1)

 (5)
> 

# the question becomes: Is T0 equal to Z0?

simplify(T0 - Z0)

0

 (6)
 

 

Download TESTE_MAPLE_sand15.mw

Branches...

sand15 1598
February 26 2026
2 1

You cannot simply use the commands solve(cot(y)=A, y) and solve(cot(y)=-A, y) to get the inverse function of cot function: arccot is made of an infinity of branches and the expression of cot(-1) depends on the branch you consider.

Using solve(cot(y)=A, y) and solve(cot(y)=-A, y) gives you the impression that cot(-1) is made of two branches, one for A > 0, the other for A < 0.
(BTW how can Maple distinguish these two cases without explicitely using assumptions 

solve(cot(y) =  A, y, useassumptions=true) assuming A > 0;
solve(cot(y) = -A, y, useassumptions=true) assuming A > 0;

?)

Use solve(cot(y)=A, y, allsolutions) instead to build cot(-1) on the branch of interest.

> 

restart
> 

kernelopts(version)

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

 (1)
> 

branches(arccot)

 
> 

# The inverse function of 'cot' depends on the branch you consider

`#msup(mo("cot"),mo("-1"))` := solve(cot(y)=A, y, allsolutions);

about(_Z1);

arccot(A)+Pi*_Z1

  

Originally _Z1, renamed _Z1~:
  is assumed to be: integer

  
> 

# Main branch

`#msup(mo("cot"),mo("-1")) ` = eval(`#msup(mo("cot"),mo("-1"))`, _Z1=0)

`#msup(mo("cot"),mo("-1")) ` = arccot(A)

 (2)
> 

# branch (-Pi, 0)

`#msup(mo("cot"),mo("-1")) ` = eval(`#msup(mo("cot"),mo("-1"))`, _Z1=-1)

`#msup(mo("cot"),mo("-1")) ` = arccot(A)-Pi

 (3)
> 

# branch (-2*Pi, -Pi)

`#msup(mo("cot"),mo("-1")) ` = eval(`#msup(mo("cot"),mo("-1"))`, _Z1=-2)

`#msup(mo("cot"),mo("-1")) ` = arccot(A)-2*Pi

 (4)
> 

# branch (Pi, 2*Pi)

`#msup(mo("cot"),mo("-1")) ` = eval(`#msup(mo("cot"),mo("-1"))`, _Z1=1)

`#msup(mo("cot"),mo("-1")) ` = arccot(A)+Pi

 (5)
 

 

Download arccot.mw

You say that other CAS answer differently? I advice you to ree you to read carefully how the inverse trig functions are (by default) defined and how to get them on any branch. 
Maybe Maple has its own strategy about what the "default" solve(cot(y)=A, y) must return,  but I doubt using Maple, MMA or SageMath give another answer if they are used cautiously.

Thanks...

sand15 1598
February 23 2026
2 8

@FDS 

A question: your second graph shows your material experiences some kind of strain-hardening and does not have an elastic regime.
As there is not linear part in the charge and discharge branches of the "stress-strain curve" (between quotes, see the end of the next paragraph), how do you define the slopes you are looking for?

By the way I'm quite surprised by the unusual concavity of the loading branches of the stress-strain curve: Shouldn't have an opposite concavity to those of the unloading branches?  What is the constitutive law which governs the behaviour of your material  (it looks like some kind of dilatant fluid)?
Would the horizonal axis represent the shear rate instead of the strain?

Please have a look to Stress_Strain_Curve.mw and tell me in what direction go.
(The reason I fit charge/discharge branches with polynomial models comes from my intuition you are working with a shear-thickening fluid, fluids whose behaviours are usually represented by powerlaw models... but maybe I am on the wrong track).

Based on MathML...

sand15 1598
February 15 2026
3 0

With "my" Maple 2015 (but this would work for more recent versions too), and using MathML:

> 

restart
> 

# MathML hexa code of symbol x

`#mo("&#x2243;")`

`#mo("&#x2243;")`

 (1)
> 

# Use this symbol as a binary operator.. version 1

`&#x2243;`(Pi, 3.14);  # not very pretty

`&#x2243;`(Pi, 3.14)

 (2)
> 

# Version 2;: define operator 'ae' (almost equal)
# see acer's comment
`print/ae` := proc(x,y)
   uses Typesetting;
   mrow(Typeset(x), mo("&InvisibleTimes;"),
        mo("&#x2243;",fontweight="bold"), mo("&InvisibleTimes;"),
        Typeset(y));
end proc:
> 

ae(Pi, 3.14);  # prettier

ae(Pi, 3.14)

 (3)
> 

`print/ra` := proc(x,y)
   uses Typesetting;
   mrow(Typeset(x), mo("&InvisibleTimes;"),
        mo("&#x2192;",fontweight="bold"), mo("&InvisibleTimes;"),
        Typeset(y));
end proc:

plot(cos(x), x=0..2, title = typeset(ra(cos(x)=0, ae(x, 1.6))), titlefont=[Times, bold, 16])

 
 

 

Download AlmostEqual.mw

to MacGyver...

sand15 1598
February 15 2026
2 0

Understandable, not "elegant" and far from being a one-liner command (unless you wrapped this intoa procedure... if course).

> 

ex:=(a+b*c)^2

(b*c+a)^2

 (1)
> 

shorten := table([_Inert_POWER="`^`", _Inert_SUM="`+`", _Inert_PROD="`*`", "("="[", ")"="]"]):
> 

inertex := ToInert(ex):
> 

shortstr := convert( eval(inertex, {_Inert_NAME = (u -> parse(u)), _Inert_INTPOS = (u -> u)}), string):
> 

shortstr := copy(shortstr):
for i in [indices(shorten, nolist)] do
  shortstr := StringTools:-SubstituteAll(shortstr, i, short[i])
end do:

cat("[", shortstr, "]")

"[`^`[`+`[`*`[b,c],a],2]]"

 (2)
 

 

Download proposal.mw

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