Re: Re: Broken DSolve + Piecewise forcing function
- To: mathgroup at smc.vnet.net
- Subject: [mg60065] Re: Re: Broken DSolve + Piecewise forcing function
- From: Maxim <ab_def at prontomail.com>
- Date: Wed, 31 Aug 2005 02:19:27 -0400 (EDT)
- References: <ddh8t9$c91$1@smc.vnet.net> <200508151050.GAA25165@smc.vnet.net> <df16rs$5ql$1@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
This is simply the convolution theorem: the Laplace transform maps the
convolution to the product. So IL[f*g] == convolution[F, G] (I used
capital letters for the originals and lower-case letters for their images
under the direct transform, and IL denotes the inverse Laplace transform).
Also we should remember that the inverse transform returns a function
which is zero for t<0, so the convolution can be reduced to the integral
over a finite range.
Refining your notation a bit, we're applying the Laplace transform to the
system
{X'[t] == A.X[t] + F[t], X[0] == 0}:
L[X'[t] - (A.X[t] + F[t])] == (s*Id - A).x[s] - f[s] == 0,
where Id is the identity matrix of the corresponding order. Then
x[s] == Inverse[s*Id - A].f[s],
X[t] == IL[x[s]] == convolution[IL[Inverse[s*Id - A]], F[t]].
Further, by analogy with the identity IL[1/(s - a)] == E^(a*t), we have
IL[Inverse[s*Id - A]] == MatrixExp[t*A];
X[t] == convolution[MatrixExp[t*A], F[t]]
== Integrate[MatrixExp[tau*A].F[t - tau],
{tau, -Infinity, Infinity}]
== Integrate[MatrixExp[tau*A].F[t - tau], {tau, 0, t}],
and in fact we need PiecewiseIntegrate only to evaluate this last step:
{exactP[2][t_], exactQ[4][t_]} =
PiecewiseIntegrate[#, {tau, 0, t}, Assumptions -> t > 0]& /@
(MatrixExp[{{0, -10000}, {2500, -500}}*tau].
{Piecewise[{{50, 2/5 > FractionalPart[20000*(t - tau)]}}], 0})
At the risk of obfuscating the things even more, I'll show another
possible approach to this problem. F[t] is periodic, so we can expand it
into the Fourier series, solve the ODEs for each of the terms, and take
the sum of the solutions. If we choose the basis
phi[k_, t_] = E^(2*Pi*I*k*t/T) /. T -> 1/20000;
then the kth coefficient in the expansion of F[t] is
c[0] = 20;
c[k_] = -25*I*(1 - E^(-4/5*I*Pi*k))/(k*Pi);
The solution for the kth term:
sol[k_] := First@ DSolve[
{P[2]'[t] == c[k]*phi[k, t] - 10000*Q[4][t],
Q[4]'[t] == 2500*P[2][t] - 500*Q[4][t],
P[2][0] == Q[4][0] == 0},
{P[2], Q[4]}, t]
Then we should be able to construct Q[4] like this:
exactQ[4][t_] = (Q[4][t] /. sol[0]) +
Sum[Q[4][t] /. sol[k], {k, -Infinity, -1}] +
Sum[Q[4][t] /. sol[k], {k, Infinity}]
But here we run into a lot of various errors in Sum and Series (when
HypergeometricPFQ is involved). It is possible to get the exact solution
this way, but it will take a lot of extra tweaking.
Maxim Rytin
m.r at inbox.ru
On Tue, 30 Aug 2005 08:54:20 +0000 (UTC), Chris Chiasson
<chris.chiasson at gmail.com> wrote:
> 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
>>
>>
>
>
- References:
- Re: Broken DSolve + Piecewise forcing function
- From: Maxim <ab_def@prontomail.com>
- Re: Broken DSolve + Piecewise forcing function