@WD0HHU 


Change

dM1 := diff~(dM1, [a, c, d]);

into

dM1 := diff~(M1, [a, c, d]);

@WD0HHU 

Use the big green up arrow in the menubar to upload your worksheet (mw file).
Be careful, the name of mile must not contain special characters to be uploaded.

@WD0HHU 

Are the outputs of dM1 in your worksheet and mine (Maple 2015) the same?

If it is so try to replace this

Sol := solve(%, [a, c, d]):

by this

Sol := solve(dM1, [a, c, d]);

and tell me what you get.
Idealy upload your worksheet and send it to me.

@Alfred_F 

Observing that b is almost equal to 1, a lot of work can be done in a formal way.
It shows you that there is at least 2 solutions [a=+A, b=1, c=C, c=D] and [a=-A, b=1, c=C, c=D]. (note the + and - signs before A): 
(
I wasn't capable to find a closed form expression for 
)

NonlinearFit_details.mw






Note that if [a=A, b=B, c=C, c=D] is a soution, then  [a=A, b=B, c=C+2ℤ𝝅, c=D] is a solution too and 
[a=-A, b=B, c=C+2ℤ𝝅+𝝅, c=D] is another solution

 and writeto and save/read (see @acer 's answer) commands work exactly the same whatever your OS Unix/Linux,  Mac OSX, and even Windows.
More of this their related help pages are very clear and I can assure you I'm not an Unix guru.

About the with(plots) (p, not P) stuff: if you don't want to do some "special" plots invoking your polynomial function, or combine several of these plots, you don't need to load the plots package.
There is no in joke here.

This is the same thing for the Digits:=20 command (default is Digits := 10). Why do you use 20 digits? Is that because of some precision requirement, or did you started from some worksheet you found somewhere without thinking about the implications of this command (it enlarges the size of the machine representation of the polynomial and maay augment the computation time)?

Without this domain you can take xend=0 (which trivialy verifies constr) and find that wn(a, ..., xend) is identically null. So max(wn) = min(wn) = 0 whatever a, b, c,d, f, g.

@erik10 

Example of a torus

Feel free to uncomment the two first lines and comment the third one in

 #r := [seq(rand(0. .. 1.)(), i=1..3)]:
 #r := 2 *~ r /~ sqrt(add(r^~2)):

  r := [1.461755141, .1237707985, 1.359394240];

ProjectedBoundary.mw

@erik10 

Orthogonal projection of the trace ("Silhouette"?) of an ellipsoid cut by a plane (key function = plots:-intersectplot)
Silhouette.mw



You can dynamically change the point of view using this worksheet Silhouette_with_Explore.mw

 

Writing simultanously k*x and k(x-L) is likely a typo and that the latter should be written instead k*(x-L). This seems consistent with the exerpt of LaTex you provide (𝝀 replaced by k?).

Another possibility is you writting k*x instead of k(x)?
In which case integrand1 (and integrand3 too) are not correctly built, look to this simple example: 

my_k := r -> exp(r):
cos(k*(L - x)) = eval(cos(k*(L - x)), k = my_k(r));
lprint(%);
              cos(k (L - x)) = cos(exp(r) (L - x))
cos(k*(L-x)) = cos(exp(r)*(L-x))

#--------------------------------------------------------------

cos(k(L - x)) = eval(cos(k(L - x)), k = my_k(r));
lprint(%);
              cos(k(L - x)) = cos((exp(r))(L - x))
cos(k(L-x)) = cos((exp(r))(L-x))

#--------------------------------------------------------------

cos(k*(L - x)) = eval(cos(k(L - x)), k = (u -> my_k(u)));
lprint(%);
                cos(k (L - x)) = cos(exp(L - x))
cos(k*(L-x)) = cos(exp(L-x))

I believe you should clarify your notations.

Nevertheless, assuming my initial guess is correct, you can look to worksheet FTCT.mw and see some technical points (for instance the integral of integrand1 seems to diverge, and it would be useful thet you provid the ranges for x and t that you want to use).
One interesting point is that the integral are almost independent of t providing it remains "small", let's say t < 10⁷.

@JoyDivisionMan 

I answered your question arround midnight yesterday and made a mistake in what comes after the attached file.
_________________________________________________________________________________________

Once said, its easier (IMO)  to think to the problem in terms of probability: a sum of iid random variables instead of a auto-convolution product.

The result I use in the attached file is the following:

• Let X a random variable with PDF f_X(t) and let Y = c*X (assuming c > 0)
  Then the PDF f_Y(t) of Y writes f_X(t/c) / c.
  Is the support of X is (0, a), then the support of Y is (0, a / c)

Application:
The function f(t) you define can be viewed as the pdf of a continuous random variable with support (0, 2).
If U and V  are iid random variables with pdf f(t), then U*c and V*c (still c > 0) are iid random variables with support (0, 2/c) and U and V pdfs are equal to f(t/c)/c.
Then U+V is a random variable whose support is (0, 4/c).

So, without changing the expression of f(t):
You can directly construct U and V while saying their distribution has PDF f(t/c) / c:

U := RandomVariable(Distribution(PDF = (t -> f(t/c) / c))):

# idem V

In your case c = "desired upper bound of the support" / "present upper bound of the support of f(t)" =  (2/sqrt(3)) / 2 = 1/sqrt(3)



Here is the corrected file 

Download SphereConvolution_sand15_corrected.mw

@Alfred_F 

Raw version:

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

ftest := (1/2)*(diff(xtest(t), t))^2+u^2*cos(xtest(t))+xtest(t)*v*sin(t):

C := collect(ftest, t, factor):

J := proc(A0, A1, B1, A2, B2, A3, B3)
  evalf(Int(eval(C, [a0, a1, b1, a2, b2, a3, b3] =~ [A0, A1, B1, A2, B2, A3, B3]), t=0..2*Pi))
end proc:


Optimization:-NLPSolve(J)

0.5568134237e-6})]

