Carl Love

## 26084 Reputation

11 years, 51 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

## Spherical coordinates...

The integral is too difficult for Maple in its current form. Convert it to spherical coordinates.

I say this with trepidation, considering how badly I erred with the limits on your previous question: Are you sure that that set R is correctly specifed? The first condition, z <= sqrt(x^2+y^2+z^2), is satisfied by any three real numbers. The second condition is strange because there is no lower limit to z. I wonder if it was supposed to be abs(z).

## 0...

The integral is easily seen to be 0 without the help of Maple.

2*y >= x^2 + y^2  -->  0 >= x^2 - 2*y + y^2 = (x-y)^2  -->  (x-y)^2 = 0  -->  x=y.

Thus D is meager, and the integral is 0.

## Variance = (Standard deviation)^2...

Variance(X)=sqrt(StandardDeviation(A)).

That should be

Variance(X) = StandardDeviation(A)^2.

Also, consider using ?KernelDensity rather than EmpiricalDistribution, which gives a discrete distribution. Then you can make a custom ?Distribution by uing the function returned by KernelDensity as the ?PDF.

## Differ by a constant...

The two answers differ by a constant, so they are both valid forms of the integral. The constant is ln(-1) = Pi*I. You can take the derivative of each and see that in both cases you get the integrand.

## parameterized colours with colour= HUE(....

It is not a silly question. The easiest way to get parameterized colours is to use colour= HUE(x) where 1 >= x >= 0. The HUE scale is the standard colour spectrum, from red=0 to violet=1.

Here's an example. I've generated example data such as you describe, but using random data since I don't have your specific vectors.

`V:= unapply(Vector(16, ()-> randpoly(x, degree= 1)), x):V||(1..4):= 'V(k)' \$ 'k'= 1..4:So now I have 4 vectors, V1, ..., V4 of 16 numeric elements each.plot([seq]([seq]([k, (V||k)[j]], k= 1..4), j= 1..16), colour= [seq](HUE(j/16), j= 1..16));Note that the argument to HUE is matched to the index into the Vectors. `

## both sides?...

I don't know what you mean by "both sides". Do you mean to solve it for each variable, T0 and t? That can be done like this (I've changed the decimals to exact fractions, but it is not necessary to do that):

