Re: Integrating a Piecewise function four times but different results
- To: mathgroup at smc.vnet.net
- Subject: [mg121791] Re: Integrating a Piecewise function four times but different results
- From: Peter Pein <petsie at dordos.net>
- Date: Sun, 2 Oct 2011 02:37:08 -0400 (EDT)
- Delivered-to: l-mathgroup@mail-archive0.wolfram.com
- References: <j66e9t$hch$1@smc.vnet.net>
Am 01.10.2011 09:08, schrieb Rob Y. H. Chai:
> Hello everyone,
>
>
>
> I tried to integrate a piecewise function 4 times but got different results
> for two supposedly identical cases. The integration was done in the Mathematica
> version 8.0.1.0. I am interested in what I did wrong, so any comments will
> be appreciated. My code is below.
>
>
>
> Sincerely
>
>
>
> Rob Y. H. Chai
>
> UC Davis
>
>
>
>
>
> (* Introduction p1[x] as piecewise funtion *)
>
> p1[x_] = Piecewise[{{2*wmax*x/L, x <= L/2}, {2*wmax*(1 - x/L), x > L/2}}];
>
>
>
>
>
> (* p2[x] is a specialized case of p1[x] with wmax = 1 and L = 1 *)
>
> p2[x_] = p1[x] /. wmax -> 1 /. L -> 1;
>
>
>
>
>
> (* now integrate p1[x] 4 times and introduce integration constants *)
>
> soln1[x_] = Integrate[Integrate[Integrate[Integrate[p1[x], x], x], x], x] +
> a*x^3 + b*x^2 + c*x + d ;
>
>
>
> (* solve for integration constants in the solution of p1[x] *)
>
> eq1a = soln1[0] == 0;
>
> eq1b = (D[soln1[x], {x, 2}] /. x -> 0) == 0;
>
> eq1c = soln1[L] == 0;
>
> eq1d = ( D[soln1[x], {x, 2}] /. x -> L) == 0;
>
> soln1A[x_] = soln1[x] /. Simplify[Solve[{eq1a, eq1b, eq1c, eq1d}, {a, b, c,
> d}], L > 0][[1]];
>
> (* check case 1 solution at x = L/2 *)
>
> soln1A[L/2] /. wmax -> 1 /. L -> 1
>
> 17/1920
>
>
>
>
>
> (* similarly integrate p2[x] 4 times and introduce integration constants *)
>
> soln2[x_] = Integrate[Integrate[Integrate[Integrate[p2[x], x], x], x], x] +
> r*x^3 + s*x^2 + u*x + v;
>
> (* solve for integration constants in the solution of p2[x] *)
>
> eq2a = soln2[0] == 0;
>
> eq2b = (D[soln2[x], {x, 2}] /. x -> 0) == 0;
>
> eq2c = soln2[1] == 0;
>
> eq2d = ( D[soln2[x], {x, 2}] /. x -> 1) == 0;
>
> soln2A[x_] = soln2[x] /. Simplify[Solve[{eq2a, eq2b, eq2c, eq2d}, {r, s, u,
> v}]][[1]];
>
> soln2A[1/2] (* case 2 solution at x = 1/2 *)
>
> 1/120
>
>
>
> Plot[{soln1A[x] /. wmax -> 1 /. L -> 1, soln2A[x]}, {x, 0, 1}]
>
>
>
> Question - why doesn't case 1 solution converge to case 2 solution when case
> 2 is simply a special case of case 1?
Case 1 has got a discontinuity at L/2 (already after the first
integration). I would prefer to iterate (Integrate w.r.t. x from 0 to z,
add a constant and replace z by x) four times:
assum = {L > 0 < wmax, Element[z, Reals]};
{p0[x_], p2[x_]} =
FoldList[FullSimplify[PiecewiseExpand[
c[#2] + Integrate[#1, {x, 0, z}, Assumptions -> assum] /. z -> x],
assum] & ,
Piecewise[{{(2*wmax*x)/L, x <= L/2}, {2*wmax*(1 - x/L), x > L/2}}],
Range[4]]
[[{-1, -3}]];
will give us the candidate function p0 and its second derivative.
the conditions are:
eqns = Simplify[Flatten[(#1 /@ {0, L} & ) /@ {p0, p2}] == 0, assum];
so we get:
In[4]:= pfunc[x_] = FullSimplify[p0[x] /. First[
Solve[eqns, Union[Cases[eqns, _c, Infinity]]]],
assum]
Out[4]= Piecewise[{
{(wmax*x*(5*L^2 - 4*x^2)^2)/(960*L), x == 0 || L >= 2*x}},
(wmax*(L - x)*(L^2 + 8*L*x - 4*x^2)^2)/(960*L)]
and the lim_{x->L/2}{pfunc(x)} is the same for both directions:
In[5]:= (Limit[pfunc[x], x -> L/2,
Assumptions -> assum, Direction -> #1] & ) /@ {1, -1}
Out[5]= {(L^4*wmax)/120, (L^4*wmax)/120}
Hope that helps,
Peter