Re: Re: Broken DSolve + Piecewise forcing function
- To: mathgroup at smc.vnet.net
- Subject: [mg60051] Re: [mg59642] Re: Broken DSolve + Piecewise forcing function
- From: Chris Chiasson <chris.chiasson at gmail.com>
- Date: Tue, 30 Aug 2005 04:43:08 -0400 (EDT)
- References: <ddh8t9$c91$1@smc.vnet.net> <200508151050.GAA25165@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
Maxim
Thank you for your excellent reply.
I have tried to understand your post, but your mathematics abilities
are far beyond mine.
In matrix notation pseudo code for an LTI system (which is what I had):
X'[t] =A X[t] + B U[t]
Taking the Lapace transform with all first derivatives at zero:
s X[s] = A X[s] + B U[s]
(s Identity[X[s]] - A) X[s] = B U[s]
X[s] = Inverse[[s Identity[X[s]] - A] B U[s]
X[t] = InverseLaplaceTransform[Inverse[[s Identity[X[s]] - A] B U[s]]
However, Mathematica has trouble with U[s] in the last command, even
if PiecewiseIntegrate, is used to derive U[s]. As you noted, the X[s]
solutions will have U[s] as a factor unless B is zero for one of them.
It would be easy to factor U[s] out.
From what I can gather, you have said the following:
X0[t] = InverseLaplaceTransform[Inverse[[s Identity[X[s]] - A] B]
X[t] = Integral[X0[blah] U[t-blah],{blah,0,t}]
Why is this true/possible? I see that your answer produces correct
numbers, I just hope you will illuminate why this works/is
mathematically valid.
On 8/15/05, Maxim <ab_def at prontomail.com> wrote:
> 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
>
>
--
Chris Chiasson
http://chrischiasson.com/
1 (810) 265-3161
- References:
- Re: Broken DSolve + Piecewise forcing function
- From: Maxim <ab_def@prontomail.com>
- Re: Broken DSolve + Piecewise forcing function