Carl Love

Carl Love

28055 Reputation

25 Badges

13 years, 3 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are replies submitted by Carl Love

@dbain1 The numerical integral is very difficult to evaluate for some values of p, such as 0.15. Your integrand looks a bit strange. Did you perhaps intend for it to be the absolute value of the expression that is to be raised to the 1/p power? If I include abs thus, there is no problem doing the numeric integration. Without abs, we end up raising negative values to fractional powers, obtaining imaginary results.

This is a topic that has been much argued on MaplePrimes. The appropriate question is How much effort (or time or expense) should be expended determing whether a given expression is identically zero when that problem is known to be in general undecidable? You may adjust that level of effort using TestZero. But no matter what you set it to, there will always be some singular matrix whose singularity can't be detected. Thus, if you call it a bug, there will necessarily always be a bug. I'd say that somthing that is mathematically impossible to fix shouldn't be called a bug. See the Wikipedia article "Richardson's Theorem" https://en.wikipedia.org/wiki/Richardson%27s_theorem

If you could provide a reference for an algorithm, I could probably implement it in Maple. I'm browsing the Wikipedia article "Forward Error Correction" https://en.wikipedia.org/wiki/Forward_error_correction, but the topic is too large.

Without going into much detail, I see that your final plot command fails because there is no numeric value for q. Is your q supposed to have a value?

@mapleatha Sorry, I misread your Question. I thought that you were getting unevaluated output and the blue was what you were expecting. My fault, totally. This Answer can just be ignored.

Is the domain of the function a finite set?

@kfli For some reason that I can't figure out, the freeze command is intentionally "crippled" to not work on names such as x[n]. This crippling is obvious in the very brief code that you can see via showstat(freeze). But it will work on x(n), which is considered a  function rather than a name. So try changing all x[n] to x(n). If you want these to display as subscripted, that can still be done. You can freeze a whole batch of expressions in a single command in a variety of ways. One example is

assign(seq(x(k)= freeze(y(k)), k= 1..9));

Let me know how that works for you.

@ecterrab I think that your impression that the OP is trying to differentiate wrt a function may be based on my example rather than on the OP's code. I just used that example because I didn't have the patience to decipher the OP's convoluted alias commands.

@macmp91 I think that you may only be looking at the eigenvalues when they are correct. In procedure Fn, immediately after the sort command, put the command

print(D[1]);

Then run the code until it errors out. Immediately before the error, I believe that you'll see the negative eigenvalue. When the sqrt is applied to this number, it becomes the imaginary value which causes the error message.

@macmp91 You should check for negative eigenvalues in procedure Fn. Do they make sense for this problem? If not, then you need to adjust something. Perhaps the input to Fn is not realistic. It is the negative eigenvalues that are causing the error message. Perhaps you intended to use singular values rather than eigenvalues. (The singular values of A are the square roots of the eigenvalues of A.A^+ and are guaranteed to be nonnegative.) See ?SingularValues. The main floating-point algorithm computes them without first computing the eigenvalues, which is more robust.

Actually sin(t)^p and (sin(t))^p  and even (sin^p)(t) mean the same thing in Maple. Personally, I prefer the form with the fewest parentheses. The operator binding between f and (x) in the Maple expression f(x) has the highest precedence and never requires extra parentheses. (But when using 2D Input, make sure not to put a space between f and (x).)

The famous mathematical constant Pi is capitalized in Maple; lowercase pi is just a variable.

Would you please attach a Maple worksheet showing the situation where you obtained a result containing Float(infinity)? You can use the green uparrow on the toolbar of the MaplePrimes editor to upload worksheets.

Often such a result comes because the answer to the problem is truly infinity. And often it comes due to round-off error causing division by a very small number. Usually we can determine which is the situation, and often the latter case can be corrected.

@acer I agree with Acer 100% on this. And I'll add that even if applyrule correctly did everything that it was supposed to do, it'd still be a fairly weak command compared to evalindets or subsindets.

@markweitzman 

Using the same technique that I used for the first problem, here are all the eigenvalues under 15 with plots of their solution curves. Note that I changed your piecewise to an equivalent Heaviside out of personal preference. I forgot to change "under 100" to "under 15" in the title (only) of the final plot
 

restart:

#Name the boundary points and boundary values:
(A, B, psiA, psiB, DpsiA):= (0, 2, 0, 0, 1):

ode:= diff(psi(u), u$2) = (2*Pi)^2*(Heaviside(u-1) - e)*psi(u):
ic:= psi(A) = psiA, D(psi)(A) = DpsiA:

Sol:= dsolve({ode, ic}, numeric, parameters= [e], abserr= 1e-13):

#The characteristic function:
F:= proc(ee)
   Sol(parameters= [e= ee]);
   eval(psi(u), Sol(B)) - psiB
end proc:

#This plot will guide us to the eigenvalues, which are its zeros.
plot(
   F, 0.1..15, numpoints= 999, view= [DEFAULT, -0.5..0.5],
   thickness= 2, size= [2000,500], color= orange, xtickmarks= 60,
   title= typeset("Zeros of characteristic function for ", ode),
   titlefont= [TIMES,24], gridlines= false,
   labels= [e, ``], labelfont= [TIMES,24]
);

 

#Candidate intervals for eigenvalues, based on the plot:
EI:= [0..0.25, 0.7..0.74, 1..1.25, 1.5..1.75, 2..2.5, 2.75..3, 3.5..3.75, 4.25..4.75, 5.5..5.75, 6.5..7, 8..8.25,
   9.25..9.75, 10.75..11.25, 12.5..13, 14.25..14.75]:

Digits:= 7:
E:= table():
for ei in EI do
   r:= fsolve(F, ei, fulldigits);
   if r::numeric then
      E[ei]:= r
   else
      printf("Failed to find the eigenvalue known to be in %a\n", ei)
   end if
end do:

#The found eigenvalues are
E:= sort(convert(E, list));

[.1843580, .7074703, 1.206517, 1.562500, 2.084468, 2.800468, 3.560780, 4.535382, 5.557413, 6.775081, 8.057593, 9.518408, 11.05834, 12.76398, 14.55906]

(1)

#Plot the corresponding solutions.
(m,M):= (min,max)(E): #Only needed to color the plots

PL:= table():
for e_k in E do
   Sol(parameters= [e= e_k]);
   PL[e_k]:= plots:-odeplot(
      Sol, [u, psi(u)], u= A..B, numpoints= 999,
      legend= [e= evalf[3](e_k)],
      color= COLOR(HUE, .8*(e_k - m)/(M - m))
   )
end do:

plots:-display(
   [entries(PL, 'indexorder', 'nolist')], gridlines= false, size= [1000,1000],
   legendstyle= [font= [TIMES,24]],
   thickness= 2,
   title= typeset("Solution curves for all eigenvalues under 15 for ", ode),
   caption= typeset("\n", e, " is the eigenvalue; hues are scaled linearly by ", e, "."),
   titlefont= [TIMES,16], axesfont= [TIMES,16],
   labelfont= [TIMES,32], captionfont= [TIMES,16]
);

 

 


 

Download ShootingForEigenvalues.mw

I think this may be handled by the Physics package, but I really don't know. So I updated your tags to include Physics so that someone who knows will notice this Question.

First 347 348 349 350 351 352 353 Last Page 349 of 709