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!

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
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!

> 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 !!
>>>
>>>
>>>
>>>
>>>
>>>
>>>
>>>
>>>
>

```

• Prev by Date: Re: Mathematica for High School Students
• Next by Date: RE: How to operate on strictly numerical functions ?
• Previous by thread: Re: Piecewise functions definition and usage
• Next by thread: Unix Sound Problems