Mathematica 9 is now available
Services & Resources / Wolfram Forums / MathGroup Archive
-----

MathGroup Archive 2012

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

Search the Archive

Re: DSolve for a real function

  • To: mathgroup at smc.vnet.net
  • Subject: [mg127992] Re: DSolve for a real function
  • From: Bob Hanlon <hanlonr357 at gmail.com>
  • Date: Fri, 7 Sep 2012 04:59:10 -0400 (EDT)
  • Delivered-to: l-mathgroup@mail-archive0.wolfram.com
  • Delivered-to: l-mathgroup@wolfram.com
  • Delivered-to: mathgroup-newout@smc.vnet.net
  • Delivered-to: mathgroup-newsend@smc.vnet.net
  • References: <k1k8np$90f$1@smc.vnet.net>

The solutions are real in the region of interest, i.e., all of the
complex components cancel out even if these cancellations are not
explicit in the representation. A trivial example of a real expression
that on the surface may look complex is  ((a+b*I) + (a-b*I)). In this
case the cancellation of the imaginaary parts can be easily done. Your
expression is sufficiently intricate that it is exceedingly difficult
if not impossible to find a representation that explicitly cancels out
all of the complex parts in a symbolic form.

sol1 = y[x] /.
   DSolve[{y''''[x] + a y[x] == 0, y''[-c] == ic0, y''[c] == ic=
0,
      y'''[-c] == ic1, y'''[c] == -ic1}, y[x], x][[1]];

sol2 = Simplify[sol1, {a > 0, Element[{ic0, ic1, c, x}, Reals]}];

sol3 = ComplexExpand[sol2];

sol4 = Simplify[sol3, {a > 0, Element[{ic0, ic1, c, x}, Reals]},
   TimeConstraint -> 1000];

For at least some exact values of the constants, the cancellations can
be done for some of these representations resulting in explicitly real
representations

sol3x = FullSimplify[sol3 /.
   {ic0 -> 2, ic1 -> 2, c -> 1, a -> 1}]

((-1 + Sqrt[2] + (1 + Sqrt[2])*
               E^(Sqrt[2]*(1 + x)))*
          Cos[(-1 + x)/Sqrt[2]] +
        ((1 + Sqrt[2])*E^Sqrt[2] +
             (-1 + Sqrt[2])*
               E^(Sqrt[2]*x))*
          Cos[(1 + x)/Sqrt[2]] +
        (1 + E^(Sqrt[2]*(1 + x)))*
          Sin[(-1 + x)/Sqrt[2]] -
        (E^Sqrt[2] + E^(Sqrt[2]*x))*
          Sin[(1 + x)/Sqrt[2]])/
     E^((1 + x)/Sqrt[2])/
   (Sin[Sqrt[2]] + Sinh[Sqrt[2]])

Plot[sol3x, {x, 0, 9}]

For specific numerical evaluations of the expressions--particularly
for large values of x, using machine precision numbers can result in
very large numerical noise which precludes the proper calculation of
the imaginary components.

{sol1, sol2, sol3, sol4} /. {
  ic0 -> RandomReal[{1, 10}],
  ic1 -> RandomReal[{1, 10}],
  c -> RandomReal[{1, 10}],
  a -> RandomReal[{1, 10}],
  x -> 100}

{1.2276490765332772*^39 +
     2.8031967442314214*^25*I,
   1.2276490765333458*^39 +
     2.0702854660900525*^25*I,
   1.2276490765333507*^39 +
     1.8728905471061013*^25*I,
   1.2276490765333449*^39 +
     1.735742182536485*^25*I}

Using much higher precision can cause the cancellation of the
imaginary components to work effectively resulting in negligible
imaginary components

{sol1, sol2, sol3,
  sol4} /. {ic0 -> RandomReal[{1, 10}, WorkingPrecision -> 75],
  ic1 -> RandomReal[{1, 10}, WorkingPrecision -> 75],
  c -> RandomReal[{1, 10}, WorkingPrecision -> 75],
  a -> RandomReal[{1, 10}, WorkingPrecision -> 75],
  x -> 100}

{3.5466406369437226924022894337705\
   5420786716057206415883532647413844\
   97269532445274677294054176313`72.0\
   723740502742*^36 +
     0``35.52384127280971*I,
   3.546640636943722692402289433770\
   5542078671605720641588353264741384\
    4972695324452746772940541761`72.62\
   40529160406*^36 +
     0``36.11439423733123*I,
   3.546640636943722692402289433770\
   5542078671605720641588353264741384\
    497269532445274677294054176096`72.\
   52533544339316*^36 +
     0``35.98812103665155*I,
   3.546640636943722692402289433770\
   5542078671605720641588353264741384\
    497269532445274677294054176072`72.\
   79515994843076*^36 +
     0``36.3762727216666*I}

The residual imaginary artifacts are removed with Chop

% // Chop

{3.5466406369437226924022894337705\
  5420786716057206415883532647413844\
  97269532445274677294054176313`72.0\
  723740502742*^36,
 3.54664063694372\
  2692402289433770554207867160572064\
  1588353264741384497269532445274677\
  2940541761`72.6240529160406*^36,
 3\
  .546640636943722692402289433770554\
  2078671605720641588353264741384497\
  269532445274677294054176096`72.525\
  33544339316*^36,
 3.546640636943722\
  6924022894337705542078671605720641\
  5883532647413844972695324452746772\
  94054176072`72.79515994843076*^36}


Bob Hanlon


On Thu, Sep 6, 2012 at 4:26 AM, Andreas Talmon l'Arm=E9e
<talmon at fsm.tu-darmstadt.de> wrote:
> I also posted  my problem at
>
> http://mathematica.stackexchange.com/questions/10253/dsolve-gives-complex=
-function-although-the-solution-is-a-real-one
>
> and got the answer to do something like:
>
> sol1 = DSolve[{y''''[x] + a y[x] == 0, y''[-c] == ic0, y''[c] =
== ic0 ,
>
>      y'''[-c] == ic1 , y'''[c] == - ic1
>     }, y[x], x];
>
> sol2 = Simplify[
>   sol1, {a > 0, ic0 \[Element] Reals, ic1 \[Element] Reals,
>    c \[Element] Reals, x \[Element] Reals}]
>
> sol3 = ComplexExpand[y[x] /. sol2[[1, 1]]]
>
>
> sol4 = Simplify[
>   sol3, {a > 0, ic0 \[Element] Reals, ic1 \[Element] Reals,
>    c \[Element] Reals, x \[Element] Reals}, TimeConstraint -> 1000]
>
> But the solution is still complex.
> In sol3 I have the solution decomposed in an Real and Imaginary part. Any
> ideas how to extract only the Real part of the solution?
>
> Re[sol3] does not work
> Simplify[Re[sol3]] also does not work
>
> Andreas Talmon l'Arm=E9e
>
>
>
>
> On 09/03/2012 06:54 PM, Bob Hanlon wrote:
>>
>> In numerical computations, the specific computations and their order
>> of computation determine the accumulation of the numerical noise. You
>> generally reduce the accumulation of numerical noise by increasing the
>> working precision.
>>
>> For example in your t2 rather than using machine precision specify
>> extended precision (say 55 places)
>>
>> t2 = Table[
>>     y4[x] /. {c -> 1, ic0 -> 2, ic1 -> 3, a -> 4}, {x, 0, 100, 10}] //N[=
#,
>> 55] &
>>
>> Similarly, when plotting y4
>>
>> pRe2 = Plot[
>>     Chop[y4[x]] /. {c -> 1, ic0 -> 2, ic1 -> 3, a -> 4}, {x, 0, 100},
>>     WorkingPrecision -> 55, PlotStyle -> Red];
>>
>> pIm2 = Plot[Im[y4[x] /. {c -> 1, ic0 -> 2, ic1 -> 3, a -> 4}],
>>    {x, 0, 100}, WorkingPrecision -> 55,
>>    Frame -> True, Axes -> False]
>>
>>
>> Bob Hanlon
>>
>>
>> On Mon, Sep 3, 2012 at 10:22 AM, Andreas Talmon l'Arm=E9e
>> <talmon at fsm.tu-darmstadt.de> wrote:
>>>
>>> You are completely right but I think there should be no real part in th=
e
>>> solution.
>>>
>>> I did the solution of the differential equation on paper and got the te=
rm
>>> for y1[x].
>>>
>>> Then I determined the constants C[5]...C[8] with the boundary condition=
s
>>> with Solve[] in  Mathematica. There I got an unusual behaviour.
>>>
>>> When using Simplify on the sol1 could generate a soltuion with only a
>>> real
>>> part.
>>> When using FullSimplify I get the imaginary part which is very small fo=
r
>>> small x but for greater x also the imaginary part starts to oscillate
>>> with a
>>> great amplitude.
>>>
>>> I do not understand how a solution can be Real or Complex depending on
>>> the
>>> Simplifications scheme I use.
>>>
>>> Shouldn't there be one solution to a problem with specified boundary
>>> conditions.?
>>>
>>> The notebook where I did all this is on:
>>> http://dl.dropbox.com/u/4920002/DGL_4th_Order_with_own_solution.nb
>>>
>>>
>>> Thanks for your great help,
>>>
>>> Andreas.
>>>
>>>
>>>
>>>
>>> On 09/03/2012 03:22 PM, Bob Hanlon wrote:
>>>>
>>>> Clear[sol];
>>>>
>>>> $Assumptions = Element[{a, c, ic0, ic1}, Reals];
>>>>
>>>> sol[a_, x_] =
>>>>    y[x] /. DSolve[{y''''[x] + a y[x] == 0, y''[-c] == ic0, y''=
[c] ==
>>>> ic0,
>>>>         y'''[-c] == ic1, y'''[c] == -ic1}, y[x], x][[1]] // Fu=
llSimplify
>>>>
>>>> ((1 + I)*a^(1/4)*
>>>>             (I*E^(I*Sqrt[2]*a^(1/4)*c) -
>>>>                I*E^(Sqrt[2]*a^(1/4)*c) +
>>>>                E^(Sqrt[2]*a^(1/4)*
>>>>                     ((1 + I)*c + I*x)) -
>>>>                E^(I*Sqrt[2]*a^(1/4)*x) -
>>>>                E^(Sqrt[2]*a^(1/4)*x) +
>>>>                E^(Sqrt[2]*a^(1/4)*
>>>>                     ((1 + I)*c + x)) +
>>>>                I*E^(Sqrt[2]*a^(1/4)*
>>>>                       (I*c + (1 + I)*x)) -
>>>>                I*E^(Sqrt[2]*a^(1/4)*
>>>>                       (c + (1 + I)*x)))*ic0 +
>>>>           Sqrt[2]*(E^(I*Sqrt[2]*a^(1/4)*
>>>>                     c) + E^(Sqrt[2]*a^(1/4)*
>>>>                     c) + E^(Sqrt[2]*a^(1/4)*
>>>>                     ((1 + I)*c + I*x)) +
>>>>                E^(I*Sqrt[2]*a^(1/4)*x) +
>>>>                E^(Sqrt[2]*a^(1/4)*x) +
>>>>                E^(Sqrt[2]*a^(1/4)*
>>>>                     ((1 + I)*c + x)) +
>>>>                E^(Sqrt[2]*a^(1/4)*
>>>>                     (I*c + (1 + I)*x)) +
>>>>                E^(Sqrt[2]*a^(1/4)*
>>>>                     (c + (1 + I)*x)))*ic1)/
>>>>        E^((-1)^(1/4)*a^(1/4)*(c + x))/
>>>>      (4*a^(3/4)*(Sin[Sqrt[2]*a^(1/4)*
>>>>                c] + Sinh[Sqrt[2]*a^(1/4)*
>>>>                c]))
>>>>
>>>> Table[sol[a, x] /. {
>>>>      c -> RandomReal[{0, 10}, WorkingPrecision -> 20],
>>>>      ic0 -> RandomReal[{0, 10}, WorkingPrecision -> 20],
>>>>      ic1 -> RandomReal[{0, 10}, WorkingPrecision -> 20]},
>>>>    {a, 1, 5}, {x, 0, 4}]
>>>>
>>>> {{0.01540650723811880994 + 0.*10^-21 I, -0.0003703423761496348 +
>>>>      0.*10^-20 I, -0.1278269715054133646 + 0.*10^-20 I,
>>>>     0.0392723836429585004 + 0.*10^-20 I,
>>>>     13.10841701158527505 + 0.*10^-18 I}, {-0.00286328430472954469 +
>>>>      0.*10^-21 I, -1.985775360248369253 + 0.*10^-19 I,
>>>> -0.311655087881410281 +
>>>>      0.*10^-19 I, -0.718055887186593765 + 0.*10^-19 I,
>>>>     21.42291506291193811 + 0.*10^-18 I}, {0.419394450195258310 +
>>>> 0.*10^-19
>>>> I,
>>>>     0.01666534839706064001 + 0.*10^-21 I,
>>>>     0.528881084742103412 + 0.*10^-19 I, -0.060937975960727530 +
>>>> 0.*10^-19
>>>> I,
>>>>     0.0086995512380544778 + 0.*10^-20 I}, {0.01847749159520439918 +
>>>> 0.*10^-21 I,
>>>>      0.00517593237021928535 + 0.*10^-21 I,
>>>>     0.0134376958759010127 + 0.*10^-20 I, -0.0679397801010012506 +
>>>>      0.*10^-20 I, -5.21491473427110550 +
>>>>      0.*10^-18 I}, {-0.00168718272676310080 +
>>>>      0.*10^-21 I, -0.00137069152420374693 + 0.*10^-21 I,
>>>>     3.051457492528748582 + 0.*10^-19 I,
>>>>     0.00026536991493989622 + 0.*10^-21 I, -59.38147365750750826 +
>>>> 0.*10^-18
>>>> I}}
>>>>
>>>> The imaginary parts cancel out; the residual imaginary parts are just
>>>> numerical noise. This can be removed with Chop
>>>>
>>>> % // Chop
>>>>
>>>> {{0.01540650723811880994, -0.0003703423761496348,
>>>> -0.1278269715054133646,
>>>>     0.0392723836429585004,
>>>>     13.10841701158527505}, {-0.00286328430472954469,
>>>> -1.985775360248369253,
>>>> \
>>>> -0.311655087881410281, -0.718055887186593765,
>>>>     21.42291506291193811}, {0.419394450195258310,
>>>> 0.01666534839706064001,
>>>>     0.528881084742103412, -0.060937975960727530,
>>>>     0.0086995512380544778}, {0.01847749159520439918,
>>>> 0.00517593237021928535,
>>>>     0.0134376958759010127, -0.0679397801010012506,
>>>> -5.21491473427110550}, \
>>>> {-0.00168718272676310080, -0.00137069152420374693, 3.05145749252874858=
2,
>>>>     0.00026536991493989622, -59.38147365750750826}}
>>>>
>>>> ?Chop
>>>>
>>>> Chop[expr] replaces approximate real numbers in expr that are close to
>>>> zero by the exact integer 0.  >>
>>>>
>>>>
>>>> Bob Hanlon
>>>>
>>>>
>>>> On Mon, Sep 3, 2012 at 3:56 AM, Andreas Talmon l'Arm=E9e
>>>> <talmon at fsm.tu-darmstadt.de> wrote:
>>>>>
>>>>> Hi,
>>>>>
>>>>> My initial conditions are the following and I am pretty sure that my
>>>>> solution consists only of a real part.
>>>>> The Solution has four eigenvalues and they are complex conjugated. Wi=
th
>>>>> all
>>>>> variables and all parameters real numbers I must be able to retrieve =
a
>>>>> real
>>>>> solution. But being a mathematica newbie,  I do not understand how to
>>>>> do
>>>>> it
>>>>> with mathematca.
>>>>>
>>>>> y''[-c] == ic0, y''[c] == ic0, y'''[-c] == ic1, y'''[c] =
== -ic1
>>>>>
>>>>> My Notebook is also ready for download at dropbox.com:
>>>>> http://dl.dropbox.com/u/4920002/DGL_4th_Order.nb
>>>>>
>>>>> Clear[sol]
>>>>>
>>>>> $Assumptions = {a \[Element] Reals, ic0 \[Element] Reals,
>>>>>      ic1 \[Element] Reals, c \[Element] Reals};
>>>>>
>>>>> sol[a_, x_] =
>>>>>    y[x] /. DSolve[{y''''[x] + a y[x] == 0, y''[-c] == ic0,
>>>>>         y''[c] == ic0, y'''[-c] == ic1, y'''[c] == -ic1},=
 y[x],
>>>>>        x][[1]] // FullSimplify
>>>>>
>>>>> Reduce[Element[sol[a, x], Reals], a, Reals]
>>>>>
>>>>>
>>>>>
>>>>> Thanks for your help,
>>>>>
>>>>> Andreas
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>> On 09/01/2012 08:27 AM, Bob Hanlon wrote:
>>>>>>
>>>>>> What are your initial conditions?
>>>>>>
>>>>>> Clear[sol]
>>>>>>
>>>>>> sol[a_, x_] = y[x] /. DSolve[
>>>>>>         {y''''[x] + a y[x] == 0,
>>>>>>          y[0] == ic0, y'[0] == ic1,
>>>>>>          y''[0] == ic2, y'''[0] == ic3},
>>>>>>         y[x], x][[1]] // FullSimplify
>>>>>>
>>>>>> (1/(2*a^(3/4)))*
>>>>>>       (Cosh[(a^(1/4)*x)/Sqrt[2]]*
>>>>>>            (2*a^(3/4)*ic0*Cos[(a^(1/4)*x)/
>>>>>>                     Sqrt[2]] + Sqrt[2]*
>>>>>>                 (Sqrt[a]*ic1 + ic3)*
>>>>>>                 Sin[(a^(1/4)*x)/Sqrt[2]]) +
>>>>>>          (Sqrt[2]*(Sqrt[a]*ic1 - ic3)*
>>>>>>                 Cos[(a^(1/4)*x)/Sqrt[2]] +
>>>>>>               2*a^(1/4)*ic2*Sin[(a^(1/4)*x)/
>>>>>>                     Sqrt[2]])*
>>>>>>            Sinh[(a^(1/4)*x)/Sqrt[2]])
>>>>>>
>>>>>> Reduce[Element[sol[a, x], Reals], a, Reals]
>>>>>>
>>>>>> a > 0
>>>>>>
>>>>>>
>>>>>> Bob Hanlon
>>>>>>
>>>>>>
>>>>>> On Fri, Aug 31, 2012 at 3:59 AM,  <"Andreas Talmon
>>>>>> l'Arm=E9e"@smc.vnet.net> wrote:
>>>>>>>
>>>>>>> Hi All
>>>>>>>
>>>>>>> Is there a way to tell mathematica to solve only for real solutions=


  • Prev by Date: Re: Exiting a Nested operation...
  • Next by Date: Re: Transforming/expanding a list
  • Previous by thread: Re: DSolve for a real function
  • Next by thread: Mathematica Prove[...] Command Possible?