MathGroup Archive 2010

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: How to interpret this integral?

  • To: mathgroup at smc.vnet.net
  • Subject: [mg112614] Re: How to interpret this integral?
  • From: Leonid Shifrin <lshifr at gmail.com>
  • Date: Thu, 23 Sep 2010 04:21:13 -0400 (EDT)

Julian,

Mathematica does not automatically handle most of singularities in the sense
you would like it to.
In your case, the non-zero aW changes the answer qualitatively. But for
Mathematica, it is not at
all obvious that the aW==0 should be "considered" by it on a separate basis.
Who knows, may be
 this is just a part of the more general expression, where the singularity
will be cancelled etc. Using
your formula at aW ==0 meant taking the limit aW->0 (and my procedure was
just a trick to make it
easier), but  Mathematica does not automatically do this kind of things
since it may not always be
 a right thing to do. What you, as a Mathematica programmer, can do, is to
define your function in a
way which reflects special points:

myFunction[0,aV_,v0_,th0_,w0_,Ts_]:==<specific answer>
myFunction[aW_,aV_,v0_,th0_,w0_,Ts_]:==<generic answer>

In this case, the more specific definitions will automatically be tried
before the more general ones.
If you have an idea about possible singular points, the process of giving
such specialized
definitions can be somewhat automated, and that is exactly the task of
applied Mathematica programmer
(a scientist using Mathematica who is able to use the domain-specific
knowledge in in implementing such automation), but you still must perform
such an analysis yourself. I actually doubt that there can be any totally
general approach to situations like this. IMO, routine stuff is for the
machine, but singularities are for humans )

Sorry if I wasn't exactly helpful.

Best,
Leonid


On Wed, Sep 22, 2010 at 7:55 PM, Julian Stoev / =D0=AE=D0=BB=D0=B8=D0=B0=D0=
=BD =D0=A1=D1=82=D0=BE=D0=B5=D0=B2 <
julian.stoev at gmail.com> wrote:

