Carl Love

Carl Love

28035 Reputation

25 Badges

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

MaplePrimes Activity


These are replies submitted by Carl Love

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

This is an updated version of the worksheet above. I show how to obtain (by the same shooting method) all the eigenvalues in an interval, and I do it and plot it for all under 100.

Note that these plots are much clearer and have more annotation when viewed in the worksheet than when they are viewed here on MaplePrimes.
 

restart:

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

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

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

#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..100, view= [DEFAULT, -100..100],
   thickness= 2, size= [1000,100], color= orange,
   title= typeset("Zeros of characteristic function for ", ode),
   titlefont= [TIMES,24],
   labels= [e, ``], labelfont= [TIMES,24]
);

 

#Candidate intervals for eigenvalues, based on the plot:
EI:= [0..2, 4..6, 10..12, 16..18, 22..25, 30..33, 39..41, 47..50, 56..59, 66..68, 76..78, 86..88, 96..98]:

E:= table():
for ei in EI do
   r:= fsolve(F, ei);
   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));

[.6679862422, 4.696795558, 10.24430877, 16.71189010, 23.88999426, 31.65945727, 39.94141775, 48.67906755, 57.82916534, 67.35751908, 77.23651840, 87.44459180, 97.96981713]

(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 100 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

@markweitzman

The name e is used for two variables in the scope of procedure F: It is a global used in the ODE, and it is the procedure parameter of F. These are two different variables that just happen to have the same name. In such situations, the global variable can be distinguished by using the prefix :-. I could just as well have used a different name for the procedure parameter, like ee. Then the statement could be Sol(parameters= [e= ee]), or Sol(parameters= [:-e= ee], or as you pointed out simply Sol(parameters= [ee]), (but not then Sol(parameters= [e])).

The algorithms for solving IVPs and BVPs are very different. As Mariusz's Answer shows, the BVP needs to be solved separately with a different approxsoln for each eigenvalue. On the other hand, the IVP only needs a single call to dsolve from which all the eigenvalues in any interval can be obtained. In a Reply below, I've updated the worksheet above, and I show how to obtain all the eigenvalues in an interval, and I actually do it and plot it for all under 100.

 

@Mariusz Iwaniuk For numeric BVPs, dsolve allows parameters in exchange for extra BCs. The BCs are counted at the syntax level, and it won't let you get away with giving the wrong number of them (even if perchance the solution satisfies an extra condition).

@_Maxim_ The nested context is what I would hope for. I can't now produce an example where the nested evalf doesn't work. I wonder if it's something that's been fixed in a recent version? I'm using Maple 2016 at the moment. Anyway, I withdrawal my first comment above: "evalf's nonstandard evaluation."

@Joe Riel wrote:

  • It would be interesting to see the score of some of the spam that used to plague this site but, thankfully, no longer appears.

Okay, here's a great candidate that I deleted from this site a few minutes ago:

spam:= "The orthopedic treatment is identified with the musculoskeletal framework for illustrations bones .the orthopedic specialists utilize surgical and non-surgical medications. The non-surgical means the turn expire innate scatters and so forth. The ortho implies straight and paidion implies tyke. In days of yore the orthopedic treatment is basically utilized as a part of military reason little damage, and conceived Brocken and so forth. Arthroplasty is the principle part of orthopedic treatment and it is supplanted the musculoskeletal is supplanted. We can supplant the joints are knee, hip, bear, elbow, wrist, lower leg, spine, and finger joints .Epidemiology is another kind of supplanting of joint. be that as it may, it is confounded for underneath 18 years or more 80.After surgery take 1 month bed rest for in clinic. The orthopedic games solution is accessible on wherever on the planet. Furthermore, every games in your city. Orthopedic games drug is one of the particular callings. It is worked in gathering and they are likewise centered just the wellbeing of a competitor.":

seq(
   X = StringTools:-Readability(spam, method= X), 
   X = op(op([1,2,1,2,-1], eval(StringTools:-Readability)))
):
evalf[3]([%]);

[GunningFog = 15.5, FRES = 43.2, FKGLF = 11.1, ARI = 16.0, ColemanLauIndex = 21.2, SMOG = 13.4]

What kind of person would let someone who writes like that perform surgery on them?

@vv That's true for functions, but it's not true for basic arithmetic. For an example, try computing the residuals from the floating-point roots of a high-degree (50 or greater) polynomial.

@John Fredsted 

typeset can be used with caption just as it can be used with annotation.

@_Maxim_ Because of evalf's nonstandard evaluation rules, it doesn't work to do

evalf[d1](evalf[d2](e))

for d1 < d2, even though it's frequently what I want to do (to do a computation at high precision and only see a few digits of the result). You need to do it as

x:= evalf[d2](e):  evalf[d1](x);

which grates on my nerves.

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