Re: DSolve for a real function
- To: mathgroup at smc.vnet.net
- Subject: [mg127950] Re: DSolve for a real function
- From: Bob Hanlon <hanlonr357 at gmail.com>
- Date: Tue, 4 Sep 2012 05:48:46 -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>
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 the
> solution.
>
> I did the solution of the differential equation on paper and got the term
> for y1[x].
>
> Then I determined the constants C[5]...C[8] with the boundary conditions
> 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 for
> 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]] // FullSimplify
>>
>> ((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.051457492528748582,
>> 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. With
>>> 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. =
My
>>>>> differential equation is of the kind
>>>>>
>>>>> y''''[x]+a y[x]==0
>>>>>
>>>>> a= constant coefficient
>>>>>
>>>>> I know that I get 4 komplex eigenvalues which are complex conjungated=