Re: Piecewise functions definition and usage

*To*: mathgroup at smc.vnet.net*Subject*: [mg24153] Re: [mg24098] Piecewise functions definition and usage*From*: Andrzej Kozlowski <andrzej at tuins.ac.jp>*Date*: Wed, 28 Jun 2000 02:12:05 -0400 (EDT)*Sender*: owner-wri-mathgroup at wolfram.com

on 6/27/00 10:54 PM, Viorel Ontalus at vio2 at Lehigh.EDU wrote: > Thanks for help, > Indeed in the first example I have a typo, it should be: > > Clear[f] > f[x_] := Sqrt[ 1 - ( (x - 3)/10)^2 ] /; > Abs[x - 3] < 10 > f[x_] := 0 /; > Abs[x - 3] > 10 ; > NIntegrate[f[x], {x, -1, 17}] > > As you noticed Mathematica gives a result, but also an error message. And it > is > no singularity with this definition ! Here I show a simple example but I used > this in more complex calculations, and there the program does not compute an > answer. Imagine that you have few integrals like this and the program will not > work It should be a way around this otherwise Mathematica has a serious flaw. > It means that it can not deal with piecewise functions properly. > My big problem is that I have to find a solution to this. The way I see it you have three problems. First you lack enough knowledge of mathematics to deal properly with the problem you are trying to solve. Secondly you do not have enough Mathematica programming skills to write a program that avoids mathematical pitfalls inherent in this type of situation. And thirdly you blame others (Mathematica) for these shortcomings! Let's start with mathematics. As I wrote in my first reply to your posting your corrected code is, in effect, no different from the original one. And the problem is still the singularity at 13. Yes, singularity. You do not seem to know what this means. Admittedly, this word is used in mathematics in several somewhat different though related meanings (sometimes dis-continuity, sometimes non-differentiability of some order, sometimes derivative = 0) but in this case it means a point where the function is not differentiable. Although it is not required that a function be differentiable for the integral to exist, the presence of such points can cause trouble in numerical integration. To see that you should be familiar with at least the simplest methods of computing integrals numerically, e.g. the trapezoid rule. To compute such an integral numerically you have to subdivide the range over which you are computing the integral into intervals, over which you construct your trapezoids. If you think about this process you can see why the presence of such a singularity will cause trouble! You have to include the singular points among the ends of the intervals in which you subdivide your domain of integration, otherwise your integrals can't be computed accurately. This is what Mathematica is warning you about. That is why the correct code, which I did send you, is: NIntegrate[f[x], {x, -1, 13, 17}] This tells Mathematica to include 13 among the points of the subdivision and gives you accurate result with no error. You seem to have completely ignired that part of my message! Next to your serious problem.: > Here is a better example of what I want to do. Please run it and you will see > what I am talking about. My program is more complex so I do not get an answer. > > PeakUp[n_] := n + 100; > G = 3; > DOSup[n_, Energy_] := (1 - ((Energy - PeakUp[n])/G)^2)^(.5) /; > (PeakUp[n] - G <= Energy && Energy <= PeakUp[n] + G) > DOSup[n_, Energy_] := 0 /; (Energy <= PeakUp[n] - G) > DOSup[n_, Energy_] := 0 /; (PeakUp[n] + G <= Energy); > > NumberDensity[n_, Evar_] := Integrate[DOSup[n, Energy], {Energy, 0, Evar}]; > Nlevel = 7; > TotalNumberDensity[Efinal_] := Sum[NumberDensity[n, Efinal], {n, 0, Nlevel}]; > Plot[TotalNumberDensity[Efinal], {Efinal, 0, 500}] > Basically, you have the same problem as earlier, but there are some programming difficulties. First of all, you need to use NIntegrate, not Integrate. However, the trouble with NIntegrate is that it can't integrate the zero function without producing an error message: In[1]:= NIntegrate[0, {t, 0, 5}] NIntegrate::"ploss": "Numerical integration stopping due to loss of \ precision. Achieved neither the requested PrecisionGoal nor AccuracyGoal; \ suspect one of the following: highly oscillatory integrand or the true value \ of the integral is 0. If your integrand is oscillatory try using the option \ Method->Oscillatory in NIntegrate." Out[1]= 0. I think this is something to do with the fact that the Precison of 0. is 0 (but I am not very expert at this numerical stuff). Anyway, in your final computation you are often integrating the 0 function so you would be getting a lot of this sort of messages unless you do something to prevent this. The other problem is just as the above: you have to tell Mathematica about the singularities. Here is my rather hastily written code that does it in your case (I am pretty sure it can be improved), with a brief commentary: In[2]:= PeakUp[n_] := n + 100; In[3]:= G = 3; In[4]:= DOSup[n_, Energy_] := (1 - ((Energy - PeakUp[n])/G)^2)^(.5) /; (PeakUp[n] - G <= Energy <= PeakUp[n] + G) In[5]:= DOSup[n_, Energy_] := 0 /; (Energy <= PeakUp[n] - G) In[6]:= DOSup[n_, Energy_] := 0 /; (PeakUp[n] + G <= Energy) So far as in your example. Now I introduce a new function: In[7]:= range[n_, Evar_] := Sort[Select[{0, PeakUp[n] - G, PeakUp[n] + G, Evar}, (0 <= # <= Evar) &]] This just makes a list {0,first singular point, second singular point, final point} We will want to consider only cases where there is at least one singular point between 0 and the final point (Evar). In[8]:= NumberDensity[n_, Evar_] := If[Length[range[n, Evar]] >= 3, NIntegrate[Evaluate[DOSup[n, Energy]], Evaluate[Prepend[range[n, Evar], Energy]]], 0]; Here is the changed definition of NumberDensity. First of all we only consider the case when the list range[n,Evar] has at least three members, because otherwise we would be NIntegrating the zero function and get all these error messages. Anyway, we know that the integral of the zero function is 0. Note also the presence of Evaluate. This is because NIntegrate has the attribute HoldAll, so it does not evaluate its arguments. Finally note that I included the singular points in the range of integration. The rest is as in your own code: In[9]:= Nlevel = 7; In[10]:= TotalNumberDensity[Efinal_] := Sum[NumberDensity[n, Efinal], {n, 0, Nlevel}]; In[10]:= Plot[TotalNumberDensity[Efinal], {Efinal, 0, 500}] here I get the expected graph with no error messages. > > > > Andrzej Kozlowski wrote: > >> You can't expect Mathematica to work if you feed it a command which vilates >> the rules of the Mathematica language (your second example). What you wrote >> is not a Mathematica program. The correct form is: >> >> f[x_] := If[Abs[x - 3] < 10, Sqrt[1 - ( (x - 3)/10)^2] , 0] >> >> However, I do not understand why you claim that Mathematica "goes nuts" in >> the first case, which is more or less correct. (Actually a bit weird: why >> do you use contradictory conditions over an overlaping range of x? >> Presumably you meant the first to be Abs[x-3]<10 and the second >> Abs[x-3]>=10. In any case, it does not metter for your second condition is >> ignored where it contradicts the first, e.g. when x=7.) What I get, even >> with your original code, is this: >> >> In[56]:= >> Clear[f] >> >> In[57]:= >> f[x_] := Sqrt[ 1 - ( (x - 3)/10)^2 ] /; >> Abs[x - 3] < 10 >> >> In[58]:= >> f[x_] := 0 /; >> Abs[x - 3] > 3 ; >> >> In[59]:= >> NIntegrate[f[x], {x, -1, 17}] >> >>> From In[59]:= >> NIntegrate::"ncvb": "NIntegrate failed to converge to prescribed accuracy >> after 7 recursive bisections in x near x = 12.9921875`." >> >> Out[59]= >> 11.7447 >> >> All this says is that Mathematica had trouble near the singularity. The >> answer it gets is correct up to the first three digits after the decimal >> point, as you can see from: >> >> In[60]:= >> NIntegrate[f[x], {x, -1, 13, 17}] >> Out[60]= >> 11.7446 >> >> This latter form of NIntegrate tells Mathematica about the singularity at 13 >> and avoids the trouble. >> In any case it is not at all surprising that a computer program finds it >> difficult to cope with singularites, so do people. If you call this "going >> nuts" I wonder what your notion of sanity is :). >> >> -- >> Andrzej Kozlowski >> Toyama International University, JAPAN >> >> For Mathematica related links and resources try: >> <http://www.sstreams.com/Mathematica/> >> >> on 6/27/00 1:51 PM, Viorel Ontalus at vio2 at mail.lehigh.edu wrote: >> >>> It seems I got into an area where Mathematica has some problems, and I >>> hope somebody can give me a hint on how to go around these problems. >>> >>> 1. I am trying to define a piecewise function and do some computations >>> with it. When I integrate mathematica does not behave. Here is an simple >>> example you can run and see what I am talking about >>> >>> Clear[f,x} >>> f[x_] := Sqrt[ 1- ( (x-3)/10)^2 ] /; >>> Abs[x-3]<10 >>> f[x_] := 0 /; >>> Abs[x-3]>3 ; (*this is a very simple piecewise function but one >>> must be sure the Sqrt is from a positive # ) >>> >>> NIntegrate[f[x],{x,-1,17}] (* Here Mathematica goes nuts !!!*) >>> >>> Of course it gives an answer but if your program is more complex, then >>> you never get an answer !! >>> ( I tried to make the upper limit a variable etc !!) >>> Does anybody know how to avoid the error , or non convergence messages I >>> get !! >>> >>> >>> >>> 2. For fun I tried the only reference from the book on piecewise >>> functions: >>> If[Abs[x-3]<10, f[x_]:=Sqrt[1- ( (x-3)/10)^2 ] , f[x_]:=0 ] >>> >>> This definition does not work !! >>> >>> >>> >>> >>> >>> >>> >>> >>> >