MathGroup Archive 2005

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

Search the Archive

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


  • Prev by Date: Re: Re: Re: export import eps | illustrator | pstoedit
  • Next by Date: Re: Re: Re: export import eps | illustrator
  • Previous by thread: Re: Broken DSolve + Piecewise forcing function
  • Next by thread: Re: Re: Broken DSolve + Piecewise forcing function