Services & ResourcesWolfram ForumsMathGroup Archive

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

Re: Broken DSolve + Piecewise forcing function

• To: mathgroup at smc.vnet.net
• Subject: [mg59636] Re: [mg59566] Broken DSolve + Piecewise forcing function
• From: Devendra Kapadia <dkapadia at wolfram.com>
• Date: Mon, 15 Aug 2005 06:50:28 -0400 (EDT)
• References: <200508120409.AAA11249@smc.vnet.net>
• Sender: owner-wri-mathgroup at wolfram.com

On Fri, 12 Aug 2005, Chris Chiasson 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,
>
Hello Chris,

Thank you for reporting the unexpected answer given by DSolve for
the system of equations with a Piecewise forcing function.

The problem appears to be caused by the fact that DSolve does not
have any information about the range of 't' and hence the condition
based on FractionalPart is not resolved into simple linear inequalities.
The solution seems to be formally correct, but is not very useful.

One method of overcoming this problem is to use $Assumptions to specify
the the range of 't'.

Here is a simple example, based on your notebook, which illustrates
the use of $Assumptions in this situation.

We first compare the results from DSolve and NDSolve over a small 
time interval without specifying $Assumptions.

=====================================================

In[1]:= $Version

Out[1]= 5.2 for Linux (June 27, 2005)

In[2]:= eqns = {Derivative[1][P[2]][t] == Piecewise[{{50, 2/5 >
                    FractionalPart[20000*t]}}, 0] - 10000*Q[4][t],
                 Derivative[1][Q[4]][t] == 2500*P[2][t] - 500*Q[4][t],
                 P[2][0] == 0,
                 Q[4][0] == 0};

In[3]:= vars = {P[2], Q[4]};

In[4]:= sol1 = NDSolve[eqns, vars, {t, 0, 7/20000}, MaxSteps -> 10000000,
                    MaxStepSize -> 1/60000, WorkingPrecision -> 20];

In[5]:= t1 = N[Table[Q[4][t] /. sol1[[1]], {t, 0, 5/20000, 1/40000}]]

Out[5]= {0., 0.0000372886, 0.000098287, 0.000194292, 0.000310843,
          0.000458414, 0.000621794, 0.000810791, 0.00100962, 0.00122761,
          0.0014486}

In[6]:= Timing[sol2 = DSolve[eqns, vars, t]; ]

Out[6]= {0.4 Second, Null}

In[7]:= t2 = N[N[Table[Q[4][t] /. sol2[[1]], {t, 0, 5/20000, 1/40000}], 20]]

Out[7]= {0., -0.00496115, 0.000154152, -0.00465683, 0.000602054,
          -0.00407408, 0.00130895, -0.00325554, 0.00222504, -0.00225731,
          0.00328891}

===============================================================

Clearly, the values given by NDSolve and DSolve (Out[5] and Out[7]
above) do not agree.

We now repeat the above calculation after specifying $Assumptions.

=============================================================

In[1]:= $Version

Out[1]= 5.2 for Linux (June 27, 2005)

In[2]:= $Assumptions = { 0 < t < 6/20000}

                    3
Out[2]= {0 < t < -----}
                  10000

In[3]:= eqns = {Derivative[1][P[2]][t] == Piecewise[{{50, 2/5 >
                   FractionalPart[20000*t]}}, 0] - 10000*Q[4][t],
                 Derivative[1][Q[4]][t] == 2500*P[2][t] - 500*Q[4][t],
                 P[2][0] == 0,
                 Q[4][0] == 0};

In[4]:= vars = {P[2], Q[4]};

In[5]:= sol1 = NDSolve[eqns, vars, {t, 0, 5/20000}, MaxSteps -> 10000000,
                   MaxStepSize -> 1/60000, WorkingPrecision -> 20];

In[6]:= t1 = N[Table[Q[4][t] /. sol1[[1]], {t, 0, 5/20000, 1/40000}]]

Out[6]= {0., 0.0000372886, 0.0000982869, 0.000194292, 0.000310843,
          0.000458414, 0.000621794, 0.000810791, 0.00100962, 0.00122761,
          0.0014486}

In[7]:= Timing[sol2 = DSolve[eqns, vars, t]; ]

Out[7]= {54.11 Second, Null}

In[8]:= t2 = N[N[Table[Q[4][t] /. sol2[[1]], {t, 0, 5/20000, 1/40000}], 20]]

Out[8]= {0., 0.0000372885, 0.0000982869, 0.000194292, 0.000310843,
          0.000458414, 0.000621794, 0.000810791, 0.00100962, 0.00122761,
          0.0014486}

=================================================================

Notice that DSolve takes much longer to return the solution (see Out[7]
above). However, the results from DSolve and NDSolve now agree quite
closely.

You could try comparing the graphs from DSolve and NDSolve in your 
notebook using $Assumptions = {0 < t < 10/20000} (the graphs should
then be plotted within this range).

For the entire interval {0 < t < 1/2}, the number of cases in Piecewise 
will require setting $MaxPiecewiseCases to some high value (such as
100000). Also, the calculation is unlikely to finish in any reasonable
amount of time. Considering the complexity of the solution in this case,
NDSolve appears to be the correct function for handling this problem.

Sorry for the inconvenience caused to you and for the delay in responding
to your email.

Sincerely,

Devendra Kapadia.
Wolfram Research, Inc.

• References:
  • Broken DSolve + Piecewise forcing function
    • From: Chris Chiasson <chris.chiasson@gmail.com>