Add variation domains for all the parameters to get a more reliabable result.
Change the starting point ('initialpoint') to check if the location of the minimizer changes or not... or use a global optimizer.

But before thinking to using a global method if you don't have any at your disposal, avery simple idea you could use to explore how J behaves is to sample (Statistics:-Sample) the 7-hypercube a0, ...b3 and observe how J behaves in this high dimensional space (there exist statistical methods to get an idea of that).

After running a quite naive algorithm J(Pi/2, 0$6) = 0.
Do this

plot(J(z, 0$6), z=-4*Pi..4*Pi)

A few important results from statistical exploration:

1. The minimum of J(a[0], A) conditionally to  A=(a[1], ..., b[3])  is obtained for a[0] = 𝜋 + 2𝜋ℤ.
    
2. The maximum of J(a[0], A) conditionally to  A=(a[1], ..., b[3])  is obtained for a[0] = 2𝜋ℤ.
    
3. J( 𝜋/2 + 2𝜋ℤ, A) = 0 whatever A.
    
4. The global minima are located at 
   a[0] = 𝜋 + 2𝜋ℤ
   b[1] = 0.09175109
   a[1] = a[2] = a[3] = b[2] = b[3] = 0
   where  J = -0.5798982790

@janhardo @jalal 

More general rotations can be easily produced using the more general Euler's angle formalism.

Example EulerAngleFormalism.mw

@callumneily 

See Wiki and reference herein.

 

Download Pi_and_random_process.mw

@callumneily

 

Download BuffonNeedle.mw

(Circa 1910)

The number N of needles which intersect the edge of the slats is proportional to the needle length a and inversely proportional to the width b of the slats. So it can be written as N = P×a/b where P is some proportionality constant of which the value is to be determined.

Instead of considering straight neeedles, Borel considers circular needles of diameter d. Their perimeters, or lengths, are equal to 𝜋×d. So, with the notation above a = 𝜋×d.
Then Borel choses the particular value d = b which leads to N = P×𝜋×b/b ⟹ N = P×𝜋.
But any circle of diameter d either cuts a given slat edge twice or is tangent to two consecutive edges... which is an event of null probability.

Thus 2 = P×𝜋 ⟹ P = 2/𝜋.