Re: DSolve for a real function
- To: mathgroup at smc.vnet.net
- Subject: [mg127942] Re: DSolve for a real function
- From: Andreas Talmon l'Armée at smc.vnet.net
- Date: Tue, 4 Sep 2012 05:46:06 -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> <20120901062746.B3D7F687B@smc.vnet.net> <504462A6.5070404@fsm.tu-darmstadt.de> <CAEtRDSfhfr=mca7_3XsfGrF_5MbAuvAF5f8oVo7y+M9=_C68Ow@mail.gmail.com> <5044BD10.5070102@fsm.tu-darmstadt.de> <CAEtRDSdHoBwvVi3ewxnnnQD_rrT_+PFDKi1sHnWct3L9eJQt0g@mail.gmail.com>
Hi, I think this is not a numerical problem. Because if I look at the outcome of the symbolic computations for y3[x] (Simplify) and y4[x] (FullSimplify) there is one solution without imaginary part (Simplify) and one with imaginary part (FullSimplify) . http://dl.dropbox.com/u/4920002/DGL_4th_Order_with_own_solution.nb I am feeling a bit uncomfortable with these outcomes because I think the solution y3[x] (Simplify) is the better one even if the real part of both solutions is the same. But I think I just have to accept it. Is there a way to suppress the calculation with complex numbers. Thanks for your time and your help Andreas Talmon 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ée > <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ée >>> <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. >>>>>> But y[x] is a real function. >>>>>> Solving this equation with DSolve always gets a complex function y[x]. >>>>>> >>>>>> Any Ideas. >>>>>> >>>>>> Thanks, Andreas >>>>>>
- References:
- Re: DSolve for a real function
- From: Bob Hanlon <hanlonr357@gmail.com>
- Re: DSolve for a real function