Carl Love

Carl Love

27656 Reputation

25 Badges

12 years, 125 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are answers submitted by Carl Love

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.

restart

q0 := 2.77:

T0 := 1850.4:

qT0:= q0/T0:

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

DIV := diff(y(x), x, x) = qT*sqrt(1+(diff(y(x), x))^2):

RV := y(0) = 0, y(100) = -1008:

Opl := dsolve({DIV, RV}, y(x));

y(x) = cosh(qT*x+ln(RootOf(((exp(100*qT))^2-exp(100*qT))*_Z^2+2016*qT*exp(100*qT)*_Z-exp(100*qT)+1)))/qT+(1/2)*(-RootOf(((exp(100*qT))^2-exp(100*qT))*_Z^2+2016*qT*exp(100*qT)*_Z-exp(100*qT)+1)*exp(500*qT)+2*RootOf(((exp(100*qT))^2-exp(100*qT))*_Z^2+2016*qT*exp(100*qT)*_Z-exp(100*qT)+1)*exp(400*qT)-2016*qT*exp(400*qT)+4032*qT*exp(300*qT)-2016*(exp(qT))^200*qT-RootOf(((exp(100*qT))^2-exp(100*qT))*_Z^2+2016*qT*exp(100*qT)*_Z-exp(100*qT)+1)*exp(200*qT)-2016*qT*exp(200*qT)+2016*(exp(qT))^100*qT-(exp(100*qT))^2*RootOf(((exp(100*qT))^2-exp(100*qT))*_Z^2+2016*qT*exp(100*qT)*_Z-exp(100*qT)+1)+2016*(exp(100*qT))^2*qT+exp(100*qT)*RootOf(((exp(100*qT))^2-exp(100*qT))*_Z^2+2016*qT*exp(100*qT)*_Z-exp(100*qT)+1)-2016*exp(100*qT)*qT)/(((exp(qT))^100-1)^2*(exp(100*qT)-1)*qT*exp(100*qT))

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

Opl2 := allvalues(Opl):

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

Sol1 := eval(Opl2[1], qT = qT0);

y(x) = 668.014440433214*cosh(0.149697362732382e-2*x-3.08007158591335)-7283.33339213385

Sol2:= eval(Opl2[2], qT= qT0);

y(x) = 668.014440433214*cosh(0.149697362732382e-2*x+2.93037422318077+3.14159265358979*I)+6275.33339186040

Sol2 has a spurious imaginary part which would disappear upon expansion in exact arithmetic. Here, we compensate for the inexact decimal arithmetic.

expand(Sol2);

y(x) = -6275.33339204111*cosh(0.149697362732382e-2*x)+(0.202069600471683e-10*I)*cosh(0.149697362732382e-2*x)-6239.67674552448*sinh(0.149697362732382e-2*x)+(0.203224327649008e-10*I)*sinh(0.149697362732382e-2*x)+6275.33339186040

fnormal(%, 11);

y(x) = -6275.3333920*cosh(0.14969736273e-2*x)+6275.3333919+0.*I-6239.6767455*sinh(0.14969736273e-2*x)

Sol2:= simplify(%, zero);

y(x) = -6275.3333920*cosh(0.14969736273e-2*x)+6275.3333919-6239.6767455*sinh(0.14969736273e-2*x)