> Thanks Leonid,
>
> I also had your first idea, but I am a little suspicious if this is a
> good general approach.
> If aW is a special point, shouldn't Mathematica provide a solution in
> terms of If[] functions?
> There must be some reason why it is not doing so.
>
> And Mathematica is not the only software giving me an answer with this
> form...
>
> Also I was not able to force it to take into account any assumptions,
> like real variables, real Sqr,...
> There must be some regular non-manual way of getting a solution?
>
> Best!
>
> --JS
>
> On Wed, Sep 22, 2010 at 5:29 PM, Leonid Shifrin <lshifr at gmail.com> wrote:
> > Hi Julian,
> > Here is your expression:
> > In[137]:== expr == Integrate[(aV*t+v0)*Cos[(aW*t^2)/2+th0+t*w0],{t,0,Ts=
}]
> > Out[137]==
> > (1/(aW^(3/2)))(Sqrt[\[Pi]] (-aW v0+aV w0) Cos[th0-w0^2/(2 aW)]
> > FresnelC[w0/(Sqrt[aW] Sqrt[\[Pi]])]+Sqrt[\[Pi]] (aW v0-aV w0)
> > Cos[th0-w0^2/(2 aW)] FresnelC[(aW Ts+w0)/(Sqrt[aW] Sqrt[\[Pi]])]+aV
> Sqrt[aW]
> > (-Sin[th0]+Sin[th0+(aW Ts^2)/2+Ts w0])+Sqrt[\[Pi]] (aW v0-aV w0)
> > (FresnelS[w0/(Sqrt[aW] Sqrt[\[Pi]])]-FresnelS[(aW Ts+w0)/(Sqrt[aW]
> > Sqrt[\[Pi]])]) Sin[th0-w0^2/(2 aW)])
> > If you want to treat the point aW == 0 specially, the simplest is to pu=
t
> it
> > to zero right when the integral
> > is computed:
> > In[138]:== exprAt0 == Integrate[(aV*t+v0)*Cos[th0+t*w0],{t,0,Ts}]
> > Out[138]== (-aV Cos[th0]+aV Cos[th0+Ts w0]-v0 w0 Sin[th0]+(aV Ts+v0) w0
> > Sin[th0+Ts w0])/w0^2
> > Alternatively, you can extract it from the full expression, but it st a
> > little more involved. One way
> > is to change variable as t == 1/aW, then expand in t around t==Infinity=
,
> then
> > take the limit t->Infinity of the
> > resulting expression. One subtlety here that you will need to expand to
> the
> > second order in t, to
> > get the correct result:
> > In[140]:== exprAt0Alt ==
> > Limit[Normal@Series[expr/.aW->1/t,{t,Infinity,2}],t->Infinity]
> > Out[140]== (-aV Cos[th0]+aV Cos[th0+Ts w0]+w0 (-v0 Sin[th0]+aV Ts
> Sin[th0+Ts
> > w0]+v0 Sin[th0+Ts w0]))/w0^2
> > You can check that the results are the same:
> > In[141]:== FullSimplify[exprAt0-exprAt0Alt]
> > Out[141]== 0
> > Hope this helps.
> > Regards,
> > Leonid
> >
> >
> > On Wed, Sep 22, 2010 at 9:57 AM, Julian <julian.stoev at gmail.com> wrote:
> >>
> >> Hello All,
> >>
> >> After long interruption with symbolic computing, I am struggling how
> >> to make a useful result out of this integral.
> >>
> >> Integrate[(aV*t + v0)*Cos[(aW*t^2)/2 + th0 + t*w0, {t, 0, Ts}]
> >>
> >> It comes from the differential equations describing the motion of a
> >> wheeled robot. aV and aW are accelerations and v0, w0, th0 are initial
> >> conditions. The numerical solution clearly exists for different real
> >> accelerations, including positive, negative and zero.
> >>
> >> However the symbolic solution of Mathematica is:
> >> (Sqrt[Pi]*(-(aW*v0) + aV*w0)*Cos[th0 - w0^2/(2*aW)]*
> >>   FresnelC[w0/(Sqrt[aW]*Sqrt[Pi])] + Sqrt[Pi]*(aW*v0 - aV*w0)*
> >>   Cos[th0 - w0^2/(2*aW)]*FresnelC[(aW*Ts + w0)/(Sqrt[aW]*Sqrt[Pi])]
> >> +
> >>  aV*Sqrt[aW]*(-Sin[th0] + Sin[th0 + (aW*Ts^2)/2 + Ts*w0]) +
> >>  Sqrt[Pi]*(aW*v0 - aV*w0)*(FresnelS[w0/(Sqrt[aW]*Sqrt[Pi])] -
> >>    FresnelS[(aW*Ts + w0)/(Sqrt[aW]*Sqrt[Pi])])*Sin[th0 - w0^2/
> >> (2*aW)])/
> >>  aW^(3/2)
> >>
> >> Note that aW can be found inside a Sqrt function and also in the
> >> denominator. While I can see that FresnelS[Infinity] is well defined,
> >> so aW====0 should not be a real problem, it is  still very problematic
> >> how to use this result at this particular point. The case of negative
> >> aW is also interesting, because it will require a complex FresnelS and
> >> FresnelC. I somehow managed to remove this problem using
> >> ComplexExpand, but I am not sure this is the good solution.
> >>
> >> I tried to give Assumptions -> Element[{aW, aV, w0, v0, Ts, th0},
> >> Reals] to the integral, but there is no change in the solution.
> >>
> >> While I understand that the solution of Mathematica is correct in the
> >> strict mathematical sense, I have to use the result to generate a C-
> >> code, which will be evaluated numerically and the solution I have now
> >> is not working for this.
> >>
> >> Can some experienced user give a good advice how use get a solution of
> >> the problem?
> >>
> >> Thank you in advance!
> >>
> >> --JS
> >>
> >
> >
>
>
>
> --
> Julian Stoev, PhD.
> Control Researcher
>

--0016e6d7ef665f64fa0490db77b2
Content-Type: text/html; charset="UTF-8"
Content-Transfer-Encoding: quoted-printable
X-Sun-Content-Length: 7720

Julian,=C2 <div><br></div><div>Mathematica does not automatically handle =
most of singularities in the sense you would like it to.=C2 </div><div>In=
 your case, the non-zero=C2 aW changes the answer qualitatively. But for =
Mathematica, it is not at=C2 </div>
<div>all obvious that the aW=0=C2 should be &quot;considered&quot; by i=
t on a separate basis. Who knows, may be</div><div>=C2 this is just a par=
t of=C2 the more general expression, where the singularity will be cancel=
led etc. Using=C2 </div>
<div>your formula at=C2 aW =0 meant taking the limit aW-&gt;0 (and my p=
rocedure was just a trick to make it=C2 </div><div>easier), but =C2 Mat=
hematica does not automatically do this kind of=C2 things since it may no=
t always be</div><div>
=C2 a right=C2 thing to do.=C2 What you, as a Mathematica=C2 progra=
mmer, can do, is to define your function in a=C2 </div><div>way which ref=
lects special points:</div><div><br></div><div>myFunction[0,aV_,v0_,th0_,w0=
_,Ts_]:=&lt;specific answer&gt;</div>
<div>myFunction[aW_,aV_,v0_,th0_,w0_,Ts_]:=&lt;generic answer&gt;</div><d=
iv><br></div><div>In this case, the more specific definitions will automati=
cally be tried before the more general ones.</div><div>If you have an idea =
about possible singular points, the process of giving such specialized</div=
>
<div>definitions can be somewhat automated, and that is exactly the task of=
 applied Mathematica programmer</div><div>(a scientist using Mathematica wh=
o is able to use the domain-specific knowledge in in implementing such auto=
mation), but you still must perform such an analysis yourself. I=C2 actua=
lly doubt that there can be any totally general approach to situations like=
 this. IMO,=C2 routine stuff is for the machine, but singularities are fo=
r humans )</div>
<div><br></div><div>Sorry if I wasn&#39;t exactly helpful.</div><div><br></=
div><div>Best,</div><div>Leonid</div><div><br></div><div><br><div class="=
gmail_quote">On Wed, Sep 22, 2010 at 7:55 PM, Julian Stoev /  =D0=AE=D0=BB=
=D0=B8=D0=B0=D0=BD =D0=A1=D1=82=D0=BE=D0=B5=D0=B2 <span dir="ltr">&lt;<a =
href="mailto:julian.stoev at gmail.com">julian.stoev at gmail.com</a>&gt;</span=
> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1p=
x #ccc solid;padding-left:1ex;">Thanks Leonid,<br>
<br>
I also had your first idea, but I am a little suspicious if this is a<br>
good general approach.<br>
If aW is a special point, shouldn&#39;t Mathematica provide a solution in<b=
r>
terms of If[] functions?<br>
There must be some reason why it is not doing so.<br>
<br>
And Mathematica is not the only software giving me an answer with this form=
...<br>
<br>
Also I was not able to force it to take into account any assumptions,<br>
like real variables, real Sqr,...<br>
There must be some regular non-manual way of getting a solution?<br>
<br>
Best!<br>
<br>
--JS<br>
<div><div></div><div class="h5"><br>
On Wed, Sep 22, 2010 at 5:29 PM, Leonid Shifrin &lt;<a href="mailto:lshif=
r at gmail.com">lshifr at gmail.com</a>&gt; wrote:<br>
&gt; Hi Julian,<br>
&gt; Here is your expression:<br>
&gt; In[137]:= expr = Integrate[(aV*t+v0)*Cos[(aW*t^2)/2+th0+t*w0],{t,0=
,Ts}]<br>
&gt; Out[137]=<br>
&gt; (1/(aW^(3/2)))(Sqrt[\[Pi]] (-aW v0+aV w0) Cos[th0-w0^2/(2 aW)]<br>
&gt; FresnelC[w0/(Sqrt[aW] Sqrt[\[Pi]])]+Sqrt[\[Pi]] (aW v0-aV w0)<br>
&gt; Cos[th0-w0^2/(2 aW)] FresnelC[(aW Ts+w0)/(Sqrt[aW] Sqrt[\[Pi]])]+aV Sq=
rt[aW]<br>
&gt; (-Sin[th0]+Sin[th0+(aW Ts^2)/2+Ts w0])+Sqrt[\[Pi]] (aW v0-aV w0)<br>
&gt; (FresnelS[w0/(Sqrt[aW] Sqrt[\[Pi]])]-FresnelS[(aW Ts+w0)/(Sqrt[aW]<br>
&gt; Sqrt[\[Pi]])]) Sin[th0-w0^2/(2 aW)])<br>
&gt; If you want to treat the point aW = 0 specially, the simplest is to =
put it<br>
&gt; to zero right when the integral<br>
&gt; is computed:<br>
&gt; In[138]:= exprAt0 = Integrate[(aV*t+v0)*Cos[th0+t*w0],{t,0,Ts}]<br=
>
&gt; Out[138]= (-aV Cos[th0]+aV Cos[th0+Ts w0]-v0 w0 Sin[th0]+(aV Ts+v0) =
w0<br>
&gt; Sin[th0+Ts w0])/w0^2<br>
&gt; Alternatively, you can extract it from the full expression, but it st =
a<br>
&gt; little more involved. One way<br>
&gt; is to change variable as t = 1/aW, then expand in t around t=Infin=
ity, then<br>
&gt; take the limit t-&gt;Infinity of the<br>
&gt; resulting expression. One subtlety here that you will need to expand t=
o the<br>
&gt; second order in t, to<br>
&gt; get the correct result:<br>
&gt; In[140]:= exprAt0Alt =<br>
&gt; Limit[Normal@Series[expr/.aW-&gt;1/t,{t,Infinity,2}],t-&gt;Infinity]<b=
r>
&gt; Out[140]= (-aV Cos[th0]+aV Cos[th0+Ts w0]+w0 (-v0 Sin[th0]+aV Ts Sin=
[th0+Ts<br>
&gt; w0]+v0 Sin[th0+Ts w0]))/w0^2<br>
&gt; You can check that the results are the same:<br>
&gt; In[141]:= FullSimplify[exprAt0-exprAt0Alt]<br>
&gt; Out[141]= 0<br>
&gt; Hope this helps.<br>
&gt; Regards,<br>
&gt; Leonid<br>
&gt;<br>
&gt;<br>
&gt; On Wed, Sep 22, 2010 at 9:57 AM, Julian &lt;<a href="mailto:julian.s=
toev at gmail.com">julian.stoev at gmail.com</a>&gt; wrote:<br>
&gt;&gt;<br>
&gt;&gt; Hello All,<br>
&gt;&gt;<br>
&gt;&gt; After long interruption with symbolic computing, I am struggling h=
ow<br>
&gt;&gt; to make a useful result out of this integral.<br>
&gt;&gt;<br>
&gt;&gt; Integrate[(aV*t + v0)*Cos[(aW*t^2)/2 + th0 + t*w0, {t, 0, Ts}]<br>
&gt;&gt;<br>
&gt;&gt; It comes from the differential equations describing the motion of =
a<br>
&gt;&gt; wheeled robot. aV and aW are accelerations and v0, w0, th0 are ini=
tial<br>
&gt;&gt; conditions. The numerical solution clearly exists for different re=
al<br>
&gt;&gt; accelerations, including positive, negative and zero.<br>
&gt;&gt;<br>
&gt;&gt; However the symbolic solution of Mathematica is:<br>
&gt;&gt; (Sqrt[Pi]*(-(aW*v0) + aV*w0)*Cos[th0 - w0^2/(2*aW)]*<br>
&gt;&gt; =C2  FresnelC[w0/(Sqrt[aW]*Sqrt[Pi])] + Sqrt[Pi]*(aW*v0 - aV*w0)=
*<br>
&gt;&gt; =C2  Cos[th0 - w0^2/(2*aW)]*FresnelC[(aW*Ts + w0)/(Sqrt[aW]*Sqrt=
[Pi])]<br>
&gt;&gt; +<br>
&gt;&gt; =C2 aV*Sqrt[aW]*(-Sin[th0] + Sin[th0 + (aW*Ts^2)/2 + Ts*w0]) +<b=
r>
&gt;&gt; =C2 Sqrt[Pi]*(aW*v0 - aV*w0)*(FresnelS[w0/(Sqrt[aW]*Sqrt[Pi])] -=
<br>
&gt;&gt; =C2  =C2 FresnelS[(aW*Ts + w0)/(Sqrt[aW]*Sqrt[Pi])])*Sin[th0 -=
 w0^2/<br>
&gt;&gt; (2*aW)])/<br>
&gt;&gt; =C2 aW^(3/2)<br>
&gt;&gt;<br>
&gt;&gt; Note that aW can be found inside a Sqrt function and also in the<b=
r>
&gt;&gt; denominator. While I can see that FresnelS[Infinity] is well defin=
ed,<br>
&gt;&gt; so aW==0 should not be a real problem, it is =C2 still very =
problematic<br>
&gt;&gt; how to use this result at this particular point. The case of negat=
ive<br>
&gt;&gt; aW is also interesting, because it will require a complex FresnelS=
 and<br>
&gt;&gt; FresnelC. I somehow managed to remove this problem using<br>
&gt;&gt; ComplexExpand, but I am not sure this is the good solution.<br>
&gt;&gt;<br>
&gt;&gt; I tried to give Assumptions -&gt; Element[{aW, aV, w0, v0, Ts, th0=
},<br>
&gt;&gt; Reals] to the integral, but there is no change in the solution.<br=
>
&gt;&gt;<br>
&gt;&gt; While I understand that the solution of Mathematica is correct in =
the<br>
&gt;&gt; strict mathematical sense, I have to use the result to generate a =
C-<br>
&gt;&gt; code, which will be evaluated numerically and the solution I have =
now<br>
&gt;&gt; is not working for this.<br>
&gt;&gt;<br>
&gt;&gt; Can some experienced user give a good advice how use get a solutio=
n of<br>
&gt;&gt; the problem?<br>
&gt;&gt;<br>
&gt;&gt; Thank you in advance!<br>
&gt;&gt;<br>
&gt;&gt; --JS<br>
&gt;&gt;<br>
&gt;<br>
&gt;<br>
<br>
<br>
<br>
</div></div><font color="#888888">--<br>
Julian Stoev, PhD.<br>
Control Researcher<br>
</font></blockquote></div><br></div>

--0016e6d7ef665f64fa0490db77b2--


  • Prev by Date: Re: How to interpret this integral?
  • Next by Date: ContourPlot under ListPlot3D
  • Previous by thread: Re: How to interpret this integral?
  • Next by thread: Calendar Recipe?