Re: Broken DSolve + Piecewise forcing function
- To: mathgroup at smc.vnet.net
- Subject: [mg59642] Re: Broken DSolve + Piecewise forcing function
- From: Maxim <ab_def at prontomail.com>
- Date: Mon, 15 Aug 2005 06:50:33 -0400 (EDT)
- References: <ddh8t9$c91$1@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
On Fri, 12 Aug 2005 04:34:49 +0000 (UTC), Chris Chiasson <chris.chiasson at gmail.com> wrote: > Dear MathGroup, > Somehow, I have managed to feed DSolve a system of equations to which > it seems to produce an unphysical answer. The problem seems to > originate from the forcing function, which is Piecewise defined. I > have correctly solved the equations with NDSolve. I have also > correctly solved the equations using a Fourier series approximation to > the forcing function and DSolve. All that remains is to figure out > what is going wrong with the Piecewise forcing function and the DSolve > command. > The notebook gives the full problem statement and some explanation on > the origins of the equation system. > http://chrischiasson.com/MathGroup/BrokenDSolvePiecewise.nb > Please help me figure out what is wrong. > Thank you for your time, The system of equations you're solving is {P[2]'[t] == Piecewise[{{50, 2/5 > FractionalPart[20000*t]}}] - 10000*Q[4][t], Q[4]'[t] == 2500*P[2][t] - 500*Q[4][t], P[2][0] == Q[4][0] == 0} The solution given by DSolve isn't continuous. This is what we would get if we solved the system up to t = 1/20000 and then picked a different initial condition P[2][1/20000] == P20 and used this initial condition when solving the system on the next interval. The resulting function satisfies the differential equations on each of the intervals, but I suppose this isn't what you wanted. It's possible to get the continuous solution by using the Laplace transform. We can find the image of the Piecewise term by using PiecewiseIntegrate ( http://library.wolfram.com/infocenter/MathSource/5117/ ): F[t_] = Piecewise[{{50, 2/5 > FractionalPart[20000*t]}}] f[s_] = PiecewiseIntegrate[F[t]*E^(-s*t), {t, 0, Infinity}] Then the transform of P[2]'[t] is s*p[2][s] - P[2][0] == s*p[2][s], and we obtain an algebraic system: sol = First@ Solve[ {s*p[2][s] == f[s] - 10000*q[4][s], s*q[4][s] == 2500*p[2][s] - 500*q[4][s]}, {p[2][s], q[4][s]}] We can't find the inverse transform directly, but we notice that the solution contains f[s] as a factor: p[2][s] == f[s]*p0[s], q[4][s] == f[s]*q0[s]. So the original function is the convolution of F[t] with the originals of p0 and q0: {P0[t_], Q0[t_]} = InverseLaplaceTransform[ {p[2][s], q[4][s]}/f[s] /. sol // Together, s, t] The convolution can be evaluated with PiecewiseIntegrate too: {exactP[2][t_], exactQ[4][t_]} = PiecewiseIntegrate[#[tau]*F[t - tau], {tau, 0, t}]& /@ {P0, Q0} This will take about 6 minutes to evaluate and will return a huge symbolic result, which is numerically the same as the NDSolve solution: Plot[{exactP[2][t], exactQ[4][t]} // Evaluate, {t, 0, .035}, PlotRange -> All] In some situations it might be necessary to use Re because exactP and exactQ contain complex quantities even though they are real-valued. A couple of additional comments. Evaluate inside Plot is not necessary but it is a more efficient way. Without Evaluate only the unevaluated expression exactP[2][t] will be placed inside Compile, and the actual huge expression that exactP[2][t] evaluates to will remain uncompiled. exactP[2][t] and exactQ[4][t] are complex enough to make the performance difference quite noticeable. The second observation is that Piecewise formatting rules are not perfect. Here are two expressions: Piecewise[{{1, x > a - g'[x]}}] Piecewise[{{1, x > a}}] - g'[x] If they're displayed in StandardForm or TraditionalForm it is next to impossible to work out whether g'[x] is included in the condition or it is a separate term. A quick workaround is to add a formatting rule to always wrap Piecewise in parentheses: Unprotect[Piecewise]; $PWparF = True; PrependTo[FormatValues[Piecewise], HoldPattern[MakeBoxes[expr_Piecewise, form_] /; $PWparF] :> Block[{$PWparF = False}, RowBox[{"(", MakeBoxes[expr, form], ")"}] ]]; Maxim Rytin m.r at inbox.ru
- Follow-Ups:
- Re: Re: Broken DSolve + Piecewise forcing function
- From: Chris Chiasson <chris.chiasson@gmail.com>
- Re: Re: Broken DSolve + Piecewise forcing function