plot(`~`[rhs]([Sol1, Sol2]), x = 0 .. 100, title =

Now you have to decide which solution you want.

 

 

Download SpuriousImag.mw

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

Int((1/10)*r*2^(1/2)*5^(1/2)*6^(1/2)*exp(-(3/5)*(r-(1/6)*z)^2)/Pi^(1/2), r = -infinity .. x)

IIV:= value(II);

(1/180)*5^(1/2)*3^(1/2)*(z*Pi^(1/2)*5^(1/2)*3^(1/2)*erf((1/5)*5^(1/2)*3^(1/2)*x-(1/30)*z*5^(1/2)*3^(1/2))+z*Pi^(1/2)*5^(1/2)*3^(1/2)-30*exp(-(3/5)*x^2+(1/5)*x*z-(1/60)*z^2))/Pi^(1/2)

It probably makes no difference to the final answer, but I'd like to do a few steps of simplifying to that.

expand(%);

(1/12)*z*erf((1/5)*5^(1/2)*3^(1/2)*x-(1/30)*z*5^(1/2)*3^(1/2))+(1/12)*z-(1/6)*5^(1/2)*3^(1/2)*exp(-(3/5)*x^2)*exp((1/5)*x*z)*exp(-(1/60)*z^2)/Pi^(1/2)

combine(%);

-(1/6)*15^(1/2)*exp(-(3/5)*x^2+(1/5)*x*z-(1/60)*z^2)/Pi^(1/2)+(1/12)*z*erf((1/5)*15^(1/2)*x-(1/30)*z*15^(1/2))+(1/12)*z

simplify(%);

-(1/12)*(-erf((1/5)*15^(1/2)*x-(1/30)*z*15^(1/2))*Pi^(1/2)*z+2*15^(1/2)*exp(-(1/60)*(6*x-z)^2)-z*Pi^(1/2))/Pi^(1/2)

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

RootOf(15^(1/2)*exp(-(3/5)*_Z^2))

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

(1/12)*2^(1/2)*6^(1/2)*exp(-(1/12)*z^2)/Pi^(1/2)

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

.282975811112918

 

 

Download FsolveInt.mw

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

with(Spread);
ssid:= CreateSpreadsheet();
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".

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

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.

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

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."
     )
);

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

f:= (x,y)-> 5/(1+x^2+y^2):
X:= t-> 2+3*cos(t):  Y:= t-> -1+2*sin(t):
plots:-display([
     plot3d(f(x,y), x= -4..7, y= -4..4),
     plots:-spacecurve(
          [X(t), Y(t), f(X(t),Y(t))], t= 0..2*Pi,
          color= black, thickness= 3
     )
]);

The coefficient is not a matrix; it is the determinant of a matrix, which is a scalar. I think that finding a closed form for this sum is well beyond the capability of any current CAS.

Here's a procedure for generating the determinant, which might help you in finding a closed form for the determinant at least. The determinants are not extremely complicated.

DetAn:= (n::nonnegint)-> LinearAlgebra:-Determinant(
     Matrix(
          n, n,
          (i,j)->
               if j >= i and (j-i)::even then
                    (j-i+1)*(j-1)!/(i-1)!*a(j-i+1)*x
               elif i-j = 1 then  -1
               else  0
               end if
     )
):

Part of the problem is that your Matrix M is created inside the i and j loops. I say "part" because when I take the creation outside the loop, there are new problems that pop up that I don't understand yet. But certainly it cannot possibly work with the Matrix being created inside the loop.

Could you isolate the relevant section of code for me and place it in a plain text file? From the final restart to the end of the file. I have a lot of trouble working with 2d input, especially with a section of code of this size.

You need some {} in your last command. Change dsolve(sys2, [x1(t), x2(t), x3(t)]) to dsolve({sys2}, [x1(t), x2(t), x3(t)]).

If ex is your expression, then simply

series(ex, c= infinity, 3);

I guess this is a bug in solve. I don't know about avoiding it, but you can easily correct if after the fact like this:

solve(f(x)=y, x);

subsindets(%, list, op);

I have worked a little with your code. I cannot back-up your claim that the numeric rank is taking any significant amount of time. I haven't found what's taking the majority of the time, but it's not the rank. The rank for a 1000 x 50 or so Matrix from your program takes about 5 millisecs.

To find exactly how the time is divided among the various parts of your program, use ?exprofile. And to count how many times each part of your program is executed, use ?excallgraph. I haven't had a chance yet to use these on your program.

Symbolic rank is a another matter. Computing the symbolic rank of a dense 6x6 matrix of rational functions is a formidable task. For matrices of your size, it is out of the question. But, do you really need an ironclad mathematical proof that the rank is a certain value? When you evaluate a symbolic matrix at random values for the symbols, the probability that the rank of the numeric matrix differs from the rank of the symbolic matrix is infintesimal---it's like the probability of picking 0 when picking a random real from the interval (-1,1). But you should randomly generate real reals, or real floats, so to speak, rather than the 2-decimal-place ones that you are using. And there's no need to check the numeric rank multiple times. Two is enough. If the numeric rank differs (which might be due to numeric instability) between the two, then try a third.

First 366 367 368 369 370 371 372 Last Page 368 of 392