```

>

eq:= (T+exp(-sqrt(t/1000)))/T = 95/100; >

solve(eq, {T}); >

solve(eq, {t}); >

## Detailed solving example...

Solving a multi-branch equation with solve.

It's a good idea to start most worksheets with a restart.

 > restart;

Balance parentheses and remove curly braces {}.  This is mostly done by eye and hand. Maple offers minimal help, like showing you where the match for a parenthesis is. This helping is done while you are typing, so I can't show it here.

Give the equation (actually it's an expression as it has no equals sign) a name.

 > eq:= 1+(25.*phi^3-297.*phi^2+1162.*phi-1500.)*(840.+332.*theta+29.*theta^2)/((phi^4-16.*phi^3+89.*phi^2-201.*phi+150.)*(3.+theta)*(28.+11.*theta+theta^2)); Convert decimals to exact numbers.

 > eq_exact:= convert(eq, rational); Solve the equation (assumed to be set to 0) for theta. (Lengthy output suppressed.)

 > Sol:= [solve](eq_exact = 0, theta):

Count the solutions. I expect three because the original would be degree three in theta if you multiplied it out.

 > n:= nops(Sol); Make the solutions a function of phi.

 > f:= unapply(Sol, phi):

Example of  the solutions evaluated at an exact numeric value.

 > f(0); Convert those expressions to decimal form.

 > evalf(%); Convert spurious imaginary parts to 0.

 > fnormal(%, Digits-2); Clean out the zeroes.

 > simplify(%, zero); Generalize the above process: Make three functions for the numeric evaluation of the three branches of the solution.

 > F||(1..n):= 'phi-> simplify(fnormal(evalf[Digits+3](f(phi)[k])), zero)' \$       'k'= 1..n:
 > plot(F1, -2..2, discont); > plot(F2, -2..2, discont); > plot(F3, -2..2, discont); A different approach---more numerically oriented with fsolve.

Convert the expression to a polynomial by taking its numerator. This technique is only valid because we assume that the expression is set equal to 0. This step is not essential, but fsolve works better with polynomials than with general equations.

 > eq1:= numer(eq_exact); Collect the coefficients with respect to the main variable. This step is also not essential.

 > eq2:= collect(eq1, theta); Make it a function of phi.

 > Eq2:= unapply(eq2, phi):
 > Sol:= phi-> fsolve(Eq2(phi), theta); > Sol(0); The fsolve technique avoids the spurious imaginary parts.

Compare with the three solutions from the solve method.

 > F||(1..3)(0); They're identical except for the order.

It is more difficult to convert the fsolve solution into a meaningful plot.

## Here's a procedure for them...

Working from the Wikipedia pages "Gauss-Hermite quadrature" and "Hermite polynomials", I whipped up this procedure for you. It takes a positive integer argument n and returns two lists of elements, the first being the nodes and the second being the weights.

`GaussHermite:= proc(n::posint)local      x,     Hm:= unapply(numapprox:-hornerform(orthopoly[H](n-1, x), x), x),      R:= [fsolve](orthopoly[H](n,x), x),     C:= evalf(2^(n-1)*n!*sqrt(Pi)/n^2);     R, map(x-> C/Hm(x)^2, R)end proc:     GaussHermite(6); [-2.35060497367449, -1.33584907401370, -0.436077411927616,    0.436077411927616, 1.33584907401370, 2.35060497367449],    [0.00453000990550894, 0.157067320322850, 0.724629595224395,    0.724629595224395, 0.157067320322850, 0.00453000990550894]`

## Bug in arbitrary-order diff---use add...

Your question is not boring, and you've done nothing wrong. You've stumbled upon an obscure bug in Maple's arbitrary-order derivatives.

Are you really trying to "solve" the formula in some way? Or do you just want to evaluate it for specific integer values of r? If the latter:

Maple tries to take the derivative of arbitrary order that occurs in your sum. The formula that it returns is apparently wrong (very wrong, in fact). We can work around that by using add instead of sum.

`Mrs:= proc(s, r::posint, x)local      i, X,     S:= s+1,     P:= i-> i+s,     T:= i-> 2*i+s,     U:= i-> 2*i+s+1,     R:= t-> (2*t)!/t!/4^t,     A:= R@P, B:= R@T, C:= R@U;     eval(          add(               (-1)^i*(A/B*C)(i)/T(i)!/i!/4^i*               diff((1-X^2)^T(i), [X \$ 2*i]),               i= 0..r-1           ),          X= x     )end proc:`

If you are trying to "solve" the formula in some way other than simple evaluation, let me know. My procedure does allow you to use symbolic or non-integer s, and some significant and interesting simplifications are possible (try s = 1/2).

## Unbalanced quotation marks...

You wrote:

While trying your solution, I noticed that my maple gives me a higher order of the complex part of the solution (E-6).

That's because I had Digits set to 15 and you had Digits set to 10. You can include the command Digits:= 15 right after your restart.

You wrote:

Furthermore I cannot plot my solution due to the following error....

Your plot command has unbalanced quotation marks, and possibly some other bad character(s). I'm surprised that it got past the parser. I recommend that you retype (not cut and paste) the whole command in a new execution group. All of the quotes need to be double quotes ("), not single quotes ('). I also recommend that you use input mode Maple Input (like the red characters that I used) rather than 2D Input, because with Maple Input both you and I can see exactly what you typed. With 2D Input it is impossible to visually distinguish a double quote from two single quotes next to each other.

You wrote:

Furthermore I get an error trying to find the slope of the catenary at x=0.

Use eval(Helling, x= 0).

You wrote:

I was wondering why maple supplies us with two solutions for the differential equation.

There is only one solution to the differential equation proper; it is only after the boundary conditions are put in that the branching occurs. You can see that the differential equation has one rather simple solution by calling dsolve without the boundary conditions:

`dsolve(DIV);                       cosh(_C1 qT + qT x)                      y(x) = ------------------- + _C2                               qT               `

You wrote:

Is it because of the square root in the system?

No. When you apply the boundary conditons to the above solution to solve for the constants of integration _C1 and _C2, you get an equation of the form cosh(x) = B. Such equations generally have two solutions (for B > 0) because horizontal lines cross cosh curves twice.

You wrote:

Furthermore, what is the origin of the complex part of the solution?

After the boundary conditions are applied, the solution contains an expression of the form cosh(A+ln(RootOf(Q))) where Q is a quadratic with two real roots, one positive and one negative. (The nature of the roots depends on your parameter values and boundary conditions. For the case we are discussing, one is positive and one is negative.) The ln of a negative number is complex; specifically, for x < 0, ln(x) = ln(|x|) + Pi*I.

It is important to note that this is not a complex solution; we are not ignoring any part of the true solution. The Pi*I turns real after the cosh is applied to it. It is only because of the inexact floating-pointing arithmetic that we see it at all. That's why I call them spurious imaginary parts. Below, I redid the worksheet in such a way that exact arithmetic is used up to the point that the complex parts disappear.

 > > Digits:= 15:
 > > >  > >  > > Opl3:= allvalues(Opl2):
 >  >  >  >  > > eval(Helling, x= 0); >

## Vector-vector product?...

Your last command has a product of two vectors, each with five elements (assuming m = 4), which you are trying to assign to a vector. But the product of two vectors ordinarily means the dot product, which yields a scalar (a single number, not a vector). Did you mean for that to be the elementwise product of the two vectors? The command for that would be

U[j+1, ..]:= U[j, ..] *~ A;

(assuming that you're using Maple 13 or later. If you have an earlier Maple, let me know.)

Also, the loop should go for j from 0 to n-1, not for j from 0 to n.

## series at infinity...

You asked that, and it was answered, within the past week. If ex is your expression, then

series(ex, c= infinity, 3);

Did that answer not work for you?

Thanks for using a shorter title.

## Array...

You cannot use an index of 0 in a Maple Matrix. Matrix indices always start at 1. But you can use 0 indices in the closely related structure called Array.

h:= 1/4:
m,n:= 4,4:
U:= Array(0..m, 0..n, datatype= float):
U[0, ..]:= < seq(evalhf(cos(Pi*i*h)), i= 0..m) >;

Note that I used evalhf instead of evalf. The usual evalf would work also, but I think that evalhf is a little better in this situation.

## Eliminating spurious imaginary parts...

The imaginary parts in your solution are spurious. They arise because the solution contains expressions in the form cosh(a+Pi*I), except that the Pi is a decimal number. In exact arithmetic, that expression simplifies to -cosh(a) (via command expand), but in floating-point arithmetic there is a residual imaginary part that is nearly zero. These residual near-zero parts can be converted to zero with fnormal, and then the zeros can be removed entirely with simplify(..., zero).

Here is your worksheet with my modifications.

 > > > > qT0:= q0/T0:

Moved T to other side of equation. Solve equation with symbolic q /T.

 > > >  These RootOfs are quadratic (see the _Z^2). Apply allvalues to get two solutions.

 > Apply the numeric value of the parameter to the two branches of the solution.

 >  > Sol2:= eval(Opl2, qT= qT0); Sol2 has a spurious imaginary part which would disappear upon expansion in exact arithmetic. Here, we compensate for the inexact decimal arithmetic.

 > expand(Sol2); > fnormal(%, 11); > Sol2:= simplify(%, zero); >  Now you have to decide which solution you want.

 >

## Integral on z=0..1. No solutions for z<0...

Here is a solution for the integral on 0..1, and some plots to verify that there are no solutions for z < 0. I now realize that this is essentially the same answer as Preben's, but I arrived at it independently before reading his. My plots are significantly different, showing the same function behaviour from a different perspective.

Important: The 3D plots are much more convincing when viewed in the worksheet than when viewed on MaplePrimes!!

 > restart:
 > with(Statistics):
 > dummy:= sqrt((1+y*e)/(y*e)):
 > dummy1:= sqrt(1/(1+y*e)):
 > Q:= RandomVariable(Normal(0, dummy)):
 > R:= RandomVariable(Normal(y*e*z/(1+y*e), dummy1)):

I changed the constants to exact rationals just to keep things neat.

 > y,e,k,l,a:= 1, 1/5, 1, 1/2, 5/4:

The inner integral:

 > II:= Int(k*r*PDF(R,r), r= -infinity..x); > IIV:= value(II); It probably makes no difference to the final answer, but I'd like to do a few steps of simplifying to that.

 > expand(%); > combine(%); > simplify(%); > ii:= %:

We intend to perform an fsolve on the above expression, specifically, to set  it to 0 for specific values of z and to solve for x. A plot3d will help us to understand where the expression is equal to 0.

 > plots:-display([      plot3d(ii, z= -1..1, x= -4..4),      plot3d(0, z= -1..1, x= -4..4, style= wireframe, color= black) ]); It looks like there are no real solutions for z < 0, and a unique real solution for z > 0.  But let's look more closely at the possibility of a solution with x < 0 for z < 0.

 > plots:-display([    plot3d(ii, z= -1..0, x= -5..-2.5),    plot3d(0, z= -1..0, x= -5..-2.5, style= wireframe, color= black) ]); It looks like the function approaches the plane very rapidly, but never crosses it. It's looking like there are no real solutions for z < 0. Let's look at some slices.

 > plot([seq](eval(ii, z= k), k= -1..-0.2, 0.2), x= -5..-2.5); We conclude that there are no real solutions for z < 0. We will still try the integral on 0..1.

The following line is the key step to removing your "z is in the equation and is not solved for" error.

 > iiu:= unapply(ii, z):

Construct the function inversion that occurs in the outer integrand. RootOf allows some analysis of the inverse; whereas fsolve has more flexibility numerically.

 > OIR:= z-> RootOf(iiu(z), x):
 > OIF:= z-> fsolve(iiu(z), x):
 > plot(OIF, 0..1); Note the possibility of a singularity at z = 0.

 > OIR(0); Clearly +infinity and -infinity are the roots of that. So there is a singularity. But since the tails of exp(-x^2) are integrable, so is the above singularity.

Let's look at the other factor in the outer integrand.

 > PDF(Q,z); > plot(PDF(Q,z), z= 0..1); It is smooth, and it's not going to change the nature of the singularity (since it is real and non-zero at z=0).

 > plot(z-> OIF(z)*PDF(Q,z), 0..1); So, a rough eyeball guess of the integral is rectangle + triangle = 1*0.2+1*(0.4-0.2)/2 = 0.3

The final integral (this takes about 30 seconds):

 > evalf(Int(z-> OIF(z)*PDF(Q,z), 0..1)); >