Carl Love

27187 Reputation

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

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

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[8]):
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[2], 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));

 >

The horizontal scroll bar is only going to be activated when you have an expression that the GUI cannot figure out how to line break, such as a fraction with a long bar.

You can put the Matrix into a Spreadsheet and then scroll in both dimensions, resize the cells, and scroll within the cells to your heart's content. Let's say that your Matrix is M. Then

SetMatrix(ssid, M);

You resize the spreadsheet window by dragging on the lower right corner. (It is not clearly marked!) You resize cells by dragging on their lower right corners (which are clearly marked). Resizing a cell changes the width of every cell in its column and changes the height of every cell in its row.

By the way, I don't understand the title of this Question, "Vertical display arrow".

Changing from Matrix to listlist...

For what it's worth, it's trivial to take an existing plot and convert the Matrix to a listlist (which is the Mapleish name for the structure is your older plot). Let's say that plot is stored at P. Then

evalindets(P, Matrix, convert, listlist);

I am curious if others think that this changes the appearence of the plot. I'm not sure. I am sure that your M12 plot has better anti-aliasing. (Look very closely at the section -1..1 on the M12 plot. Depending on the quality of your eyesight, you may need to export the plot to a JPEG viewer and zoom in.). The M17 plot may not have any anti-aliasing. Can you check the anti-aliasing setting on both versions of Maple? But, alas, I don't think that anti-aliasing is the complete answer to the difference. The M17 plot is too nasty for that to be the only problem. (And we discussed this at length with respect to plot(sin(x), x= 0..2*Pi) a month or two ago.)

RootOf(..., index= 1)...

Take your expression newpar (the first, unerroneous one) and apply the following transformation:

evalindets(newpar, 'RootOf'(algebraic), allvalues@(e-> RootOf(op(e), index= 1)));

Let me know how that goes. I couldn't test it on your actual code because it wasn't attached to your post.

Histogram, Sample, binbounds...

To make an actual histogram, we take a random sample of the Poisson distribution and use Statistics:-Histogram. I chose a sample size of 2^13. The default is that the bars (bins) are bounded by the integer values. That doesn't make sense to me. I want the bars centered on the integer values. So I set the binbounds as n-1/2..n+1/2.

Take the previous code, and simply change the line defining H[n] to

`H[n]:= Histogram(Sample(Poisson(n), 2^13), binbounds= [seq](-1/2+k, k= 0..25));`

DensityPlot and typeset...

Statistics:-DensityPlot will plot an ordinary plot of the PDF for a continuous distribution, and it will plot something like a histogram for a discrete distribution. So, it can be used for both plots. If you want an actual histogram, we can draw a random sample from the Poisson and use Statistics:-Histogram. (Let me know.) I also updated your string plotting to modernity with typeset.

Note that the second parameter to Normal is not the variance; rather, it is the square root of the variance. Also, the expressions "NormalPDF", "PoissonPDF", and "ProbHist" are meaningless in Maple.

`with(plots): with(Statistics):for n from 1 to 15 do      tracker[n]:= textplot([18,0.3,typeset(lambda = n)], color= blue);      H[n]:= DensityPlot(Poisson(n));      N[n]:= DensityPlot(          Normal(n, sqrt(n)), range = 0..25, color= "Niagara Red"     );     P[n]:= display([H[n], N[n], tracker[n]]) od: display(     [seq](P[n], n= 1..15), insequence=true,     title= typeset(          "Normal approx. to the Poisson. ",           lambda,          " is increasing from 1 to 15."     ));`

Make assumptions...

Try solve(q,z) assuming z>0 and solve(q,z) assuming z<0. I get two solutions for each.

 First 362 363 364 365 366 367 368 Last Page 364 of 389
﻿